Effects of urbanization variables on model parameters for watershed divisions

Abstract

A semi-distributed model with a parallel connection was applied to examine the effects of urbanization variables. Data were obtained from watershed divisions that were characterized by various degrees of urbanization. The mean rainfall was calculated using the kriging method. The model inputs were obtained by subtracting mean rainfall from Ф-index values, based on the spatially uniform loss assumption. Regression analysis was applied to determine the relationship between the parameters of 64 calibrations and urbanization variables among the divisions. The results showed that overland parameters produced more consistent change in response to imperviousness than to population. Conversely, the channel parameter was unaffected by changes in urbanization. The verification results of 46 cases showed that power linkage was a potential option for linking division parameters with the corresponding imperviousness based on four evaluation criteria. The changes in imperviousness on overland parameters show the hydrological effects of division urbanizations.

Editor D. Koutsoyiannis; Associate editor T. Wagener

Citation Chen, R., Chuang, W.-N., and Cheng, S., 2014. Effects of urbanization variables on model parameters for watershed divisions. Hydrological Sciences Journal, 59 (6), 1167–1183. http://dx.doi.org/10.1080/02626667.2014.910305

Résumé

Un modèle semi-distribué à connexion parallèle a été utilisé pour examiner les effets des variables d’urbanisation. Les données ont été obtenues pour des parties de bassins versants caractérisées par différents degrés d’urbanisation. Les précipitations moyennes ont été calculées en utilisant le krigeage. Les entrées du modèle ont été obtenues en soustrayant les précipitations moyennes des valeurs d’indice Ф, en se basant sur l’hypothèse d’une perte spatialement uniforme. L’analyse par régression a été utilisée pour déterminer la relation entre les paramètres de 64 calages et les variables d’urbanisation des différentes parties considérées. Les résultats ont montré que les paramètres de surface ont manifesté des changements plus cohérents en réponse à l’imperméabilisation qu’à la population. En revanche, le paramètre du chenal n’a pas été affecté par des changements de l’urbanisation. Les résultats de validation sur 46 cas ont montré qu’une fonction puissance était une option possible pour relier les paramètres des parties considérées avec l’imperméabilisation correspondante sur la base de quatre critères d’évaluation. Les changements de l’imperméabilisation sur les paramètres de surface montrent les effets hydrologiques de l’urbanisation des parties considérées.

Key words:

• block kriging
• urbanization variables
• watershed divisions
• semi-distributed model
• runoff hydrograph
• watershed characteristics

Mots clefs:

• krigeage par bloc
• variables d’urbanisation
• parties de bassin versant
• modèle semi-distribué
• hydrogramme de débit
• caractéristiques des bassins versants

INTRODUCTION

Since people first inhabited the Earth, they have clustered in specific areas. Urban areas contained increasingly concentrated populations, and gradually developed into societies or cities. The migration of people results in urban development. The demand for human quality of life and a sophisticated environment resulted in the development of numerous impervious surfaces over the ground in populated areas. These impervious areas include schools, railroad lines, streets, roofs, parking lots, shopping malls, waterway, highways, and commercial and industrial buildings.

There are differences between the impervious areas found in rural spaces and those found in urban spaces. These varying degrees of urbanization may exert distinct influences on the upstream and downstream areas of a basin. The surface imperviousness caused by urbanization might alter the localization in a hydrological cycle. In the hydrological context, urbanization occurs when people increase the number of impervious surfaces in a geographical area. This increase primarily affects the infiltration mechanism of the Earth’s surface; that is, the infiltration of rainwater into the soil layer just below the surface is diminished. Simultaneously, rainwater loss and surface water levels are altered. Thus, the hydrographical characteristics of a developed watershed vary in conjunction with the degree of urbanization (Kliment and Matoušková 2009). Previous studies have shown that the affected factors include rainwater loss (Gremillion et al. 2000, Cheng et al. 2008b), watershed function (Krug 1996, Kang et al. 1998, Aronica and Cannarozzo 2000, Cheng and Wang 2002), surface runoff (Boyd et al. 1994, Junil et al. 1999, Rodriguez et al. 2003), runoff volume (Arnell 1982), peak discharge (Huang et al. 2008a, 2008b), time to peak (Huang et al. 2012) and baseflow (Simmons and Reynolds 1982). However, the hydrological effects of urbanization may be particularly severe in watershed divisions that are vulnerable to the breaking of links in the hydrological cycle. These variations in vulnerability are caused by non-uniform spatial changes in the migration of people, as well as varying degrees of impervious paving.

Hydrological modelling is not the only approach for analysing changes in urbanization (Olivera and DeFee 2007). Numerous studies have shown that hydrological routing remains a valuable tool for determining how changes in an outlet hydrograph affect a watershed system over time (O’Connell and Todini 1996, Melone et al. 1998, Bhadra et al. 2010). Hydrological models that include system conception (Dooge 1959) generally treat a watershed as an input–output transformation system. The system inputs and outputs are analogous to natural flows as rainfall and runoff. Mathematical formulations based on this approach have appeared in numerous models: for example, the Nash model (Nash 1957, Agirre et al. 2005), Clarke model (Clarke 1973, Ahmad et al. 2009), geomorphologic model (Jin 1992, Franchini and O’Connell 1996, Nourani et al. 2009), distributed models based on parallel-type reservoirs (Hsieh and Wang 1999), a model of three serial linear reservoirs (Cheng 2010a, 2010b, 2010c) and a linear cascade model of three serial reservoirs with one parallel reservoir (Yue and Hashino 2000, Li et al. 2012).

Researchers seeking to understand and develop solutions to complex urban problems frequently use a systems analysis approach towards runoff simulation. System analysis can facilitate urban water resource management by calculating volumes, flow rates and water quality (Liu et al. 2012), as well as urban storm runoff. The problems encountered in urban water systems, which are inherently distributed systems, should be analysed for both spatial and temporal variations. However, deriving a mathematical formulation for urban water systems distributed in both time and space is a complex task. Spatial variation can be removed by dividing an entire watershed system into several subsystems, which are considered as lumped system models. These models are subsequently linked to produce a complete system model (Agirre et al. 2005).

The goal of this study was to investigate how changes in two urbanization variables (impervious percentage and population density) affect the hydrograph parameters of four watershed divisions. These two variables were subsequently used as an index of watershed development. This study proposes a semi-distributed model that can be applied to assess a parallel-type cascade of linear reservoirs. The model is a conceptual unit hydrograph (UH)-based model, with storage constants K_o and K_c representing the overland and channel storages of watershed divisions, respectively. The mean values of rainfall on watershed divisions were calculated by using the block kriging technique. Rainfall losses were assumed to be spatially uniform, and were estimated based on the mean rainfall and direct runoff of an entire watershed by using the Φ-index method. The distributed effective rainfall was input to the model to complete the related calibrated model parameters with urbanization variables, and to verify the usability of these relationships derived from model parameters and urbanization variables. The results of this study show that the model parameters are affected more consistently by changes in imperviousness than by changes in population density. The overland storage parameter shows a marked change related to changes in imperviousness, whereas the channel storage parameter is independent of the urbanization process. Thus, the calibrated overland parameter values can be related to the imperviousness percentages, and the channel parameters are constants on watershed divisions. The relationships, which have the form of a power equation between model parameters and urbanization variables on watershed divisions, were confirmed by performing model verification. This confirmation shows that these linkages are appropriate to urbanization changes among the examined watershed divisions. Changes in hydrographic characteristics among these watershed divisions can be identified to represent various degrees of urbanization.

WATERSHED DESCRIPTION

Geographical features

One of the main tributaries of the Tamshui River, which is the third longest river in Taiwan, is the Kee-Lung River (Fig. 1(a)). The Wu-Tu watershed is located upstream of the Kee-Lung River, and was selected as the research site for this study. The Wu-Tu watershed covers nearly 204 km² and surrounds Taipei City in northern Taiwan (Fig. 1(b)). The mean annual precipitation and runoff of the entire Wu-Tu watershed are 2865 and 2177 mm, respectively. The watershed comprises a large pervious area (high mountains) and a smaller impervious area (watershed downstream), with the majority of the runoff flowing from the pervious area. The rugged topography of the watershed implies that the runoff path-lines are short and steep; in addition, rainfall is non-uniform in time and space. Large floods occur rapidly in the mid-to-downstream reaches of the watershed, causing substantial damage during the summer.

Fig. 1 Map of the Wu-Tu watershed and its four watershed divisions: (a) Tamshui River Basin; (b) Kee-Lung River watershed; and (c) Wu-Tu watershed divisions and their grids.

Data collection

The data used in this study were obtained from 14 raingauges located along the Tamshui River (Fig. 1(a)). Three of the raingauges (Jui-Fang, Wu-Tu and Huo-Shao-Liao) and one discharge site (Wu-Tu) are located within the Wu-Tu watershed (Fig. 1(b)). The study data comprised records for 110 rainfall–runoff events between 1966 and 2008. The degree of urbanization in the research area was gleaned from annual data on population density and imperviousness percentages. This study used data from 64 events to calibrate the model parameters of four watershed divisions by using relative urbanization variables, and employed the remaining events to test the model.

Two urbanization variables (impervious portion and population density) were used as the primary references of urbanization. All rainfall generates surface runoff in impervious areas. The annual imperviousness percentage for each year was obtained based on areas that include streets, roads, railroad lines, highways, roofs, buildings, parking lots, ponds, lakes and waterways. As shown in Figs 2 and 3, changes in population density and percentage of impervious area were plotted for four divisions of the Wu-Tu watershed between 1966 and 2008. The relationships between these urbanization variables and various other parameters were investigated to examine the hydrological consequences of the urbanization of the Wu-Tu watershed.

Fig. 2 Changes of population density in Wu-Tu watershed and its divisions.

Fig. 3 Changes of imperviousness percentage in Wu-Tu watershed and its divisions.

METHODS

This study applied a semi-distributed parallel-type hydrological model to describe the storm waters of the watershed and its four divisions. The suitability of potential models varies substantially, depending on the level of accuracy desired. Certain calculations, such as the calculation of mean rainfall and rainfall losses on the lumped watershed divisions, are necessary before rainfall–runoff routing can be determined. The univariate geostatistical method (block kriging technology) was used to estimate the areal rainfalls of the watershed and its divisions. The assumed spatially uniform losses of the watershed were calculated using the Ф-index method. The division inputs of the model for an event (i.e. the effective rainfalls of the relevant watershed divisions) represent the mathematical difference between the mean rainfall (per division) and the identical Ф-index computation. The optimal parameters for each division were calibrated using the semi-distributed model with an optimal seeking method and a reliable objective function. The parameters of the watershed divisions were used to observe the urbanization processes in each division.

Block kriging technique

Rainstorms vary greatly in space and time. Traditional methods, such as the Thiessen polygon method, are typically applied to calculate the areal rainfall based on data obtained from raingauge sites located in a watershed. These methods involve employing a linear combination of n raingauge observations Z(x_i) located at position x_i, and with weights, λ_i, which can be expressed as follows:

(1)

The block kriging method and its variants (Xie et al. 2011) also involve using equation (1) to obtain areal and point calculations to estimate the rainfall in a specific area. These methods have numerous applications in various research fields. Typical cases have included raingauge network design (Bastin et al. 1984, Cheng et al. 2008a), raingauge evaluation (Cheng 2011a, Cheng et al. 2012), spatial interpolation of rainfall (Goovaerts 2000, Syed et al. 2003, Basistha et al. 2008), and space-time rainfall interpolation (Cheng et al. 2007). The primary difference between the kriging method and traditional methods is the weighting computation. The kriging method entails using a semivariogram γ(t,h_ij) (Lebel and Bastin 1985)  to illustrate the spatial relationship of a rainfall process and describe rainfall variation in space, and can be expressed as follows:

(2)

where h_ij represents the distance between arbitrary raingauges x_i and x_j; T denotes the duration of all rainfall events; and p(t, x_i) defines the rainfall depth measured by the ith raingauge at the tth period.

The kriging system obtains the optimal weights using a specific semivariogram of rainfall. Thus, the system is obtained by applying Lagrange’s multiplier method, and the point or areal rainfall is calculated using the block kriging system:

(3)

(4)

where γ(x_i,x_j) is the semivariogram of raingauge x_i and raingauge x_j; represents the mean semivariogram of estimated area V and raingauge x_i; λ_j is the weighting of each raingauge; is the kriging estimated variance; and μ is the Lagrange multiplier.

Bastin et al. (1984) proposed a basic semivariogram for the quick computation of an hourly semivariogram. This model is referred to as the scaled climatological mean semivariogram, and is represented by . This approach can be applied to establish an hourly semivariogram γ(t,h_ij) by using dimensionless rainfall data obtained from a project basin. Equation (5) shows the relationship between the hourly semivariogram and the scaled climatological mean semivariogram:

(5)

where denotes the sill of the semivariogram for period t and is time-variant, a represents the range of the scaled climatological mean semivariogram and is time-invariant, and s(t) denotes the standard deviation of rainfall (mm) for all raingauges for period t. The basic semivariogram is expressed as:

(6)

A basic experimental semivariogram can be calculated using equation (6). This semivariogram is always derived from discontinuous point-observations; thus it is not spatially continuous. A realistic application of the block kriging method involves using a semivariogram model to obtain the spatial continuity of rainfall variations.

When applying the block kriging method to estimate the hourly mean rainfall during storm events within specific watershed divisions, the estimated area V in equation (3) must be divided into M grids. Therefore, equation (3) should be revised as follows:

(7)

where V_m is the mth grid of the estimated area V, and γ(V_m,x_i) represents the semivariogram of the mth grid of V_m and raingauge x_i.

Semi-distributed model of a parallel-type cascade of linear reservoirs

The model parameters are derived primarily by acquiring the results of simulated and observed rainfall–runoff events. Several conceptual models and their parameters have originated from the instantaneous unit hydrograph (IUH). This approach is referred to as UH-based modelling. Because UH-based conceptual models typically possess specific definitions for their parameters, any parameters with physical significance can be used to represent the hydrological status of urbanized watersheds during various periods. This study applied a semi-distributed model with a parallel connection to examine the effects of urbanization variables in watershed divisions that displayed a range of characteristics.

A watershed is frequently conceptualized as a multi-cascade reservoir. Reservoirs are typically used to represent rainfall–runoff transformations of a lumped watershed or combined watershed divisions. The general form of the IUH, U_n, from the nth cascade of linear reservoirs that has various storage constants k_n and period t can be derived as follows (Hsieh and Wang 1999):

(8)

A lumped watershed system having identical values for linear reservoir storage is a special case of equation (8). The IUH form of a cascade of n linear reservoirs for one SI (Système International) unit of effective rainfall (Nash 1957) is calculated as:

(9)

where Γ( · ) denotes the gamma function and .

This study considered each watershed division as an independent subsystem. The hydrological status of each subsystem was indicated by two significant parameters: K_o and K_c. The overland storage of a watershed division is represented by K_o, whereas K_c represents the channel storage. A connection from a watershed division to the watershed outlet was constructed using a specific flow path. The Wu-Tu watershed was divided into four subsystems; thus, four parallel flow paths were determined to route the outflow hydrograph at an outlet of the Wu-Tu watershed. The flow paths derived from the four divisions and their downstream channels at the watershed outlet were classified as follows:

Flow Path 1 from watershed division 1:

Flow Path 2 from watershed division 2:

Flow Path 3 from watershed division 3:

Flow Path 4 from watershed division 4:

where o_i denotes the overland storage of the ith watershed division, and c_i represents channel storage of the ith division.

The total direct runoff at the watershed outlet can be obtained by using the convolution integral, in which the spatially averaged rainfall excess I_i(τ) of each flow path is operated on by U_i(t–τ) and integrated over time t to yield outlet-runoff Q(t). For example, the IUH form of Flow Path 4 derived from equation (8) was determined as follows:

(10)

The convolution integral formula was obtained as follows:

(11)

where function U_i(t–τ) is the IUH derived from equation (8), the mean rainfall of the lumped divisions was obtained by using the block kriging method, the mean excess values were calculated by subtracting the mean rainfall from the corresponding Ф-index values by assuming spatially uniform losses, and N is the number of divisions in the watershed, set at N = 4 in this study.

CRITERIA FOR EVALUATING MODELS

To measure the suitability of the model parameters for the discussed basin, the following four criteria were used to analyse the goodness of fit.

Coefficient of efficiency

The coefficient of efficiency (CE) is given by:

(12)

where Q_est(t) denotes the discharge of the simulated hydrograph for period t (m³ s^-1), Q_obs(t) is the discharge of the observed hydrograph for period t (m³ s^-1), and represents the mean discharge of the observed hydrograph (m³ s^-1). A CE value close to 1 indicates a good fit.

Error of peak discharge (%)

The error of peak discharge (EQ_p) is given by:

(13)

where Q_p,est is the peak discharge of the simulated hydrograph (m³ s^-1), and Q_p,obs is the peak discharge of the observed hydrograph (m³ s^-1).

Error of time to peak

The error of the time for the peak () to arrive is defined as:

(14)

where T_p,est denotes the time (h) required for the peak to occur in the simulated hydrographs, and T_p,obs represents the actual time (h) required for the peak to occur in the observed hydrographs.

Error of runoff volume (%)

The error of runoff volume (VER) is given by:

(15)

RESULTS AND DISCUSSION

Urbanization variables are frequently used as an index of watershed development. Two examples of these variables are impervious percentage and population density. This study investigated the effects of changes in these two urbanization indices on hydrograph parameters among the four watershed divisions discussed. Identifying the relationships between the urbanization variables and model parameters enables the comparison of space–time continuous variations in the hydrological status of the various watershed divisions. This study examined data that represented parameters from the hydrological model across various periods. The proposed model is a semi-distributed model that can be applied to assess a parallel-type cascade of linear reservoirs. The model is a UH-based conceptual model of watershed divisions, with storage constants K_o and K_c representing separate storages of overland flow and channel flow, respectively.

Mean rainfall and effective rainfall in watershed divisions

A set of discontinuous point-rainfall depths (p(t,x)) during a specific period can be used to generate two-dimensional random fields. Considering n raingauges in a river basin, for each period t, a realization π(t) of the random n vector can be expressed as:

(16)

The hourly semivariogram of rainfall is a time function of period t, isotropy, and a time mean form with non-zero and T intervals. Rainfall measurements were obtained from 14 raingauges located along the Tamshui River between 1966 and 2008. The scaled climatological mean semivariogram and its power form (Fig. 4) applied for fitting were calculated as follows:

(17)

Fig. 4 Climatological mean semivariogram of the research watershed.

where ω₀ denotes the scaled parameter of the scaled climatological mean semivariogram (mm²).

The hourly variance s²(t) of each period t can be calculated using hourly measurements of rainfall that occurred in the identical realization π(t). Based on equations (5) and (17), each hourly semivariogram of rainfall is calculated directly from the hourly variance and scaled climatological mean semivariogram. Equation (7) was subsequently used to obtain the mean rainfall of the entire watershed and its four divisions, based on the hourly semivariograms, data from three raingauges located in the Wu-Tu watershed (Jui-Fang, Wu-Tu and Huo-Shao-Liao), and the grid divisions (Fig. 1(c)).

The lack of records for division-outlet flows implies that rainfall losses in each watershed division could not be measured directly using the principle of volume equilibrium. This principle states that the volume of effective rainfall on a watershed division is equal to that of the direct runoff from the division’s outlet. In other words, the value of rainfall loss is calculated as the difference between rainfall and effective rainfall (or direct runoff). To solve this problem, this study assumed that rainfall losses of rainfall–runoff events were spatially uniformly distributed throughout the entire watershed. These losses were represented by the Ф-index values, which were computed using the mean of measured rainfall and the separated direct runoff of events throughout the entire watershed. The hyetographs of an event over a divided area can be obtained by subtracting the Ф-index constant from the mean hourly rainfall for each division. The event-based effective rainfalls for the four watershed divisions were subsequently used as the inputs of the model for parameter calibration and verification in this study.

Calibration of representative parameters

The direct runoff of a watershed outlet results from a constant separation in baseflow. This study considered the baseflow as groundwater flow before the commencement of rainfall, or the lowest discharge of a rising limb in a streamflow hydrograph (if such discharge was considered a constant). The effective rainfall in the four watershed divisions was calculated using the Ф-index method, under the assumption of spatially uniform losses throughout the watershed.

The division parameters were obtained from 64 samples collected from 110 rainfall–runoff events between 1966 and 2008. The parameter optimization process adopted the shuffled complex evolution algorithm (Duan et al. 1993) to obtain values that were representative of each division’s parameters. The parameters that achieved significance were considered adequate to describe the urbanization status of the watershed and its divisions. Table 1 shows a comparison of simulated and observed runoff hydrographs for the four criteria (CE, EQ_p, ET_p and VER), and Fig. 5 shows four cases drawn from the 64 calibrated events.

Fig. 5 Calibration of observed and simulated hydrographs of typhoons and storms.

Table 1 Calibrated results of the selected events (dates given as dd-mm-yyyy) with the four evaluation criteria; CE: coefficient of efficiency; EQ_p: error of peak discharge; ET_p: error of time for peak to arrive; VER: error of runoff volume.

Event name (date)	Evaluation criteria				Event name (date)	Evaluation criteria
	CE	EQp (%)	ETp (h)	VER (%)		CE	EQp (%)	ETp (h)	VER (%)
Alice (2 September 1966)	0.85	–23.31	1	–4.62	Storm (4 November 2000)	0.93	–6.83	–1	–7.05
Jean (19 July 1974)	0.94	–18.67	0	–5.82	Bebinca (8 November 2000)	0.93	–6.74	1	1.98
Betty (22 September 1975)	0.87	–25.32	1	–4.37	Storm (13 December 2000)	0.97	–5.63	0	–2.87
Storm (16 September 1976)	0.87	–19.61	0	–4.13	Storm (19 December 2000)	0.92	–10.58	2	–5.17
Norris (27 August 1980)	0.87	–28.38	0	–5.53	Storm (7 June 2001)	0.73	–22.66	2	–7.97
Storm (19 November 1980)	0.93	–24.32	–1	–6.78	Storm (3 September 2001)	0.93	–12.81	0	–4.63
Storm (2 June 1984)	0.90	–30.47	0	–5.02	Nari (15 September 2001)	0.96	–15.03	–1	–3.45
Freda (7 August 1984)	0.90	–27.83	0	–6.19	Lekima (25 September 2001)	0.91	–24.91	3	–4.27
Storm (14 August 1984)	0.85	–33.14	1	–2.23	Mindulle (1 July 2004)	0.85	–11.80	–1	–6.27
Storm (8 February 1985)	0.96	–12.69	0	–0.17	Ranani (12 August 2004)	0.89	–29.72	–1	–14.11
Jeff (29 July 1985)	0.95	–20.07	–2	–0.06	Storm (9 September 2004)	0.85	–14.07	2	–8.45
Nelson (22 August 1985)	0.91	–29.04	0	–1.72	Storm (1 October 2004)	0.94	–16.81	0	–7.96
Brenda (3 October 1985)	0.91	–31.57	0	–2.16	Nanmadol (2 December 2004)	0.95	–13.48	0	–3.59
Wayne (22 August 1986)	0.81	–24.24	1	–6.38	Storm (7 December 2004)	0.98	–2.85	0	4.93
Alex (27 July 1987)	0.94	–22.31	0	–7.55	Storm (2 May 2005)	0.92	–25.72	–1	–9.71
Storm (24 September 1988)	0.93	–13.75	1	–1.94	Storm (9 May 2005)	0.94	–24.00	1	–2.44
Storm (28 July 1989)	0.90	–15.89	–3	–5.59	Storm (29 May 2005)	0.93	–21.39	1	–2.06
Sarah (11 September 1989)	0.93	–8.99	2	–1.84	Storm (10 September 2005)	0.86	–30.98	1	–9.31
Storm (1 September 1990)	0.96	–13.21	0	–5.11	Longwang (1 October 2005)	0.95	–20.98	1	–4.42
Storm (3 September 1990)	0.89	–26.87	1	–4.58	Storm (11 December 2005)	0.88	–22.69	4	4.29
Nat (22 September 1991)	0.94	–13.48	–2	–7.51	Storm (11 January 2006)	0.96	–16.16	1	7.23
Ruth (28 October 1991)	0.92	–17.34	1	–0.85	Storm (8 June 2006)	0.95	–17.67	0	–3.62
Polly (27 August 1992)	0.86	–21.83	0	–3.67	Storm (8 September 2006)	0.86	–28.63	0	–0.17
Storm (18 June 1994)	0.89	–24.76	1	–4.68	Storm (15 June 2007)	0.96	–22.52	0	–2.76
Doug (7 August 1994)	0.87	–33.42	2	–2.68	Wipha (17 September 2007)	0.95	–13.87	–4	–5.65
Fred (20 August 1994)	0.97	–13.66	1	–4.73	Krosa (6 October 2007)	0.98	–12.26	1	–5.84
Seth (19 October 1994)	0.92	–19.02	1	–5.22	Mitag (18 November 2007)	0.93	–25.10	1	0.94
Zeb (15 October 1998)	0.97	–3.52	0	–3.34	Storm (29 January 2008)	0.95	–13.30	1	2.21
Storm (23 April 2000)	0.94	–21.96	4	–1.01	Storm (4 February 2008)	0.85	–22.06	–2	1.55
Storm (12 June 2000)	0.94	–11.18	–2	–0.51	Storm (27 May 2008)	0.92	–19.29	1	–3.04
Kai-tak (8 July 2000)	0.90	–26.25	1	0.06	Storm (30 May 2008)	0.94	–10.35	1	3.18
Bilis (22 August 2000)	0.92	–24.37	0	–7.63	Storm (10 October 2008)	0.98	–12.64	1	–4.11

This study employed an identical set of four criteria to evaluate the suitability of the proposed model for the discussed watershed. CE is a tool for measuring the performance of a model in hydrology. This study also examined the differences between observed and simulated data for peak quantity and time to peak based on the EQ_p and ET_p, respectively. However, loss values could be calculated only for the entire watershed and not for its divisions, because of unavailable division flow data. Therefore, spatial uniformity of rainfall loss was assumed, and the VER was used to examine the volume equilibrium between the observed and simulated hydrographs of the watershed outlet.

In model calibration, 47 of the 64 rainfall–runoff events obtained CE values greater than 0.9, 16 events obtained values that ranged between 0.8 and 0.9, and only one event (7 June 2001) yielded a value below 0.8. For EQ_p, 24 cases showed less than 20%, and all remaining cases ranged between 20% and 30%. The five samples that produced approximately 30% were: two storms of 2 January 1984 and 14 August 1984; Typhoon Brenda (13 October 1985); Typhoon Doug (7 August 1994); and the storm of 10 September 2005. For ET_p, the values for all events were below 3 h, except for three events (two storms on 23 April 2000 and 11 December 2005, and Typhoon Wipha on 17 September 2007). Only one case (Typhoon Ranani, 12 August 2004) gave VER slightly larger than the VER criterion of 10%. These results show that the simulation values were similar to those that were observed, despite the assumption of spatially uniform losses. The model calibrations for the three previous criteria showed that the parameters of the semi-distributed model adequately represented the watershed conditions during urbanization.

Division parameters and corresponding urbanization variables

In this study, the outlet-hydrograph simulations were evaluated using four criteria. The model parameters that were based on the results of the hydrograph simulations may adequately represent the hydrological conditions of a watershed and its divisions. However, hydrograph parameters are affected not only by urbanization variables but also by weather factors, antecedent moisture conditions and other unknown variables. These unknown variables are referred to as hydrological uncertainties, which frequently result in irregular and unpredictable variations in parameters; thus, applying valid methods is necessary for examining obvious or visible tendencies. Previous methods that involve employing the annual average and optimal interval are used to solve this problem; applying such methods smoothes the disordered, unsystematic variations in division parameters across selected hourly-based cases. The use of the annual average and optimal interval methods has been shown to effectively smooth parameter behaviours for varying degrees of urbanization (Cheng and Wang 2002, Huang et al. 2008a).

The annual average method was used to obtain annual parameters for watershed divisions that displayed varying degrees of each urbanization case. An annual parameter was computed by obtaining the mean for various calibrated parameters for the same year. By contrast, the optimal interval method was used to estimate a suitable interval based on variations in the model parameters and urbanization variables. The values of these calibrated parameters in the same interval were considered a fixed value; thus, each interval had an identical computation. The optimal magnitude of the interval was finally determined until an obvious tendency between averaged values of model parameters and urbanization variables was observed. The optimal interval method is capable of providing greater clarity on the tendencies of parameters related to urbanization variables, although it is more complex than the annual average method (Huang et al. 2012). This study applied the optimal interval method for the correlation analysis to identify changes in parameters that are related to complex urbanization processes in the watershed divisions. Figure 6 shows the correlations between hydrological parameters and imperviousness percentages, and Fig. 7 shows the correlations between hydrological parameters and population density.

Fig. 6 Changes of parameter tendencies in relation to imperviousness percentages.

Fig. 7 Changes of parameter tendencies in relation to population densities.

Correlativity of parameter behaviour and urbanization variables

The development of an urban area within a basin is caused by an increase in population density, followed by urbanization. This process is inevitably accompanied by an increase in the number of impervious surfaces that prevent rainwater from being absorbed into the soil surface. The correlations of variables in an urbanizing watershed and its divisions can be analysed through the study of urban hydrology. From this perspective, the simulated results of outlet-hydrographs for the Wu-Tu watershed were completed for parameter calibration. The calibrated parameters indicated the possible effects of urbanization on the rainfall–runoff relationship (Figs 6 and 7). The parameters among the watershed divisions varied according to the degrees of urbanization, and changes in the calibrated parameters were related to urbanization variables, including imperviousness and population density.

Regression analysis was used to further examine the identified relationships between model parameters and urbanization variables in the watershed divisions. Before establishing correlations, the applicability of the urbanization variables and two division parameters was evaluated. This study concerns the applicability of urbanization variables for regular variations of hydrographic parameters. Previous research showed that changes related to urbanization in watershed divisions may be associated with imperviousness or population density (Cheng et al. 2010). Relationships may exist between model parameters and urbanization variables. Regarding hydrological effects, the results of this study show that imperviousness directly affected the model simulation because of the separation effect of rainwater, whereas population density was unrelated. As shown in Figs 2 and 3, the demands of population growth cause increased levels of imperviousness, whereas imperviousness decreased the uncertainly related to a population decrease. These results may be compared with the findings shown in Figs 6 and 7. The parameter changes related to imperviousness appeared to be more regular than those related to population density. Thus, the hydrological analyses conducted in this study confirm that imperviousness was directly related to the effect of urbanization.

In addition, the applicability of parameters having physical significance was related to the sensitivity of parameter variation for various changes throughout the urbanization process. The proposed parallel-type cascade of linear reservoirs model included two parameters for all watershed divisions: these two parameters represented the physical significance of overland and channel storages. Within a watershed or a division thereof, the area of overland storage is substantially larger than that of the channel storage. Regarding the area contributing to the runoff, the contribution to a watershed outlet of an overland storage must be larger than that of the channel storage, with a significant difference between overland and channel areas. Thus, parameter values representing the channel storage should vary with a greater level of consistency than those describing the overland storage, which is independent of urbanization changes. In response to imperviousness and population, K_o exhibited a greater degree of variance than K_c, as shown in Figs 6 and 7. Therefore, this study confirmed that the sensitivity to urbanization processes among the four watershed divisions was greater for the significance of overland storage (represented by K_o) than for the significance of channel storage (represented by K_c), regardless of imperviousness percentage or population density.

The two discussed applicability evaluations were completed previously. In the first type, the applicability of imperviousness as a direct variable affecting urbanization in hydrology was evaluated. It was shown that the overland storage parameter exhibited greater sensitivity than the channel storage parameter. This study used imperviousness as a major variable and assumed the channel parameter K_c as constant across each watershed division. The smoothed values of the channel parameters for each division were averaged independently to obtain the respective constants. The average results for the channel parameters in each watershed division are shown in Table 2. These fixed values denoted the channel storage of their corresponding divisions, and were independent of the urbanization process. Furthermore, the discrete values of overland division parameter K_o were related separately to the imperviousness percentage of each corresponding division. The natural logarithmic form of regression analysis was used to identify any salient relationships that may be suitable for further study (Cheng 2011b).

Table 2 Averaged constants of channel parameters (K_c) and regression results of overland parameters (K_o) related to corresponding imperviousness percentages (Im).

Division	Channel parameter, Kc	Overland parameter, Ko
1	2.947	ln(Ko1) = –0.085 · ln(Im) + 2.118; R^2 = 0.85
2	2.003	ln(Ko2) = –0.254 · ln(Im) + 1.442; R^2 = 0.66
3	2.665	ln(Ko3) = –0.386 · ln(Im) + 2.679; R^2 = 0.63
4	2.868	ln(Ko4) = –0.243 · ln(Im) + 3.514; R^2 = 0.69

The regression results for the overland parameters that were related to their corresponding imperviousness percentages are shown in Table 2 and Fig. 6. The coefficient of determination (R²) was used to evaluate whether any power relationships were favourable. The correlations shown in Table 2 and Fig. 6 provided convenient data to examine any changes in overland parameter K_o that were in response to continuous changes in impervious paving.

Relationship verification

The applicability of the urbanization variables and model parameters of two types has been evaluated previously. The imperviousness variable had a direct effect on the hydrograph modelling. Furthermore, K_o was more sensitive to changes in imperviousness than K_c. This study examined the data from 46 rainfall–runoff events that occurred between 2002 and 2008 to further examine reported correlations between K_o and imperviousness percentage.

The baseflow of an event was assumed to be a constant that was derived either from the groundwater flow before the start of rainfall, or from the lowest discharge of a rising limb in a streamflow hydrograph. The mean rainfall of the entire watershed and its divisions were estimated using the block kriging method. This estimation was also based on the assumption in calibration that rainfall loss was spatially uniformly distributed throughout the watershed. The assumed-uniform Ф-index of an hourly-based event was computed according to the entire mean rainfall and outlet direct runoff. For each division, the uniform loss value was subtracted from the mean rainfall to yield the effective rainfall of that division, which allowed the rainfall–runoff model to be verified. The second column of Table 2 shows that the values of parameter K_c were fixed division constants; the K_o values obtained through regression equations of watershed divisions are also shown. Parameter K_o exhibited non-linear variation with respect to changes in annual imperviousness percentages. The effective rainfall constants K_c and varying parameters K_o were subsequently used to obtain direct runoff estimations. The simulated runoffs were compared with actual runoff observations using the four evaluation criteria.

Table 3 shows the results for the comparison among 46 cases, and Fig. 8 shows plots of four verified cases among the 46 rainfall–runoff events. The CE values for the model verification exceeded 0.70 in all but five events, the EQ_p values were less than below 25% in all but nine events. The ET_p values for all but six events were below 3 h, and the VER values of all examined events were less than 10%.

Fig. 8 Verification between observations and simulations of typhoons and storms.

Table 3 Verification results of the selected events with the four evaluation criteria.

Event name (date)	Evaluation criteria				Event name (date)	Evaluation criteria
	CE	EQp (%)	ETp (h)	VER (%)		CE	EQp (%)	ETp (h)	VER (%)
Gilda (17 November 1967)	0.86	–36.78	1	–1.45	Storm (5 September 2001)	0.84	–10.40	3	1.95
Storm (25 July 1968)	0.72	–10.51	2	–1.95	Rammasun (3 July 2002)	0.90	–6.19	0	–1.59
Elaine (28 September 1968)	0.77	–16.53	1	–0.10	Storm (26 March 2004)	0.87	19.08	0	–1.25
Elsie (26 September 1969)	0.83	–25.57	3	–2.79	Storm (8 July 2004)	0.85	–20.03	1	–1.65
Fran (5 September 1970)	0.89	–24.19	0	–1.31	Storm (18 October 2004)	0.79	–22.96	0	2.42
Agnes (17 September 1971)	0.84	–16.15	2	–0.75	Storm (23 December 2004)	0.77	26.22	1	2.00
Bess (22 September 1971)	0.70	–38.25	3	–1.80	Storm (24 February 2005)	0.69	–3.67	1	–0.12
Storm (22 October 1974)	0.82	–38.81	–1	–9.41	Haitang (17 July 2005)	0.79	12.76	4	–0.32
Storm (3 July 1976)	0.74	23.40	1	–1.01	Sanvu (11 August 2005)	0.78	11.21	4	–8.70
Ora (12 October 1978)	0.79	–38.50	4	–5.96	Talim (30 August 2005)	0.85	20.98	0	0.28
Cecil (10 August 1982)	0.64	–36.44	2	–7.17	Storm (6 October 2005)	0.88	6.65	2	–0.35
Storm (18 November 1984)	0.82	–29.86	2	–3.94	Storm (20 January 2006)	0.93	–6.46	0	2.75
Gerald (9 September 1987)	0.86	–7.02	0	–3.37	Chanchu (13 May 2006)	0.87	–6.41	3	0.84
Lynn (23 October 1987)	0.81	–16.82	–2	–2.13	Storm (6 June 2006)	0.65	2.43	2	–1.68
Storm (16 September 1988)	0.74	–1.63	2	–3.13	Storm (13 December 2006)	0.87	–9.30	0	1.58
Storm (28 September 1988)	0.91	–13.58	4	–0.32	Storm (18 December 2006)	0.74	–16.57	0	2.87
Storm (5 June 1993)	0.84	13.87	3	–0.06	Storm (26 June 2007)	0.66	22.60	1	–0.45
Gladys (1 September 1994)	0.91	–23.64	1	–0.87	Storm (4 September 2007)	0.87	–24.35	2	–4.49
Zane (27 September 1996)	0.88	–22.97	2	–1.69	Storm (21 January 2008)	0.74	14.97	–4	1.24
Storm (26 November 1998)	0.87	9.91	2	–0.62	Storm (17 February 2008)	0.78	5.73	–1	2.60
Dan (3 October 1999)	0.77	22.25	1	–2.13	Jangmi (28 September 2008)	0.81	–2.65	4	–3.16
Storm (17 June 2000)	0.84	0.00	–1	–5.10	Storm (6 October 2008)	0.89	–15.30	2	2.60
Storm (16 November 2000)	0.64	29.84	2	7.34	Storm (8 December 2008)	0.81	–25.68	0	–3.30

Regarding the two discussed model parameters, the variation shown by K_o exhibited a greater level of sensitivity to imperviousness than K_c. Furthermore, the R² values clearly showed a non-linear correlation between K_o and imperviousness. The verification analysis (based on 46 rainfall–runoff events) also showed a clear correlation between changes in K_o and changes in imperviousness. These results confirm that the equations listed in Table 2 and Fig. 6 reflect the hydrological effect of urbanization. Thus, the impervious area is a major variable indicating the ongoing urbanization of the Wu-Tu watershed. Furthermore, this analysis can be applied further in relevant studies.

CONCLUSIONS

The applicability of urbanization variables and parameter types on watershed divisions were successfully evaluated based on the 64 calibrated and 46 verified events, as well as the integrated approaches. These methods were used to establish the correlation between hydrograph parameters and urbanization variables among four watershed divisions. The approaches employed in this study did not require detailed hydrological data.

The block kriging method was used to estimate the mean rainfall across the watershed and its divisions. Rainfall losses were assumed to have spatially uniform Ф-index values, whereas the distributed effective rainfall within divisions was input to the semi-distributed model of a parallel-type cascade of linear reservoirs. The results for the error of runoff volume criterion (VER) VER show that the principal volume equilibrium was unaffected by the assumption of spatially uniform loss distribution because the simulation volumes were close to those that were actually measured. The remaining three evaluation criteria (CE, EQ_p and ET_p) also confirmed that simulations of the watershed-outlet hydrographs supported the simulation performance well and that the calibrated parameters were reliable for relating them with urbanization variables. Changes in these division parameters resulting from the optimal interval method reflected various degrees of urbanization variables among the divisions.

This study showed that the discussed model parameters consistently exhibited a greater level of sensitivity to changes in imperviousness than to changes in population density. The parameter representing overland storage showed a marked response to changes in imperviousness, whereas the parameter denoting channel storage is independent of the urbanization process. Based on these two confirmations, the averaged values, which were derived from the optimal interval method of the overland parameter, were related to the imperviousness percentages, and the channel parameters were constant within each watershed division. The power linkage provides a method for linking continuous relationships between division parameters and their corresponding imperviousness levels. The verification results for the four discussed evaluation criteria further confirm that the power relationship between the overland parameter and imperviousness responded appropriately to changes in urbanization throughout the watershed divisions. These power equations yielded parameter values that are representative of overland storages for each watershed division based on the observed imperviousness data. Changes in hydrograph characteristics among the watershed divisions were identified to represent various degrees of urbanization among these divisions.