Long-term effects of land use/land cover change on surface runoff in urban areas of Beijing, China

Abstract The objective of this paper is to present a case study to derive land use/land cover (LULC) maps and investigate the long-term effects of LULC change on surface runoff in the fast urbanizing Beijing city. The LULC maps were derived from Landsat TM/ETM+ imagery (acquired in 1992, 1999, 2006, and 2009) using support vector machine method. A long-term hydrologic impact assessment model was applied to assess the impact of LULC change on surface runoff. Results indicated that the selected study area experienced rapid urbanization from 1992 to 2009. Because of urbanization, from 1992 to 2009, modeled runoff increased 30% for the whole area and 35% for the urban portion. Our results also indicated that the runoff increase was highly correlated with urban expansion. A strong relationship ( R 2 = 0.849 ) was observed between the impervious surface percent and the modeled runoff depth in the study area. In addition, a strong positive relationship was observed between runoff increase and percentage of urban areas ( R 2 = 0.997 for the whole area and R 2 = 0.930 for the urban portion). This research can provide a simple method for policy makers to assess potential hydrological impacts of future urban planning and development activities.

Long-term effects of land use/land cover change on surface runoff in urban areas of Beijing, China

Introduction
Land use/land cover (LULC) changes have direct impacts on the hydrological cycle and stream quality. There are also indirect impacts on climate and the subsequent impact of the altered climate on the waters. 1 Therefore, LULC changes have been treated as one of the most important and sensitive factors for global environmental change. Urbanization is a major force that is driving LULC changes. Although urbanization due to urban sprawl provides social and economic benefits to the community, the detrimental consequences of urbanization to the urban environment are widespread, especially in most developing countries. Increased impervious surface percentage in urbanized areas can reduce the time of runoff and increase the peak discharges of stream flow, resulting in large and more frequent incidents of flooding. 1,2 In addition, because of the pollutants in runoff and sediments, the increase of impervious surfaces also degrades stream quality. 2,3 Therefore, urbanized areas are increasingly susceptible to flooding under high precipitation conditions. 4 In recent decades, with a constant increase in urban population, urbanization has become a significant urban environmental and ecological concern, especially in most developing countries. Therefore, it is of great significance to accurately map LULC changes in urbanized areas and to evaluate its impacts on surface runoff and stream quality for urban planning policy makers and water/land resource management. Due to their relatively low cost and suitability for large area mapping, satellite remote sensing (RS) images have been widely applied for detecting urban growth and accurate and timely mapping of LULC changes over a long period. 5 Furthermore, based on the RS, LULC changes mapping also allows for a rapid assessment of its impacts on surface runoff and water quality. 2,3 Recently, numerous researchers have assessed the impact of LULC changes induced by urbanization on surface runoff. [6][7][8][9] Ng and Marsalek 10 used hydrologic simulation program fortran (HSPF) model to simulate the effects of urbanization on stream flow for the Waterford River Basin near St. John's, Newfoundland, Canada. Kang et al. 11 investigated the effect of urbanization on runoff from the On-Cheon Stream watershed in Pusan, Korea, using the linear reservoir model; results revealed that the peak discharge increased and the mean lag time decreased owing to urbanization. Brun and Band 12 also used HSPF to assess the effects of land use change on urbanizing watershed behavior. Brun and Band determined a runoff ratio relationship indicating the existence of a threshold percent of impervious cover (∼20%) above which the runoff ratio changes dramatically. Weng 1 used the soil conservation service (SCS) model for estimating the effects of urban growth on surface runoff with the integration of RS and geographic information system (GIS). The results revealed that urbanization lowered potential maximum storage and thus increased runoff coefficient values. Shuster et al. 13 presented a literature review about impacts of impervious surface on watershed hydrology and found that impervious surfaces have a major effect on watershed hydrology. Shi et al. 4 also used SCS model to assess the effect of LULC change on surface runoff in Shenzhen region, China. The results showed that urbanization intensified the flood process. Li and Wang 2 used a long-term hydrologic impact assessment (L-THIA) model to evaluate the effect of LULC change on surface runoff in Dardenne Creek watershed, St. Louis, Missouri. Results indicated that the runoff increase is highly correlated with urban expansion. Franczyk and Chang 14 used the ArcView Soil and Water Assessment Tool hydrological model to evaluate the impacts of climate change and urbanization on the runoff in a fast urbanizing region. Hamdi et al. 15 used an urban soil-vegetation-atmosphere-transfer model to examine how the surface runoff of the Brussels Capital Region responded to historical urbanization (1960 to 1999).
To simulate the land use change impacts on watersheds, many hydrologic and water quality models have been developed and integrated with a GIS, such as the United States Department of Agriculture (USDA) Technical Release 55 (TR-55), 16 EPA Storm Water Management Model, 17 and Soil and Water Assessment Tool. 18 Although some of these models can be used to assess the long-term impacts of LULC changes on hydrology and stream quality, they require extensive data inputs and parameters that are usually not available to land-use planners. 2 Therefore, there has been a need for a much simpler-to-use model to assess relative hydrological impacts of LULC change at the watershed scale. The L-THIA model was developed and integrated with GIS to estimate direct runoff from very basic input data, such as daily rainfall, land use, and the hydrologic soil group. 19,20 In China, rapid urban sprawl has had a profound influence on runoff and resulted in larger and more frequent incidents of flooding in many urban areas. Over the past three decades, Beijing has been one of the fastest-growing urban areas in China and has undergone intense urbanization that has significantly impacted the urban surface runoff. The objective of this paper is to derive LULC classification maps and further quantitatively assess the effects of LULC change on surface runoff in Beijing. The LULC maps are derived from Landsat TM/ ETM+ imagery (acquired in 1992, 1999, 2006, and 2009) using a support vector machine (SVM) method. To estimate the effects of LULC change on runoff, the L-THIA model is used. The rest of this paper is organized as follows. Section 2 describes the study area and the datasets used. Section 3 discusses methodology. The results and analysis are presented in Sec. 4. Section 5 presents the conclusions.

Description of Study Area
Beijing, the capital of China, was selected as the study area to evaluate the impact of urbanization on surface runoff (Fig. 1). It is located within 39.28°N to 41.05°N and 115.25°E to 117.35°E, covering ∼16;386.30 km 2 . There are five large water systems, including Juma River, Rongding River, Beiyun River, Chaobai River, and Jiyun River. With a warm temperate-zone semihumid continental monsoon climate, it is hot and rainy in the summer, cold and dry in the winter, and short seasons in the spring and autumn. Average annual temperature is 10 to 12°C, the temperature in January is from −7 to −4°C, and the temperature in July is 25 to 26°C. The average annual precipitation is 600 mm, and rainfall in the summer season from June to August comprises 75% of the annual total. According to the soil horizontal zonality, the main soil types are cinnamon soil, moisture soil, and brown soil. The study area contains a number of LULC types, including central business district (CBD), high-density residential, low-density residential, agricultural land, forest, exposed soil, and water bodies. Forest covers 50% of the land area. Over the past three decades, with economic development, Beijing has been one of the fastest-growing urban areas in China and has undergone intense urbanization. Beijing's development pattern is a typical concentric expansion, showing a ring-shaped pattern as you move from the inner city to the outskirts.

Datasets and Preprocessing
Four Landsat TM/ETM+ scenes, 13 WorldView scenes, soil-type classification data, daily precipitation data, and daily runoff data were used in this research. Their major characteristics are summarized in Table 1.

Remotely sensing imagery
LULC data were derived from the multispectral Landsat TM/ETM+ imagery (acquired in 1992, 1999, 2006, and 2009). These images were rectified to a Universal Transverse Mercator (UTM) coordinate system, with the pixel size of 30 m. Thirteen WorldView scenes (acquired on February 11, 2009) with a spatial resolution of 0.5 m were registered to the TM image (UTM) and used for LULC accuracy assessment. The high-resolution WorldView scenes did not contain snow, ice, or cloud.
The digital number of the TM scenes was converted to unitless planetary reflectance following the method proposed by Chander and Markham. 21 We assumed homogeneous atmospheric conditions within the scene, so atmospheric corrections were not performed.

Soil data
The soil-type classification map (1:200,000) of the study area was obtained from the Soil and Fertilizer Workstations of Beijing. Curve number (CN) values were used to determine soil hydrological input data. According to four hydrologic soil groups defined by the U.S. Soil Conservation Service (SCS), 16 the soil-type classification data were reclassified to output the soil hydrological map using ARCVIEW software [ Fig. 2(a)]. For the purpose of overlaying with same pixel data, the soil hydrological map was rectified to a UTM coordinate system, with a spatial resolution of 30 m.

Hydrological data
There are two weather stations and one gauge station in the selected area ( Fig. 1). Thirty-two years (1978 to 2009) of daily precipitation data were obtained from the meteorological department [ Fig. 2(b)]. In addition, the observed daily runoff data were obtained from the Beijing gauge station. These data were used for the L-THIA model parameter calibration.

Methodology
With the development of RS, GIS, and hydrological modeling, it is possible to investigate longterm impacts of LULC change on urban surface runoff. It is difficult for the distributed hydrological models to assess the long-term impacts of LULC change on hydrology because they require extensive data inputs and parameters (e.g., meteorological data, hydrological data, geological data, terrain data, soil data, and so on) that are usually not available for land-use planners. In our research, the simpler-to-use L-THIA model integrated with GIS and RS was used to assess the impacts of LULC changes on surface runoff in an urbanizing area. The L-THIA model requires very basic input data, such as daily precipitation data, LULC, and hydrologic soil group data. LULC classification maps were derived from long-time-series Landsat TM/ETM + imagery. Figure 3 shows the data processing framework of research on hydrological response to LULC changes. In Fig. 3, the left dashed box is the data processing flow chart of LULC change, while the right dashed box is the processing flow chart of L-THIA model simulation. In our research, the robust algorithm, SVM, was selected to perform a supervised classification and further to evaluate the temporal change of LULC in the study area. In the following section, a description of the SVM algorithm and L-THIA model is provided.

Support Vector Machine
The concept of the SVM originated from a nonparametric machine learning methodology based on Vapnik's 22 structural risk minimization (SRM) principle. The basic idea of SVM is to map data into a high-dimensional space and to find the hyperplane of the different classes with the maximum margin between them. SVMs have been used for image classification, [23][24][25][26] soil moisture estimation, 27 image retrieval, 28 and impervious surface estimation 29,30 for remote sensing studies in recent years. The theory of SVM has been extensively described in the literature. 31,32 Therefore, only a brief description of the concept of SVM in the framework of classification is given here. Suppose that training data ðX i ; y i Þ (X i ∈ R K belongs to a K-dimensional space, i ¼ 1; : : : ; L, L is the number of the training samples), where (y i ∈ fþ1; −1g) is the class label of the training vector X i for a two-class problem and can be linearly separated by a hyperplane [Eq. (1)].
where W is a weight vector and b is a scalar, often referred to as a bias.
To describe the separating hyperplane, let us use the following form: The separating hyperplane that creates the maximum margin is called the optimal separating hyperplane (OSH). The goal of the learning process based on SVM is to find the OSH to separate the training data by solving the following quadratic optimization problem: subject to the following constraints: : : : ; L ; (4) Fig. 3 The data processing framework of research on hydrological response to land use/land cover (LULC) changes.
where C denotes a penalty parameter on the training error and ξ i is called the slack variable. By converting the problem with Karush-Kuhn-Tucker conditions into an equivalent Lagrange dual problem, the optimization problem becomes Subject to∶ The final decision function is given by In Eqs. (6) and (7), α i denotes Lagrange multipliers; the parameter N (usually N ≪ L) stands for the number of the selected points or support vectors; and the symbol KðX; X i Þ is the kernel function that measures nonlinear dependence between the two input variables X and X i . Three admissible kernel functions include the polynomial kernel function, the radial basis kernel function (RBF), and the sigmoid kernel function.
In the theoretical development, the SVM has been developed for a two-class problem. Strategies are needed to adapt this method to multiclass cases. For example, one against all or one against one (OAO) algorithms have been proposed to adapt the SVM to multiclass problems. 33 However, support vector classification predicts only class labels but not probability information. In our research, percent impervious surface was obtained from the RS imagery using an SVM method. Therefore, it is essential to briefly describe how we extend SVM for probability estimates. More details of SVM have been extensively described in the literature. 34,35 Given k classes of the data, for any X, the goal is to estimate p i ¼ pðy ¼ ijXÞ; i¼ 1; 2; · · · ; k: Following the setting of the OAO approach for multiclass classification, the pairwise class probabilities can be estimated with the following equation: 35 where A and B are estimated by minimizing the negative log-likelihood function. Then, the probabilities p i can be derived from all r ij 's by solving the following optimization problem: 34 subject to the following constraints: As described above, the steps in SVM modeling are (1) selecting a suitable kernel function and kernel parameter (kernel width G) and (2) specifying the penalty parameter C. In our paper, an RBF kernel is used for constructing the classifier. In addition, specifying parameters G and C is the key step in the SVM because their combined values determine the boundary complexity and thus the classification performance. 36 In our research, the kernel width G was set to 0.15 and the penalty parameter was set to 100. 30 For the implementation of the training and modeling procedure, we employed the existing SVM library LIBSVM presented by Chang and Lin. 35 The LIBSVM represents integrated software for SVM classification, regression, distribution estimation, and multiclass classification. In order to keep the classification runtime low, the maximum number of training and testing samples for a given class did not exceed 3,000 and 5,000, respectively.

L-THIA Model Description
The L-THIA model was developed to estimate direct runoff using the CN method from daily rainfall depth, LULC, and hydrologic soil group data. 19 It provides a quick and simple impact analysis for initial plan review and for community awareness of potential long-term problems. 37 For runoff calculation, L-THIA uses the CN method developed by the USDA SCS (now Natural Resources Conservation Service), which is a core component of many hydrology models. 38 The runoff equation is where Q is the runoff (inches), P is the rainfall (inches), S is the potential maximum retention after runoff begins (inches), and I a is the initial abstraction (inches), which can be approximately estimated by Using Eq. (13), Eq. (12) can be rewritten as S, in mm, is related to the soil and LULC conditions and can be described by CN through  where the CN value is the main parameter of the L-THIA model and is used to describe the relationship between rainfall and runoff. CN is determined through the combination of LULC and soil (hydrological soil group) information for each cell. Theoretically, CN values are between 0 and 100. However, practically, CN values change between 30 and 100. Then, L-THIA model uses daily rainfall series and the CNs to calculate the mean annual total surface runoff in terms of depth and volume. 38 Figure 4 shows the data processing flow chart of the L-THIA model. 20,38 As illustrated in Eq. (14), the L-THIA model requires long-term daily precipitation data, soils, and LULC data. Daily precipitation series for the period of 1978 to 2009 were obtained from the meteorological department. Soils GIS data are required as a hydrological soilgroup. The soil data were obtained from the Soil and Fertilizer Workstations of Beijing. The LULC data were derived from Landsat TM/ETM+ imagery (acquired in 1992, 1999, 2006, and 2009) using SVM method.

LULC Change
In our research, we grouped the LULC types into seven categories: (1) high-impervious surface (>90%, high-density urban, including CBD, roads, squares, and so on), (2) medium-impervious surface (50 to 90%, high-density residential land), (3) low-impervious surface (<50%, low-density residential land), (4) arable land (agricultural land), (5) forest land, (6) bare land, and (7) water body. Because we focused on urbanization, only urban-related classes, including high-/medium-/low-impervious surfaces, were utilized for accuracy assessment. Sun et al. proved that the estimation accuracy of impervious surfaces was superior to Landsat TM/ ETM+ imagery (acquired on June 3, 2009) using an SVM method. 30 Figures 5 and 6 illustrated the percentages and spatial distributions of LULC classes in the Beijing region from 1992 to 2009. Forest and arable lands dominated in the whole area but decreased with the continuous increase of urban areas, especially agricultural land. The urbanization level (high-impervious surface), the proportion of urbanized area to the total land area, was 4.28, 6.03, 8.34, and 12.78% in 1992, 1999, 2006, and 2009, respectively. The proportion of residential area to the total land area was 3.21, 3.15, 3.83, and 6.66% in 1992, 1999, 2006, and 2009, respectively. Urban areas increased from 4.28% in 1992 to 12.78% in 2009, and residential areas increased from 3.21% in 1992 to 6.66% in 2009, while agricultural land areas decreased from 27.20% in 1992 to 11.44% in 2009. From Fig. 6, we can observe that Beijing's development pattern is a typical concentric expansion, showing a ring-shaped pattern. In the LULC maps of 1992 [ Fig. 6(a) In summary, the selected areas experienced rapid LULC change, especially urbanization, from 1992 to 2009. Temporally, the urbanization rate was relatively low in the beginning of 1990s and accelerated after 1999. Fig. 9 The statistical relationship between the impervious surface percent and runoff depth.

Parameter Determination of the L-THIA Model
According to the L-THIA model, the CN value is determined by the combination of LULC types and soil (hydrological soil group) information for each cell. The LULC types were derived from Landsat TM/ETM+ imagery as described earlier. According to soil classification categories in the L-THIA model, the soil data obtained from the Soil and Fertilizer Workstations of Beijing were reclassified to create the soil hydrological map using ARCVIEW software. Based on CN values of the L-THIA model calibrated by the observed daily runoff data, local soil infiltration properties, and other published results, CN values in the Beijing region were estimated and are presented in Table 2. Figure 7 shows the CN value distribution of the study area in 1992, 1999, 2006, and 2009. Land surface infiltration categories: 16 (a) high infiltration rate, (b) moderate infiltration rate, (c) low infiltration rate, and (d) very low infiltration rate.

Modeled Runoff Variations
In our research, the L-THIA model was used to simulate annual average runoff from daily rainfall data, soil data, and LULC data derived from long-term TM/ETM+ imagery (1992, 1999, 2006, and 2009). Figure 8 shows the simulated annual average runoff using the L-THIA model for 1992, 1999, 2006, and 2009, respectively. From the figure, we can observe that urban area corresponds to the highest runoff depth and bare land tends to higher runoff depth, while forest land corresponds to the lowest runoff depth followed by agricultural land. With the expansion of urban areas, the runoff depth significantly corresponded to urban area increase from 1992 to 2009. In order to analyze the relationship between impervious surface percent (ISP) and runoff depth, 4,000 sample points were randomly selected from the ISP image and modeled runoff depth acquired in 2009 (Fig. 9). Figure 9(a) shows the ISP derived from Landsat TM/ETM+ using SVM method; Fig. 9(b) shows the simulated annual average runoff using the L-THIA model; and Fig. 9(c) illustrates the statistical relationship between the ISP and runoff depth. From the figures, a strong relationship (R 2 ¼ 0.849) was observed between the ISP and the modeled runoff depth in the study area.
To further evaluate the effects of urbanization on surface runoff, using the simulated runoff for 1992 as the reference, a runoff change ratio γ for a certain year (such as 1999) can be calculated using the runoff simulated for that year divided by the runoff simulated for 1992. This runoff change ratio can be used to assess the spatial distribution of runoff variations due to the LULC change. 2 Then, the ratios were divided into four categories: (1) no runoff increase (γ ≤ 1.0); (2) small runoff increase (1.0 < γ ≤ 2.0); (3) moderate runoff increase (2.0 < γ ≤ 3.5); and (4) significant runoff increase (γ > 3.5). Figure 10 shows the spatial distribution of modeled runoff change in 1999, 2006, and 2009 compared to the reference in 1992. Because of the temporal difference, bare land could not be differentiated from the vegetation area in the optical imagery. Bare land corresponded to higher runoff depth, while vegetation corresponded to lower runoff depth. Therefore, LULC changes (bare soil and forest land) resulted in significant runoff changes in the north of the study area [Figs. 10(a) and 10(c)]. In Fig. 10(b), because of short temporal difference, LULC change was small in vegetated areas. Therefore, runoff depth change was very small, and no runoff change and small runoff change were prevalent in forested land. The temporal difference had no influence on the urban area, so we could observe the following from the urban portion (as shown in Fig. 10): (1) Compared to the modeled runoff for 1992, spatial variation of runoff change in 1999 was limited; however, a moderate-to-significant runoff increase became apparent around the urban area [ Fig. 10(a)]. (2) In 2006, with the expansion of the urban area, runoff change significantly increased around the city, and a moderate-to-significant runoff increase dominated gradually in the urban portion [ Fig. 10 , we can see that runoff depth increases with urban expansion. Therefore, the runoff increase was highly correlated with urban expansion. Figure 11(a) expressed the modeled runoff increase percentage for 1999, 2006, and 2009 compared to 1992. In this figure, modeled runoff increased 30% from 1992 to 2009, while modeled runoff increased 35% for the urban portion (Fig. 1). The runoff increase was relatively low (<5%) from 1992 to 1999 and much higher (>20%) after 1999. Figure 11(b) illustrates the relationship between the runoff increase and the percentage of urban area. From the figure, we can see a strong positive relationship. The whole area produced a correlation coefficient R 2 ¼ 0.997, while the urban portion produced a correlation coefficient R 2 ¼ 0.930. Our results indicated that the runoff increase was highly correlated with urban expansion.

Conclusions
This paper presented a case study to investigate the long-term impacts of LULC change on surface runoff in the fast urbanizing Beijing city, China. First, the LULC maps were derived from Landsat TM/ETM+ imagery (acquired in 1992, 1999, 2006, and 2009) using an SVM method. Then, a simpler-to-use L-THIA model was applied to simulate surface runoff variations in the study area from 1992 to 2009. Finally, the long-term impacts of urbanization on surface runoff were assessed using LULC classification maps and modeled surface runoff variations.
Our results indicated that the study area experienced rapid urbanization from 1992 to 2009. Urban areas increased from 4.28% in 1992 to 12.78% in 2009 and residential areas increased from 3.21% in 1992 to 6.66% in 2009, while agricultural land areas decreased from 27.20% in 1992 to 11.44% in 2009. To analyze the relationship between ISP and runoff depth, 4000 sample points were randomly selected from the ISP image and modeled runoff depth acquired in 2009. A strong relationship (R 2 ¼ 0.849) was observed between the ISP and the modeled runoff depth in the study area. As a direct result of the urbanization from 1992 to 2009, the long-term surface runoff increased 30% for the whole area, while modeled runoff increased 35% for the urban portion. Our results also indicated a strong positive relationship between runoff increase and the percentage of urban areas. For the whole area, the correlation coefficient R 2 is 0.997, while for the urban portion, the correlation coefficient R 2 is 0.930. Therefore, the runoff increase was highly correlated with urban expansion. This research provides a simple method for policy makers to assess potential hydrological impacts of future planning and development activities.