Evaluation of remotely sensed precipitation product in a hydrological model of the Bednja watershed

In this paper, a statistical and spatial analysis of precipitation for the period 2000-2018 for the Bednja basin was performed, were the measured data from meteorological and/or rainfall stations of Croatian Meteorological and Hydrological Service (DHMZ) were compared with the data in form of remotely sensed precipitation product - CHIRPS (Climate Hazards Group InfraRed Precipitation with Station). The results of the analysis in the form of the annual sum, monthly distribution within the year and the spatial distribution and input data ratio over the basin show a good correlation between the measured and remotely sensed precipitation. In order to further evaluate the quality of the remotely sensed product, a SWAT hydrological runoff model was created.


Introduction
In recent times, we have increasingly been faced with extreme weather conditions, which may result in sharp rise in river water levels. Some of the recent events relating to Croatia, characterized by high water levels, are those experienced from May/June 2010 on the Danube, in November 2012 on the Drava and Mura, and in May 2014 on the Sava River. The theme of this paper is the Bednja River basin, which has already been the subject of research with the aim of making flood hazard maps, due to the presence of numerous local communities along its watercourse [1]. However, in this paper, the water wave volume was crucial for the conducted research, as an alternative approach in hydrological modelling due to the characteristics of smaller catchments, where water waves are often caused by precipitation over the entire catchment area [2,3]. Furthermore, due to negligible runoff coefficient variations, resulting in a constant ratio of effective to total precipitation [4], the water wave could in this case be a good indicator when comparing the precipitation over the basin. Several parameters are important for the formation of water waves, the most notable being: high-intensity precipitation, inflow water wave, spatial distribution of rain over the catchment area, type of land cover and soil type within the basin, relief, etc. The water wave is defined by several key characteristics: peak flow, water wave duration, water wave rise and fall time, and/or water wave volume [5,6]. In the analysis of precipitation over the basin or the development of hydrological models, the existence of a sufficient number of measuring points is very important, but their distribution within the basin is also highly significant [7]. However, the number of measuring stations is often insufficient in rural and hilly areas -partly due to sparse population density, and partly due to inaccessible terrain -which is a serious challenge for accurate hydrological modelling. As the Bednja River basin is not affected by the lack of measuring stations, it is quite suitable for assessing the possibility of applying the remotely sensed precipitation product (CHIRPS -Climate Hazards Group InfraRed Precipitation with Station) [8] in hydrological modelling using the SWAT (Soil and Water Assessment Tool) hydrological model [9]. The SWAT hydrological model has been widely applied for the purposes of hydrological modelling of basins [9][10][11], for measuring the impact of climate change on basins [12], and also for comparing different precipitation input data [13,14]. The aim of this paper is to compare precipitation traditionally measured at weather (rainfall) stations of DHMZ with the remotely sensed precipitation product -CHIRPS in the form of statistical analysis, time course, intra-year distribution and spatial distribution above the Bednja basin. The application of the SWAT hydrological model has not been widely investigated in Croatia, i.e. only one paper [15] deals with the SWAT model implementation in Croatia, while the application of remotely observed precipitation products in hydrological modelling is completely unexplored in Croatia. Therefore, the possibility of applying the SWAT hydrological model with the measured and remotely sensed precipitation input data was tested as part of this research on the example of the Bednja River basin.

Hydrogeological characteristics
The investigated area is situated in north-west of Croatia and includes the Bednja River basin (Figure 1), which is a part of the Drava river basin. In a broader context it is a part of the Danube basin. The main course of the Bednja River is about 103 km long [16]. The Bednja River runs in the west-east direction throughout its course, except for a very small deviation near Novi Marof where it first runs southward and then turns toward the north. Regarding the hydrogeological features, three basic parts of the basin can be distinguished [17,18]. The first part is the area of the underlying mountains (Ivanščica, Kalničko gorje and Ravna gora) which is built of carbonate rocks -mostly dolomite of secondary porosity dating back to the Mesozoic Era. The thickness of the carbonate aquifer is estimated at several hundred meters [17], and the permeability of these rocks is estimated at 3 to 25 % [19]. Other hydrogeological parameters have not so far been sufficiently studied, but it is known that folds and faults give basic structural features to this part of the catchment area. The second part is the area composed of sedimentary complexes of the Tertiary Period, which are the result of leaching and erosion of older rocks. These sediments are associated with hilly hills and smaller stream valleys on the slopes of the underlying mountains. The third part is the area composed of the Quaternary age deposits in the lowland area along the Bednja River. Today, this area is a fundamental feature of the Bednja River because gravel, sand and clays of Holocene age are deposited in the valley of the river and its tributaries, and they represent a significant alluvium of Varaždin County. This alluvium is primarily formed by geological weathering of older rocks of Mesozoic and Tertiary age and through accumulation of this material in river valleys. The base is dominated by gravel and sand, and the granules are predominantly made of quartz; there are also traces of eruptive rocks and fragments of limestone and dolomite [20]. From the hydrogeological aspect, alluvial sediments of the Bednja River form aquifers of intergranular porosity, of small thickness, and of laterally and vertically heterogeneous composition. Since fine-grained particles dominate, the permeability of these deposits is low, and groundwater recharging occurs primarily by infiltration of precipitation and seepage from watercourses [21]. The general flow of groundwater within the Bednja River basin follows the direction of the Bednja River flow and the groundwater flow is on a regional scale controlled by recent structural features, i.e., by faults and geological structures [22].

Hydrological characteristics
The Bednja River basin is defined by a topographic catchment. The catchment area of 596 km 2 and the mean river basin width of about 5.8 km were defined in previous investigations [17]. The catchment area is dominantly elongated and spreads from Trakošćan in the west to Ludbreg in the east. It is a predominantly lowland area bordered on the west by Maceljsko gorje mountains, on the south by Ivanščica mountain and Kalničko gorje mountains, and on the north by smaller hills that represent a water divide towards the Plitvice watercourse. To the east is the confluence with the Drava River. The very beginning of the Bednja River course is not precisely defined because there are various literary references -a stream that emerges at Brezova gora, a source in Ravna gora near the municipality of Bednjica, a source in Kamena gorica, an outflow of water from Lake Trakošćan. The confluence of the Bednja and the Drava River is located near the settlement of Mali Bukovec. In its longitudinal profile, the Bednja River is predominantly a lowland river with a very small longitudinal slope. However, the associated basin is almost exclusively hilly [16]. Along the course of the Bednja, one can distinguish the upper course (to the mouth of the Željeznica), the middle course (to the mouth of the Velinečki creek), and the lower course to the mouth of the Drava River. The upper course is fan-shaped, which together with the existing rainfall intensity results in a large and rapid concentration of inflow that regularly forms flood waves [1]. In the entire upper and middle course, the Bednja riverbed is formed of clay material with added mixtures of sand and mud. The material of the riverbed is highly susceptible to erosion, and so the faster flow was solved by widening the profile, as well as by cutting meanders. In its lower course the Bednja runs through a distinctly flat area, mostly composed of arable land and meadows. Erosion processes are present throughout the basin. There are five active hydrological stations in the Bednja River basin (Table 1) According to hydrological research conducted so far [1], the Bednja is a very torrential watercourse, and the flow directly depends on the volume of precipitation, but also on the melting of snow in the area around the source of the Bednja. The highest runoff occurs in spring (March and April) when the snow melts and spring storms are at their strongest [23].
Storm showers (usually occurring in August and September) are also common, causing sudden extreme increase in flow, as well as prolonged and heavy rains in September and October. By analysing the flow at hydrological stations along the Bednja River in the 2000 -2018 period (Table 1) The conducted analysis reveals that the number of extreme events in autumn has been increasing in recent years, which partly shows a somewhat greater impact of rain, and somewhat lower impact of snow on the runoff of the Bednja River.

Weather data and characteristics
Weather data for the period from 2000 to 2018 were obtained by the Croatian Meteorological and Hydrological Service [24]

SWAT model and input data
SWAT (Soil & Water Assessment Tool) [9,25] is a physics-based hydrological model developed for the purpose of evaluating changes in hydrological processes, river erosion and water quality in large basins [10,26]. Model simulation is possible in annual, monthly, daily and hourly time steps. Structurally, the model divides the basin into sub-basins and further into hydrological units -HRU (Hydrologic Response Units), which are characterized by unique properties related to land use and cover (Land Use / Land Cover), soil type (Soil), and topographic characteristics -Slope. The hydrological cycle within each HRU is determined by water balance (Eq. 1), which is controlled by weather parameters, precipitation, maximum and minimum daily air temperature, relative humidity, wind speed, and solar radiation. (1) where: SW t -the volume of water in the unit volume of soil at time t, SW -the volume of water in soil at time t-1 t -time R -the precipitation Q -the runoff ET -the evapotranspiration P -the percolation into deep aquifer, QR -the seepage from shallow aquifer to the watercourse. All parameters are in mm [9].
While creating the model, the SWAT2012 version was used together with the ArcSWAT 2012.10.23 plugin for ArcMap, which enabled creation of a SWAT model of the basin in a graphical interface. The list of input data and sources is given in Table 3.
A digital elevation model with a spatial resolution of 25 x 25 m was used as an input parameter for determining the direction of flow and division into sub-basins. Land cover / Land use data (Corine Land Cover) and the data on soil type according to FAO85 classification were used to determine the HRUs. In addition to the land cover and soil type data, SWAT uses a digital elevation model for the classification of terrain slope when dividing into HRUs. Given the predominantly lowland character of the basin and a moderate difference between the highest and lowest points in the basin (approx. 900 m), the slopes of the terrain were divided into three classes. Thus, the Bednja river basin was divided into 12 sub-basins and 148 HRUs.

Filling missing time-series data
Since SWAT simulations require continuous input of data, it was necessary to find a sufficiently long period without significant gaps in the input of data and, when necessary, to fill in the missing data, either by applying appropriate interpolation procedures or by selecting other available data.  Ludbreg, it was established that no other station can be used as its replacement. However, an additional analysis showed a very good correlation with the nearby weather station Varaždin (correlation coefficient 0.95) and so the data from that station were used to supplement and interpolate the data missing at the Ludbreg station ( Figure 2). Solar radiation data (W/m 2 ) were not available (insolation data were obtained, which is not enough to produce a SWAT model), and so the option of simulating data according to the database was used. Namely, the entry of the required climate data into the SWAT model is possible as a direct time series or as a simulation. SWAT simulations of climate data are performed based on databases and, in this case, the Climate Forecast System Reanalysis [29] (hereinafter CFSR) database was used, as it had been successfully used in several studies [14,[30][31][32] as the only input parameter in the SWAT model.

SWAT-CUP calibration tool
The SWAT-CUP (SWAT -Calibration Uncertainty Programs) [33] is a program designed to calibrate and assess sensitivity parameters (Sensitivity Analysis) of the SWAT model. The SUFI-2 (Sequential Uncertainty Fitting 2) method was used in this paper. The SUFI-2 offers up to 10 different statistical indicators for the calibration of quality. in this paper, the coincidence of modelled and observed flows was used, and the coefficient of determination R 2 (Eq. 2) and coefficient NSE (Eq. 3) (Nash-Sutcliffe efficiency) were selected as the main indicators [34]. R 2 [35] significantly depends on peak flows, and so it was chosen as a significant indicator of the quality of torrent models such as that involving the Bednja.
where: Q m, i -the measured runoff over time i -the mean observed runoff Q s, i -the modelled runoff over time i -the mean modelled runoff. NSE (Eq. 3) determines the importance of residual variance between the modelled and observed runoff [34,35].
where: Q m, i = Q m -the observed runoff over time i Q s -the modelled runoff over time i m Q -the mean observed runoff.
All R 2 [35] and NSE [36] values above 0.5 are considered satisfactory, while NSE values between 0.75 point to an extremely reliable model [37]. The p-factor and the r-factorparameters directly related to SWAT-CUP -were also evaluated. The P-factor is a measure of the amount of data within the 95PPU band (simulated flow values between 2.5 and 97.5 % confidence), while the r-factor indicates the width of the 95PPU band [10]. According to [10], p-factor values between 0.7 and 0.75 are considered adequate in runoff modelling [10] while, according to the same author [38,39], r-factor values of around 1.5 are considered desirable.

Statistical precipitation analysis
Mean values of remotely observed precipitation for the Bednja basin were compared (arithmetic mean of 75 raster points within and along the basin) with the precipitation measured at GMP Varaždin, and at DHMZ climatological and precipitation gauge stations (arithmetic mean of 14 stations (Figure 1) for the period 2000 -2018) in order to assess the possibility of applying the remotely sensed precipitation (CHIRPS product). Statistical daily precipitation values (Table 2), mean values of measured precipitation at the rain gauge stations (METEO_AVG), and mean values of CHIRPS remotely observed precipitation (CHIRPS_AVG) were also determined.      (Table 4). Furthermore, the measured precipitation has extremes (minimum and maximum values) that are more pronounced, while the quartiles are very similar for both modes of observation. Graphical representation of the total annual precipitation (Figure 3) additionally confirms the good match between the mean remotely observed precipitation (CHIRPS product) and the arithmetic mean of precipitation observed at weather stations and/or precipitation stations. It should be noted that there is a certain deviation of remotely sensed precipitation in rainier years (2010,2013,2014). The median amount of mean monthly precipitation at the measuring stations is almost identical to the median mean monthly precipitation of the remotely observed product (Figure 4), while the arithmetic mean is slightly lower for the remotely sensed precipitation compared to the observed precipitation. The analysis of monthly precipitation values throughout the year (Figure 4) shows that the hydrological year in the observed basin is characterized by two periods, i.e. the wet period (May -October) and the dry period (November -April). The wet part of the year is characterized by the months with total monthly precipitation in excess of 70 mm, locally reaching as many as 120 mm in September at some stations, while average monthly values during the dry part of the year range between 60-65 mm per month. The month with the highest precipitation is September (P_m_September = 113 mm), while the lowest values are recorded in January (P_m_January = 51 mm) and March (P_m_March = 58 mm).

Spatial precipitation analysis
A spatial analysis was conducted in order to show more clearly the variation of precipitation values over the basin. The spatial distribution analysis consists of the analysis of the mean daily precipitation in the observed period, the analysis of the mean daily precipitation during the dry season, and the analysis of the mean daily precipitation during the wet part of the year ( Figure  5) for each of the observation/measurement methods, which makes up a total of six components. For this purpose, the average daily precipitation in the observed period was determined on an annual basis, during the dry part of the year, and during the wet part of the year. The results were calculated using the Python 3 script, while the visualization was conducted in ArcGIS using spatial interpolation and the kriging method [40]. Visual interpretation shows a similar spatial distribution of precipitation values above the basin, with a slightly larger range of minimum and maximum mean daily precipitation values for the remotely sensed product (1.74 mm -3.48 mm) compared to climatological measurements (1.88 mm -3.42 mm). Furthermore, the mean precipitation on an annual basis ( Figure 5 -top) and during the dry part of the year (Figure 5 middle) generated with the remotely observed product resulted in a much more uniform spatial distribution of precipitation measured above the basin, which may be due to more points in the grid during interpolation. In meteorological measurements, spatial variations are somewhat more pronounced along the western hilly parts and along the eastern lowland part around Ludbreg. In the wet part of the year, the situation changes, i.e., the spatial distribution of the remotely sensed product is very similar to that of meteorological measurements -with relatively greater precipitation along the western hilly part of the basin, while downstream precipitation decreases linearly with both types of input data. In the final step, the ratio of the remotely sensed product (CHIRPS) to precipitation measurements was analysed, so that the raster ratio of remotely sensed product CHIRPS and precipitation measured at DHMZ stations was calculated. The coincidence of the remotely sensed product with the measured precipitation is better when the ratio is closer to 1. A ratio greater than 1 points to the overestimation of the remotely sensed product in relation to weather measurements, while a ratio smaller than 1 indicates areas where the remotely observed product shows lower precipitation values compared to weather measurements. With average daily precipitation at the annual level and when observing the dry period, a very good overlap can be observed in the central part of the basin, around the area of Novi Marof, while along the marginal parts, the remotely sensed product shows lower precipitation values. When observing the wet part of the year, the situation is reversed, and so in the central part the overestimation of (blue) remotely sensed product (CHIRPS) is visible, while in peripheral parts of the basin the ratio is around one, which indicates a very good overlap between remote observations and weather measurements. (Figure 6).

Hydrological runoff model results and discussion
The possibility of applying the remotely sensed precipitation (CHIRPS product) in the observation of the water wave volume was tested using the SWAT hydrological model. As very often the model does not give a satisfactory result in its initial form, it is important to calibrate the model [10,41,42], and the result is then validated on an independent time period that is not covered by previous calibration. Before performing such calibration/validation, it is important to select appropriate time periods. Therefore, the hydrograph of the most downstream hydrological station in the hydrological model (Ludbreg) was analysed. For the high quality and reliability of the calibration/validation process, it is very important to select periods that contain a similar number of wet and dry years, i.e., they should be statistically as similar as possible [5,20] (and personal communication with Dr. K. Abbaspour). Wet years are characterized by hydrographs with the pronounced peak flow and a pronounced duration (higher volume of the water wave), while dry years are characterized by hydrographs with a less pronounced peak flow and, in general, with a very low base flow [2].   Evaluation of remotely sensed precipitation product in a hydrological model of the Bednja watershed To ensure optimum results, model calibration was performed at three points (measuring stations): Lepoglava, Željeznica, and Ludbreg. The reason for this is to consider the lower flows in the upstream part of the basin (Lepoglava station), the boundary between the upper and middle part of the river (Željeznica), and the last measuring station on the Bednja River (Ludbreg), which is also the outlet of the hydrological model. The process of calibrating the model with SWAT-CUP begins with Sensitivity Analysis [44]. In this step, the model is simulated with several parameters that include the processes observed in the basin, considering the physical characteristics of the basin, and especially the ranges of selected parameters to find the most sensitive ones (Table 5). This information is then used to calibrate the model [10,38,44,45]. Calibration of SUFI-2 parameters by the SWAT-CUP software algorithm can be performed in one of three ways: relative change of r__, absolute change of a__, and replacement of parameters with a new value of v__ [33,45].  (Table 6) and the parameters of the simulation that gives the smallest statistical error are selected as the "best simulation" (red line - Figure 7), while the whole range of possible runoff values depends on the parameters shown by the area 95PPU (green - Figure 7.). After analysis of statistical calibration and validation parameters presented in Table 6. and according to the proposed values of the coefficients R 2 , NSE, p-factor and r-factor for SWAT hydrological model [10,[35][36][37][38][39], it can be concluded that the model provides a reliable result. This is supported by the fact that the NSE coefficient is greater than 0.5 [35].
The assumption is also that extending of the calibration period would average the results,  Karlo Leskovar, Petar Mrakužić, Hrvoje Meaški but this is beyond the scope of this paper. It is interesting to note that in the calibration process, the model with weatherly observed precipitation gives a slightly worse result, while during the validation period the situation is reversed.
Compared with similar studies in other areas [13], the obtained results show that the remotely sensed product CHIRPS could be a valid input to the SWAT hydrological model for small basins in Croatia, in the cases when a low density of precipitation gauging stations, or interruptions in time series, are present. However, this approach to modelling has certain limitations, which are noticeable in the visual interpretation of the runoff hydrograph of the SWAT hydrological model, as shown in Figure 7. It can be seen that the meteorologically measured precipitation as an input parameter of the SWAT hydrological model better describes peak flows in wet years (i.e., 2013 or January 2015), where CHIRPS, due to an error in the observed precipitation volume (Figure 3.), does not achieve the actual peak flows recorded in that month. If we look at the years in which the CHIRPS product nicely overlaps with volumes of meteorologically measured precipitation, the error recorded in the model is also lower. The situation changes in dry years (periods), when the CHIRPS product, due to a larger number of observation points (raster), describes much better the minimum amounts of precipitation, which is reflected in a better representation of the base flow. This feature is also an advantage of remotely observed products in areas with extremely low density of weather stations. Due to that, they can give much more reliable results than the meteorologically measured precipitation, because they can transfer the spatial distribution of the precipitation over the basin in a more realistic manner.
The obtained data show that the CHIRPS product describes well the time of rising and falling of the water wave in both wet and dry years, which is evident from the comparison of the mean monthly precipitation above the basin (Figure 4 -graphs 15 and 16).

Conclusion
An alternative approach to hydrological modelling is presented in this paper, for cases in which the number of precipitation stations in the catchment area or in the immediate vicinity is insufficient, or when the existing time series have interruptions in measurements. In such complicated situations, a remotely sensed product with global coverage, i.e. CHIRPS, could be a valid replacement for meteorologically measured precipitation. It should be noted that the reliability of such models is slightly lower compared to the model in which precipitation observed at weather stations and/or precipitation stations is used as an input parameter, bearing in mind that the Bednja basin is very well covered by the DHMZ gauge network. The question of the situation in the basins where the density of precipitation gauge stations is lower remains open, which is a topic for further research. According to the obtained results, it can be assumed that the accuracy could be further improved by calibrating the original model over a longer time period of precipitation observed at weather stations and/or precipitation stations.