Analysis of a Distributed Hydrological Model Response by Using Spatial Input Data with Variable Resolution F. Soria1, S. Kazama2, M.

Sawamoto1 Graduate School of Engineering, Tohoku University, Japan; 2 Graduate School of Environmental Sciences, Tohoku University, Japan soria@kaigan.civil.tohoku.ac.jp
1

1. Introduction The Stanford Watershed Model was in the 1960s one of the first programs developed to predict streamflow through simplifications and idealizations towards understanding the components of the hydrological cycle, the physical environment, and all interactive and complex interactions that link them together. From the engineering point of view, the interest may be focused in the temporal variation and spatial distribution of water quantity, water quality and sediment load, because of its importance in water management policies development and currently increasing concerns about climate and land use changes and its effects in water resources availability. Distributed hydrological models can be considered as a useful tool to analyze and quantify catchments response because its capacity to account spatial variability of watershed processes, inputs, boundary conditions and physical characteristics, only if spatial input data demand is adequately fulfilled. Such issues become difficult to accomplish in poorly and ungauged basins specially in developing countries, where measurement networks density have decreased in the last decades mainly because the high cost of maintenance and monitoring, and where geographical, geological and physical information is also scarce to accomplish certain modeling resolution requirements. In such frame, some initiatives such as the PUB (prediction of ungauged basins) of the International Association of Hydrological Sciences have been considering situations to better represent the catchment response in poor data situations. Heterogeneity of the land surface condition, soil structure, vegetation, land use, etc, and the spacetime variability of climatic inputs occurring over a wide range of space and time scales, and which never seems to disappear at whatever scale we observe (Sivapalan, 2003) could be seen as the main problem for hydrologic predictions. Its representation in distributed models is difficult to overcome and the best suited modeling scale should theoretically follow the definition of Singh (1995) about the term scale itself as the size of a unit or a subwatershed within which the hydrological response can be treated as homogeneous . Many papers analyze the effects of such heterogeneity by introducing variation in the scale of input data. Kalin et al. (2003) uses the Kinematic Runoff and Erosion Model KINEROS (Woolhiser et al., 1990) to show that as the geomorphologic resolution decreases the peak runoff also decrease for all simulated runoff hydrographs; moreover, the author shows that time to peak does not change appreciably with increasing resolution, except for the pure overland flow case where time to peak is larger than expected.. Since it is clear that heterogeneity becomes more evident as smaller the scale considered, resolution issues need more attention for smaller basins (Shresta et al., 2006). Grid size effects can also be evaluated on catchment parameters. Yin and Whang (1999) estimate some drainage basin parameters for 20 basins for the US Geological Survey (USGS) 1:24 000 (30m) Digital Elevation Models DEMs and the USGS 1:250 000 (approximately 92m) DEMs. They found

that for the coarser resolution considered parameters as the mean maximum basin elevation, the total stream length and main stream length, the drainage density and in general slope parameters were underestimated; the opposite was observed for the mean minimum basin elevation and the relief ratio. Jetten et al. (2003) analyze the effect of 5m, 10m, 20m and 50m grid size in simulated discharge by the application of the Limburg soil erosion model LISEM. There it is shown that increasing grid size has indirect relationship to both peak discharge and runoff volume which decrease. The author points to the change in slope as the main responsible for lower water depth and velocity values obtained. LISEM model is also studied by Hessel (2005), who analyzes modeling results for grid cells ranging from 5 to 100 for a single time step, and for time steps ranging from 2 to 120 s with single grid size. Simulated discharge results show similar behavior in peak discharge as cited from previous researchers; besides, the author establishes that the numerical errors effect in the kinematic wave solution might be an additional cause for flood dispersion and infiltration time increasing. Predictions of hydrological response as river discharge are linked to issues cited above. Seems clear that higher data resolution would represent better watersheds heterogeneity than the opposite one, although is also clear that such high-quality data may not be always available. This paper aims to evaluate the variations that could be observed in runoff simulation as a way to contribute in the analysis of modeling uncertainty in various scenarios. A distributed hydrological model is applied for different spatial input data resolutions. Then, observed and simulated hydrographs are analyzed through visual and statistical comparison. No considerations have been focused on analyzing climatic input data resolution or time step calculation interval, which are maintained invariable. Final conclusions and recommendations are based on such frame. 2. Model description The model is developed under the structure suggested by Tsuchida et.al (2003). It aims to describe spatial variability of hydrological processes in three reservoirs, as a way to estimate river discharge in a watershed. Snowmelt contribution to total discharge is calculated using the degree day method, where the relation snowmelt - temperature is
M = M f (Ti Tb )

(1)

where M is the depth of meltwater in a period of time, Ti and Tb are the index and base air temperature respectively, and Mf is the melt factor. Commonly used values for Ti and Tb are the mean daily temperature and 0 C respectively. Mf is determined empirically by Totsuka (2003); such values are taken in this study. To describe overland flow routing of runoff is used the classical kinematic wave equation. It is derived from the continuity equation (equation 2), considering lateral inflow or outflow per lineal distance along the watercourse q as precipitation and snowmelt, besides water loss as evaporation. Then, the equation can be expressed in terms of finite-differences as follows

A Q + q = 0 t x Bh Q + = (r S m E ) B t x

(2) (3)

where A is the flow cross-sectional area with constant surface width of flow B, depth of flow h, time t, discharge Q, rainfall rate r, snowmelt rate Sm, evapotranspiration rate E, and x as the distance along the longitudinal axis of the watercourse. Infiltration estimation simply considers constant parameters estimated assuming a homogeneous soil structure for the entire basin. Infiltration capacity varies with a soil saturation rate controlled with the number of days with rain. The volume of storage S at any time in each cell of the watershed is related to a representative discharge q in the element at the same time as in equation 1. Moreover, Yoshikawa (1966) suggests empirical values of 40.3 for k and 0.5 for m. Then, storage depth is estimated with a storage function method (Yoshikawa, 1966): S = kq m (4) where k is a dimensional parameter and m a dimensionless.

dS o = r q dt
where So is the storage depth, and q the discharge rate.

(5)

Dynamic wave model is used for channel flow routing. The method employs St. Venant equations describing the laws of mass conservation (equation 2) and momentum conservation (equation 6). Later, the numerical solutions for both are written in equations 7 and 8.
2 1 A 1 v 2 H n v v + + + 4 =0 g t 2 g x x h 3
j +1

(6)
t vi 1 B q

hi = j 1 hi

t x

j 1

hi + j 1 hi +2 2
j

vi +1

j 1

hi 2 + j 1 hi 2
j

(7)

j +1

vi =

j 1

vi (1

gn 2 j vi t 2h
4/3

t t 2 2 ( j 1 vi +2 j 1 vi 2 ) g ( j H i +1 j H i ) 4x x gn 2 j vi t 1+ 2h 4 / 3

(8) where g is acceleration due to gravity, H water level, v flow velocity, and n Mannings friction factor. Besides, in the x-t plane, i refers to x spatial axis and j refers to time axis.

3. Studied watershed and input data The study is developed in the Natori River basin, center of Miyagi Prefecture, northeastern Japan (Fig. 1). The contributing area is 1015km 2, somehow bigger than the 970km2 showed in other studies since it is considered a major proportion of plain land. The main stream length is 55km (Natori River). Two large-scale dams (Kamafusa Dam and Ookura Dam) are located upstream Goishi River and Ookura River respectively. The basin is approximately 72% a mountainous region and 26% plain-land, with elevations that vary from 0MSL to 1470MSL. The average river slope is 40% in the mountainous region. Meteorological data is available hourly for three stations: Kawasaki, Nikkawa and Sendai, from the Automated Meteorological Data Acquisition System (AMeDAS) database. Inverse distance method is used to estimate spatial distribution of precipitation and temperature. Discharge records in daily frequency are available at the outlet of Kamafusa Dam, Ookura Dam and Yokata gauging station. Input grid data consists of elevation, river network, flow direction for north-south and west-east directions, flow accumulation, land use and mean monthly evapotranspiration raster maps. Land use is characterized by the Geographical Survey Institute and considers mainly urban areas, forests and crop fields. Evapotranspiration grid map is calculated indirectly from the Normalized Difference Vegetation Index NDVI, estimated from the National Oceanic and Atmospheric Administration-Advanced Very High Resolution Radiometer NOAA-AVHRR sensor (Nourbaeva et al, 2003). The 250m grid size input data is used to derive the correspondent 500m and 1000m inputs. Elevation data, land use and evapotranspiration data are derived from the correspondent 250m grid by bilinear interpolation, and then manually corrected for non-sense values. Flow direction grid is derived from the correspondent corrected elevation grid, and determined by the D-8 method finding the steepest direction of each cell and its eight neighbors; later, flow direction grid is manually corrected to force and allow the main elements of the original network to appear in coarser resolution maps. Flow accumulation is derived from corrected flow direction grids. Figure 2, 3 and 4 shows variations in grid size for the elevation and river network maps. The model is applied for 250m, 500m, and 1000m grid resolution, with simulations carried out for the period July 1999 to June 2000, only for simplicity. Time step calculation is fixed to 10s, and time step of meteorological input data is fixed hourly in the modeling process. Calibration is developed in the original 250m grid size based on comparisons with observed discharge at Yokata station. Calibrated parameters of infiltration, groundwater storage function, surface flow and snow melt are maintained invariable for the simulations of 500m and 1000m grid resolution. Initial conditions for mean water surface velocity, surface and groundwater water depth are estimated from results of an initial testing simulation. The results from a second simulation are the ones considered for analysis.

15 km

Figure (Shiraiwa et al., 2005)

1. Theof Natori River basin 1.Figure Location Natori River

basin

Figure 2. Natori River basin system. 250 m grid size

Figure 3. Elevation grid for a) 500m resolution, b)1000m resolution

4. Results and discussion The distributed model is calibrated using the Yokata gauging station discharge records for the period July 1999 to June 2000, assuming (later verified) that the best model performance would be for the 250m grid size resolution and time step calculation of 10s. The calibration considers additional runoff contribution from discharge of Kamafusa Dam and Okura Dam to allow the model describe an extraordinary event observed in August 1999 (Fig. 7). Other periods do not consider such in the simulation. Simulations for the 500m and 1000m grid size were developed with calibrated values of the 250m grid size simulation. Model performance examination is made for the period July 1999 to November 1999, to simplify the analysis avoiding at this time the evaluation for winter season and the correspondent snow melt contribution in spring.

Table 1. Catchment and river network

general characteristics for different resolutions

a) Hirose River
350 300 Altitude [m sl] 250 200 150 100 50 0 0 5 10

b) Natori River

0 5 15 20 25 30 Horizontal longitude [km ]

10 35

15 40

20 45

250 m grid size 500 m grid size 1000 m grid size

Catchment characteristics Contributing Elevation[msl] Grid Slope[%] Area [km2] Max. Mean Max. Mean 1014.88 1470 327 79 13 984.25 1412 322 47 8 943.00 1386 324 27 5

River network characteristics Stream Elevation [msl] Slope [%] length [km] Total 82 79 70 Max. Mean 250 41 278 47 300 55 Max. Mean 4.0 0.3 8.6 0.6 8.2 1.0

Altitude [m sl]

400 300 200 100 0 0

250m grid 500m grid 1000m grid 250m grid 500m grid Hirose river 1000m grid

Natori river

10

20

30

40

Horizontal longitude [km ]

Figure 5. Longitudinal section profile in a) Hirose River and b) Natori River for different resolutions

Analysis of general characteristics of the catchment and main streams might be useful to deepen the analysis of variations in discharge. The method used to create low grid resolution inputs has direct effect in some variables as catchment area, average and maximum elevation, and corresponding average and maximum slope values (table 1). Expected effect in runoff response might be linked to a reduction of total discharge as a consequence of catchment area reduction, besides a decreasing effect in runoff velocity and for instance in peak times arrival. The method used for deriving coarser grid maps has important effect when delineating the stream network (Fig. 2, 3, 4), which length generally decreases. In the 1000m grid size elevation and flow direction maps main rivers are more difficult to

identify automatically and requires manual corrections which are applied only to main network vicinity to avoid further deterioration of geographical features, which cause effects shown in Figure 5. Elevation and slope values have a tendency to increase, contrary to the effect observed by Jetten (2003) which supports his final conclusions. Errors related to re-location of required river network may have influenced in the appreciations of elevation and slope values. As stated later, slope reduction influence is not as expected, perhaps because the grid size considered for the analysis is too big, or because an underestimation in other streams of the network (e.g. Fig. 5, right).
700

Observed discharge Qi Sim ulated discharge Pi, 250 m grid Sim ulated discharge Pi, 500 m grid Sim ulated discharge Pi, 1000 m grid

600

Discharge Q(m / s)

500

400

300

200

100

1- Sep

1- Aug

Aug-99

Sep-99

Oct-99

1- Nov

1- Oct

Nov-99

Figure 6. Simulated discharge at Yokata gauging station. Period: August 1999 to November 99 Table 2. Monthly average values of total discharge at the outlet of Yokata gauging station for different resolutions
Jul 99-Nov 99 Total Observed 5952.3 Discharge [m3/s], 250 m 5353.9 Discharge [m3/s], 500 m 5478.3 Discharge [m3/s], 1000 m2303.7 Mean 38.9 35.0 35.8 15.1 Jul Total 1648.6 1246.5 1175.3 446.9 Mean 53.2 40.2 37.9 14.4 Aug Total 1842.9 1363.2 1772.6 1704.6 Mean 59.4 44.0 57.2 27.9 Sep Total 1347.7 1412.6 1291.8 314.9 Mean 44.9 47.1 43.1 10.5 Oct Total Mean 686.8 668.8 588.5 123.3 22.2 21.6 19.0 4.0 Nov Total 426.3 662.8 650.1 151.7 Mean 14.2 22.1 21.7 5.1

Simulation results show a decreasing tendency of discharge estimation, except for August 1999 where an extraordinary event occurred (Figure 2). The influence of topographical features in model results can be supported when comparing total and mean monthly values. Hessel (2005) works in a small catchment (3-5km 2) with higher resolutions than those used in the current study, where effects cited are observed; besides he establishes that the infiltration would also be beneficiated since more time is available for the process. Although such effect is not studied in detail here, infiltration contribution to baseflow does not seem to be important in the 1000m grid simulation where is observed an important tendency to underestimate total discharge. Visual comparison between observed and simulated hydrographs can be made (Fig. 6) together with statistical indices to further analyze time series behavior (Table 3). The Pearson Moment Correlation coefficient PM, the Nash-Sutcliffe Coefficient NS, the Index of Agreement IA and the Root Mean Square error

RMSE are now considered. Values of PM, NS and IA indices close to one and/or RMSE close to zero would mean time series coincidence.

PM =

[(O
[

(O
i

Om )( Pi Pm )
0.5 i

Om ) 2

] [( P P ) ]
m

2 0.5

(1)

IA = 1 -

Pi Om + Oi Om ]2

(O

Pi ) 2

(3)

(O P ) (O O ) (O P ) RMSE =
2

NS = 1 -

(2)
2

(4)

where Oi is the observed value, Pi the simulated value, and O m and Pm are the average values.
Table 3. Statistical indices comparison of discharge simulated at Yokata gauging station for different resolutions
Jul 99 PM NS IA 250m grid 0.675 0.334 0.669 500m grid 0.335 0.005 0.500 1000m grid 0.427 -0.47 0.484 Aug 99 PM NS IA 0.967 0.548 0.752 0.477 0.136 0.653 0.483 0.206 0.617 Sep 99 PM NS IA 0.936 0.848 0.951 0.549 0.268 0.723 0.572 -0.39 0.481 Oct 99 PM NS IA 0.948 0.783 0.915 0.261 -0.03 0.458 0.278 -0.19 0.314 Nov 99 PM NS IA 0.780 -3.69 0.574 0.642 -4.87 0.527 0.915 -4.05 0.462

General performance of simulation for coarser grids tends to be towards underestimating discharge. An exception is found when simulating the event of August 1999, where the response to direct contribution from dams is more notorious in larger grids for the peak observed. Despite such anomalies, the simulation for the calibrated grid size has an overall description much closer than the simulation results for coarser grids, as should have been expected.
Table 4. Maximum discharge and time to peak at Yokata gauging station for different resolutions
Jul Max. Time to discharge peak [day] [m 3 /s] Observed 183.8 12 250 m grid 110.0 12 500 m grid 106.9 13 1000 m grid 41.7 no peak Aug Max. Time to discharge peak [day] [m 3 /s] 672.2 14 224.4 14 519.5 15 388.4 15 Sep Max. Time to discharge peak [day] [m 3 /s] 178.1 15 176.3 15 162.9 16 33.8 16 Oct Max. Time to discharge peak [day] [m 3 /s] 205.9 28 130.7 28 120.2 29 26.2 29 Nov Max. Time to discharge peak [day] [m 3 /s] 28.1 no peak 35.8 no peak 39.0 no peak 16.8 no peak

Peak discharge and time to peak follow the same direction than the ones analyzed before. The coarsest grid simulation tends to have difficulties when describing discharge peaks, with an exception in the case of August 1999. Both 500m and 1000m grid size simulations present delay in peak arrival. That could be explained remembering what seen before, related mainly to the general tendency to underestimate discharge perhaps because flow spread caused by the increasing size of individual channels of every grid cell as grid size is increased. 5. Conclusions Results show that the effect of grid size in distributed hydrological modeling has the expected influence in model performance. Main characteristics observed lead to establish that as the gird increases, discharge values decrease, and peak flow estimations are every time more difficult to overcome. The affirmation of Shresta (2006) respect to the need of more attention when working in smaller basins seems to be evident.

Jetten (2003) establishes additional factors to explain discharge decreasing tendency as grid size resolution grows based on known practical facts referring to the limitations of the kinematic wave theory and the kind of numerical method applied for its solution (Maidment, 1993). The author cites the effect that would have the kinematic wave theory and the numerical method used for its solution in flow routing description. Although no further analysis is shown at the time, his affirmation should lead to further work to quantify the influence of kinematic wave numerical solution in the simulation of total discharge values. DEMs preparation should get special attention to overcome issues related to the inaccurate description of catchments characteristics, as shown in the results of the current where despite the overestimation of elevation and slopes values in the river network discharge values remained below observed ones. Prediction in ungauged basins may find in this point one of the main sources that leads to uncertainty in hydrological response predictions, for which further work to overcome it is needed. Meteorological input resolution is an additional issue where the minimum requirements and the influence of the variations in the time step should also be analyzed to contribute in the construction of the solution of this interesting field. References
Chaplot V. Impact of DEM mesh size and soil map scale on SWAT runoff, sediment, and NO3N loads predictions, Journal of Hydrology 312 (2005), pp. 207222. Hessel R. Effects of grid cell size and time step length on simulation results of the Limburg soil erosion model (LISEM), Hydrol. Process. 19 (2005), pp. 30373049. Jetten V., Govers G., Hessel1 R. Erosion models: quality of spatial predictions, Hydrol. Process. 17 (2003), pp. 887-900. Kalin, L., Govindarajua R.S. and Hantush M.M. Effect of geomorphologic resolution on modeling of runoff hydrograph and sedimentograph over small watersheds, Journal of Hydrology 276 (2003), pp. 89111 Maidment D. Handbook of Hydrology. Mc Graw Hill, USA (1993), pp.10.9 10.11, 7.25. Nourbaeva G., Kazama S., Sawamoto M. Assessment of Daily Evapotranspiration Using Remote Sensing Data, Environmental Informatics Archives, 1 (2003), pp. 421-427. Shiraiwa J, Kazama S, Sawamoto M. Analysis of river temperature based on a hydrological model.XXXI IAHR Congress, Seoul, Korea (2005). Shresta , R., Tachikawa, Y.,Takara, K.. Input data resolution for distributed hydrological modeling. J. of Hydrology 319 (2006), 36-50. Singh, V.P.. Computer Models of Watershed Hydrology. Water Resources Publications, USA (1995), pp. 1-12. Sivapalan M. Prediction in ungauged basins: a grand challenge for theoretical hydrology. Hydrol. Process. 17 (2003), pp. 31633170. Totsuka T., et al. The Analysis of the Snow Water Equivalent Distribution Using Snow Model and Satellite Information of Snow in Tohoku District in Japan, Proceedings of 2003 Annual conference, JSHWR (2003), 230-231. Woolhiser, D.A., Smith, R.E., Goodrich, D.C. KINEROSa kinematic runoff and erosion model: documentation and user manual. US Department of Agriculture, Agricultural Research Service (1990), ARS-77, 130. Yin Z-Y, Wang X. A cross-scale comparison of drainage basin characteristics derived from digital elevation models. Earth Surface Processes and Landforms 24 (1999), pp. 557562. Yoshikawa H. River Engineering. Asakura shoten editions (Japanese version) (1966), pp. 45-48.