Bias Correction of Satellite-Based Rainfall Estimates for Modeling Flash Floods in Semi-Arid regions : Application to Karpuz River , Turkey

Bias Correction of Satellite-Based Rainfall Estimates for Modeling Flash Floods in Semi-Arid regions: Application to Karpuz River, Turkey Mohamed Saber 1, 2, , Koray K. Yilmaz Geology Department, Faculty of Science, Assiut University, Assiut 71516, Egypt 5 Department of Geological Engineering, Middle East Technical University, 06800, Ankara, Turkey Water Resources Research Center, DPRI, Kyoto University, Goka-sho, Uji City, Kyoto 611-0011, Japan; Correspondence to: Mohamed Saber (mohamedmd.saber.3u@kyoto-u.ac.jp)


Introduction
Spatio-temporal variation of rainfall is important to understand the hydrological and climatic characteristics of watersheds, as well as for planning effective water management and hazard mitigation strategies for water-related disasters such as flash floods and droughts.The influence of rainfall representation on the modelling of the hydrologic response is expected to depend on complex interactions between the rainfall space-time variability, the variability of the catchment soil and landscape properties, and the spatial scale (i.e.catchment area) of the problem (Obled et al., 1994;Woods and Sivapalan, 1999;Bell and Moore, 2000;Smith et al., 2004).River hydrograph forecasts are highly dependent upon the input rainfall data used.
Streamflow in arid and semi-arid regions is characterized by rapid response to intense rainfall events.Such events frequently have a high degree of spatial variability.Coupled with poorly gauged rainfall data this situation sets a fundamental limit on the capacity of the rainfall-runoff models to reproduce the observed flow (Wheater et al., 2008), thus hampering prediction efforts.It has been widely stated that the major limitation of the development of arid-zone hydrology is the lack of high quality observations (McMahon and Greene, 1979;Pilgrim et al., 1988).
Accurate temporal knowledge of global precipitation is definitely important for understanding the multi-scale interactions among weather, climate and ecological systems, as well as for improving the ability to manage the available water resources and predict high-impact weather events including hurricanes, floods, droughts and landslides (Hou et al., 2008).Rain gauge observations yield relatively accurate point measurements of precipitation but are not well distributed and not available over most oceanic and unpopulated land areas (Xie and Arkin, 1996;Petty and Krajewski, 1996).In particular, the arid and semiarid regions are suffering from poor rain gauge network coverage and the lack of continuous observations, especially at the hourly timescale, which will be consequently a challenge for the flash floods mitigation efforts.Therefore, algorithms developed for bias correction of the satellite-based rainfall data are desperately needed in such regions.
Recent decades are marked by increasing flash flood hazards in many regions over the world due to increase in the frequency of flash floods as well as rapid urbanization.The term 'flash flood' identifies a rapid hydrological response, with water levels reaching a peak within less than one hour to a few hours after the onset of the generating rain event (Creutin et al., 2013;Collier, 2007;Younis et al., 2008).The most important challenge to model implementation and calibration to simulate flash floods in semi-arid and arid regions is the lack of necessary observational networks of both rainfall and discharge.
Responsible factors for the short duration of the flash flood include intense rains that persist on an area for a few hours, steep slope, impermeable surfaces, and sudden release of impounded water over small basins (Georgakakos, 1986).
Flood occurrences are complex since they depend on interactions between many geological and morphological characteristics of the basins, including rock types, elevation, slope, sediment transport, and flood plain area.Moreover, hydrological phenomena, such as rainfall, runoff, evaporation, and surface and groundwater storage (Farquharson et al., 1992;Flerchinger and Cooley, 2000;Nouh, 2006) affect floods.According to (Few et al., 2004), each flood acquires some particular and inherent characteristics of the occurrence locality, such as flow velocity and height, duration, and rate of water-level rise.
In many countries and regions of the world, flash floods are the most costly natural hazards in terms of both loss of human lives and material damage (Fattorelli et al., 1999;Creutin et al., 2013).In particular, arid and semi-arid regions have become more vulnerable to flash floods than before possibly due to the global climate change and rapid population growth.The main obstacle to study flash floods is clearly the lack of reliable observations in most of the flash flood prone basins.Furthermore, the danger also comes from the rarity of the phenomenon, which demands new observation strategies, as well as new forecasting methodologies. . Hazards Earth Syst. Sci. Discuss., doi:10.5194/nhess-2016-339, 2016 Manuscript under review for journal Nat.Hazards Earth Syst.Sci.Published: 11 November 2016 c Author(s) 2016.CC-BY 3.0 License.

Nat
Various problems associated with forecasting flash floods caused by convective storms over semi-arid basins have been studied by (Michaud and Sorooshian, 1994).Rapidly increasing availability of good quality weather radar observations is greatly expanding our ability to measure and monitor rainfall distribution at the space and time scales which characterize the flashflood events (Borga et al., 2007).Moreover, some hydrological approaches and understanding of runoff characteristics in arid environment have been developed by (Saber et al., 2010b;Saber et al., 2013), and a method for estimating flash flood peak discharge, hydrograph, and volume has been presented by (Koutroulis and Tsanis, 2010).Hence an integrated and comprehensive research regarding flash flood modeling and forecasting approaches as well as mitigation strategies are desperately needed in semi-arid and arid regions.
A suite of models is available to represent rainfall-runoff relationships, but they have limitations in the hydrologic parameters that are used to describe the rainfall-runoff process in semi-arid and arid systems (Wheater et al., 1993).It has been widely stated that the major limitation of the development of arid-zone hydrology is the lack of high quality observations (McMahon and Greene, 1979;Pilgrim et al., 1988).Thus, implementation of a hydrological model driven by locally corrected satellitebased rainfall estimates could be useful in overcoming majority of the problems in simulating flash floods, thus potentially serve as a tool for hazard mitigation and sustainable development of the target basin.
A comparative study has been done (Saber, 2010;Saber et al., 2013) between GSMaP and the Global Precipitation Climatology Center (GPCC) in the arid and semi-arid regions over the globe to characterize the GSMaP product bias as compared to observed GPCC product.Additional analysis using the local raingauge network at a semi-arid region will be considered in this research for further validation of GSMaP product.Consequently, the forecasting models driven by the bias-corrected satellitebased rainfall datasets are expected to be more powerful and reliable.This study aims to compare GSMaP product with the gauge-based precipitation estimates in Karpuz River located in Antalya, Turkey in an effort to devise a correction methodology for the GSMaP product to drive a hydrological model for flash-flood simulation.Due to rapid occurrence of flash floods at sub-daily time scales generally hourly spatio-temporal rainfall data is used for flash flood simulation studies.Thus, satellitebased rainfall data at the hourly timescale with continuous availability in time provides an alternative to ground-based Nat.Hazards Earth Syst.Sci. Discuss., doi:10.5194/nhess-2016-339, 2016 Manuscript under review for journal Nat.Hazards Earth Syst.Sci.Published: 11 November 2016 c Author(s) 2016.CC-BY 3.0 License.observations for flash flood simulation studies.The paper consists of two main parts.First, data comparison and the procedure for correction of Satellite-based rainfall data (GSMaP) is introduced.Second, flash flood modeling through a hydrological model with calibrated and validated parameters for Karpuz River -a semi-arid basin in Antalya, Turkey is provided.

Study Area and Datasets
The study area is the Karpuz River Basin located to the west of the city of Antalya situated in the Mediterranean region of 5 Turkey.The study area lies between 30.50E-32.50E longitude bands and 36.00N-37.50NLatitude bands with a total area of about 1920 km2, where the land data pixels were only considered for the analysis (Fig. 1).

Comparison of GSMaP product with Rain gauges
The main objective of this comparison was to validate GSMaP data products with the available rain gauge data in the vicinity of the Karpuz River Basin (Antalya, Turkey), in order to contribute to the enhancement of hydrological forecasting of flash floods in semi-arid regions.The time period for the evaluation was selected as the years 2007 through 2013 based on the availability of the rain gauges.Statistical analysis was performed for the GSMaP data in order to estimate the data bias and to examine its feasibility to be used for the flash floods simulations.The gauge-based rainfall data were interpolated to a 0.1° x 0.1° resolution grid to be consistent with the GSMaP product resolution using the automated Thiessen polygons (Han and Bray, 2006), based on the distance formula: where D is the distance between the rain gauge and the target pixel (x, y), i refers to the number of station and j for the pixel number.In the procedure, the pixels are assigned the same rainfall rate with the nearest station.
Monthly rainfall averages over the target area were calculated both for hourly GSMaP data, and daily rain gauges with the same spatial resolution.The analysis for comparison consisted of various scenarios including different time series and different rainfall thresholds, areal average and pixel to pixel comparisons.Four thresholds (0 mm, 1 mm, 5 mm, and 10 mm) were examined in order to conduct the comparison between rain gauges and GSMaP rainfalls data properly.Statistical parameters considered in the comparison include: (i) Bias (Eq.2) the ratio of spatially averaged monthly rainfall obtained from GSMaP and rain gauges, (ii) Percent bias (PBIAS; Eq. 3) measures the average tendency of the satellite data to the rain gauge estimates based rainfall and gauge-based rainfall.R2 ranges between 0 and 1, with higher values indicating less error variance, and typically values greater than 0.5 are considered acceptable (Santhi et al., 2001;Van Liew et al., 2003).
Where    is the gauge precipitation values at time i,    is the Satellite precipitation values The results of the statistical analysis, as summarized in Table 3, show that monthly rainfall estimates derived from GSMaP data are linearly correlated with the rain gauges but significant underestimation of bias is evident.It can also be seen that the bias magnitude is slightly changing as a function of season (monthly complete time series, April-September, or October to March), as well as the rainfall threshold (0 mm, 1 mm, 5 mm, and 10 mm) as shown in Table .Scatter plots and hyetographs provided in Figs 2, 3, 4, and 5 show a reasonable and acceptable linear correlation but with an underestimated bias in GSMaP data in most of the discussed cases, especially for the September-April time period.The total average bias factors for the monthly data comparison indicate underestimation for 1mm threshold case with values 0.78, 0.50, and 0.98 for the three time periods, namely monthly complete time series, October-March and April-September, respectively (Table 3).We found that the best reasonable correlation recorded in case of 1mm threshold is for the time period of April-September.Also, the best PBias value also corresponds to the same threshold and time period, implying that the underestimation bias is higher in the rainy season.The statistical analysis show that both 1mm and 5 mm thresholds are better than the 0mm and 10 mm thresholds.
This situation might be due to the fact that 0mm threshold include all the pixels in the analysis whereas 10 mm threshold excluding most of the pixels.The GSMaP data was corrected using the bias ratio of 1 mm threshold case, which in turn used to drive the hydrological model for effective flash floods simulation.Focusing on the total monthly average precipitation for the 2007-2013 time period (Table 4), the analysis show that the bias and PBias are slightly changed for each precipitation threshold and time period compared to the previous case (Table 3).
Correlation coefficient and NSE exhibit a significant enhancement for total rainfall average estimates derived from GSMaP data.Scatter plots and hyetographs (Figs 6, 7, 8, and 9) show a reasonable and acceptable linear correlation but with a 5 significant underestimation bias in GSMaP data in most of the discussed cases, especially for the September-April time period.

Point vs. grid comparison (rain gauges with the corresponding satellite pixels)
Point-scale rainfall measurements at Antalya meteorological station (Lat.36.92,Lon.30.8) was compared with the overlying 5 GSMaP-based grid rainfall estimates.The results of comparison show very good correlation between GSMaP data and rain gauge data (Table 5, and Figs 12, 13).Different thresholds were also considered including 0 mm and 5 mm.The different scenarios that we have discussed show that GSMaP and rain gauges show good linear correlation but significant underestimation bias persist in GSMaP estimates which depend on the seasons and the selected rainfall thresholds.Therefore, in order to select the appropriate threshold bias results for the GSMaP data correction, we calculated the number of daily rainfall occurrences above each threshold at different locations (Fig. 14).This analysis focusing on the number of daily occurrences of above-threshold rainfall for GSMaP and gauges for the whole basins and selected rain gauge stations indicate that threshold value of 1 mm is the best choice for the bias correction between both GSMaP and rain gauges at all stations.
The counting of 10 mm thresholds is not a good option as it excludes most of the rainy days and most of the pixels.Also, 0 mm threshold is considering most of the pixels and days but the problem is that the difference between the two products is too large -the number of days in GSMaP is higher than those of rain gauges about 20% for the whole daily time series at all stations.Therefore, the best choice based on this analysis as well as the previous statistical analysis is 1mm threshold because it is exhibiting reasonable daily rainfall occurrence and also the difference between days of GSMaP and rain gauges is not remarkable -approximately about 2% difference in all cases.Another interesting issue was found in the difference between the three stations.In case of 5mm and 10 mm thresholds, GSMaP shows underestimated numbers of days at the three stations, but at Ibradi rain gauge station where the elevation is 1036 m, the difference is more significant compared with the other two stations at low elevations.For instance, at Ibradi rain gauge, the difference is 141 days, but at Antalya rain gauge the difference is only 8 days.This result implies that the GSMaP product has a more tendency to underestimate rainfall at high elevations compared to low elevations.This may be due to the snow effect, as stated by (Derin and Yilmaz, 2014) that there are major challenges to satellite-based precipitation estimation algorithms over complex topography such as those related to orographic precipitation and precipitation estimation over cold surfaces.For instance, their analysis at the daily time scale revealed that the CMORPH product suffered from daily precipitation detection problems specifically in the cold season and the windward region, this might be due to the surface snow and ice screening procedure embedded in the algorithm (Joyce et al., 2004) .The important finding here is that for bias correction of the satellite based data, it is recommended to select the appropriate threshold for the bias correction.The second finding is that GSMaP data is showing an underestimation bias in all the discussed 5 cases, however it is showing good linear correlations with rain gauge data.The third issue is that GSMaP data shows an elevation dependent bias.The bias factors estimated from this analysis were used for the correction of GSMaP data in order to use for the flash floods simulation in the second part of the study.Recent studies (Özgüler, 2003) indicated that the frequency floods and flash floods has dramatically increased in the region (Fig. 17).In Turkey, flooding is the second important natural hazard after the earthquakes, with 22 floods and 19 deaths per year on average (Özgüler, 2003) with ever increasing economic loss (Fig. 17).Turkey is the fourth country in terms of the greatest losses from flash floods after Italy, France, and Romania (Llasat et al., 2010) .

Model Description
Hydrological River Basin Environmental Assessment Model (Hydro-BEAM) was selected to simulate the flash flood events in the study basin.Hydro-BEAM (Fig. 18) is a physically-based distributed hydrological model originally developed (Kojiri et al., 1998).Hydro-BEAM model has been adopted and implemented in the arid regions by (Saber et al., 2010a;Saber, 2010;Saber et al., 2013) .The components of the Hydro-BEAM model utilized in this study include a GIS interface for data input and visualization, surface runoff and stream routing modeling based on the kinematic wave approximation, the initial and transmission losses modeling is estimated via the Curve Number approach (SCS, 1997) and Walter's equation (Walters, 1990) respectively, and groundwater modeling based on the linear storage model.However, understanding of hydrological process resulting in flash floods is hampered by the observational difficulties.In this study, various remotely sensed datasets was used as input to the hydrological model.Hydro-BEAM is a distributed model consisting of horizontal spatial discretization the scale of which could be adjusted based on the basin size.Vertically, each pixel is represented by a combination of one surface layer and three subsurface layers (Fig. 18).The surface and subsurface layers are named as A, B, C and D. A-Layer and the river channel is governed by the kinematic wave model for the overland flow calculation and the layers C and D (subsurface layers) are modeled using the linear storage model.Climate data input is one of the most important factors for hydrological models.Therefore, spatio-temporal distribution of rainfall and evaporation are needed.Hourly satellite-based GSMaP rainfall data bias corrected following the methodology described in Section 2.1was used as input to the model.In the procedure, calculated bias for the daily timescale was assumed to be also valid at the hourly time scale.
Thornthwaite method was adopted to calculate daily mean potential evapotranspiration for each grid because of its simplicity and data availability.This method only requires mean air temperature and duration of possible sunshine for each grid as meteorological input.Thorntwaite method was applied with the following equations (6-9): = .  (9 where Ep is the daily-averaged potential evapotranspiration in i-th month (mmd-1), Do is the possible sunshine duration (h/12h), J is the heat index, Ti is the a monthly-averaged temperature in i-th month, Ea is the daily-averaged actual evapotranspiration in i-th month (mmd-1), and M is a evapotranspiration coefficient (reduction coefficient, vapor effective parameter).
Next, main Hydro-BEAM model calculates the streamflow discharge using the kinematic wave model and linear storage model.

Flash floods Simulation at Karpuz River
The streamflow data for the Karpuz River was available only for the period Oct. 2007-Sept. 2013, with a gap between Feb. 2011-Sept. 2012.Therefore, simulation time period was selected based on data availability.The advantage of using satellite- based rainfall data compared to rain gauges include better spatial and temporal (hourly) resolution of the former without gaps which would be a better choice for hydrologic modeling studies including flash floods simulation.Moreover, we utilized the bias correction procedure discussed in Section 2.2.for GSMaP product, before driving the HydroBEAM model for simulation of flash floods events within the time period from 2007-2013.

Model parameter calibration and validation
The model parameters (Table 6) were calibrated using the flash flood events from Oct. to Dec. 2007 and then validated by the flash flood events from Oct to Dec. 2009.The performance of the model during calibration and validation process were assessed with statistical measures including correlation coefficient (R), standard deviation, bias, and the Kling-Gupta efficiency (KGE).The Kling-Gupta efficiency (KGE), is an alternative model performance criterion developed by (Gupta et al., 2009) and simply combines the above measures into a single criterion as follows: = √ ( − 1) 2 + (∝ −1) 2 + ( − 1) 2 ( 10) Where ED is the Euclidian distance from the ideal point:  is the ratio between the mean simulated and mean observed flows, i.e.  represents the bias; r: is linear correlation coefficient between simulated and observed flows, ∝ is the ratio between standard deviations of simulated and observed flows ( a measure of the relative variability of flows).
The results indicated an acceptable performance (Table 7) of the hydrological model with good results for bias (0.91 and 0.

Spatial distribution maps of bias-corrected GSMaP rainfall and simulated discharge
Spatial variability of simulated discharge (Fig. 21) were investigated to identify the regions that are most prone to the flash flood events and hence to enable planning appropriate mitigations strategies.Moreover, GSMaP-based rainfall distribution maps were also visualized to exhibit the variability of rainfall over the target basin (Fig. 22).These maps represent four snapshots during the onset of a flood event with increase in discharge rate at the downstream outlet from 13 m3/sec to 80 5 m3/sec within 10 hours (see these snapshots on hydrograph in Fig. 20a).This confirms the reality of the short time to reach the maximum peak during the flash flood which resulting in a short time for warning and evacuation.The produced distributions maps will be helpful for any local planning for the water management and disaster risk reduction by implementing the appropriate mitigation measures.The total available water during the flash floods events can be easily estimated at the target basin and consequently, appropriate mitigation and water management strategies could be put in action.

Figure 1 :
Figure 1: The study area selected for the comparison between GSMaP and rain gauge data.
with a best value of zero.Negative values indicate an overestimation by GSMaP, while positive values indicate an underestimation.(iii) The Nash-Sutcliffe efficiency (Eq.4) is a normalized statistic that determines the relative magnitude of the residual variance ("noise") compared to the measured data variance ("information")(Nash and Sutcliffe, 1970) .NSE indicates how well the Satellite data matches the gauge data and it ranges between negative infinity and unity; the latter indicating a perfect agreement.(iv) Coefficient of determination (R2) describe the degree of collinearity between satellite-Nat.Hazards Earth Syst.Sci.Discuss., doi:10.5194/nhess-2016-339,2016   Manuscript under review for journal Nat.Hazards Earth Syst.Sci.Published: 11 November 2016 c Author(s) 2016.CC-BY 3.0 License.

Figure 16 :
Figure 16: Karpuz River Basin and the digital elevation model at 1Km X 1Km spatial resolution matching the hydrologic model resolution.

Figure 17 :
Figure 17: (a) the number of floods in last decade and (b) flood-induced monetary loses in Turkey (Özgüler, 2003) .

Figure 19 :
Figure 19: Flow chart of data processing and calculations in Hydro-BEAM.
96)and relative variability of flows (1.08 and 0.98) with somewhat satisfactory performance for correlation coefficient (0.68 and 0.63).The KGE values were 0.66 and 0.62 for calibration and validation period, respectively.The simulated and observed hydrographs for the downstream point representing calibration and validation periods are shown in Fig.20.The hydrograph during the calibration period (Fig.20a) shows that the magnitudes of high and low flows are simulated well with the model however there is a minor shift in the timing of the events, as also indicated by the correlation coefficient values.During the validation period (Fig.20b), although the timing of the first high flow event matched well by the model, the magnitude of the Nat.Hazards Earth Syst.Sci.Discuss., doi:10.5194/nhess-2016-339,2016 Manuscript under review for journal Nat.Hazards Earth Syst.Sci.Published: 11 November 2016 c Author(s) 2016.CC-BY 3.0 License.

Figure 20 :
Figure 20: Discharge hydrograph of Calibration (a) and Validation (b) at Karpuz River.The square markers labeled (a, b, c, and d)

Table 1 . Selected rain gauges around the study area 10 Station Name Station ID Start Date End date Latitude Longitude Elevation
Global Satellite Remote Sensing Datasets of GSMaP were compared and validated with rain gauge-based precipitation data in the study area.Daily precipitation data from five meteorological stations located in the study area were obtained from the