*Article* **Spatial Pattern of Functional Urban Land Conversion and Expansion under Rapid Urbanization: A Case Study of Changchun, China**

**Guolei Zhou 1,2,\*, Jing Zhang 1,2, Chenggu Li 1,2 and Yanjun Liu 1,2**


**Abstract:** As populations continue to be concentrated in cities, the world will become entirely urbanized, and urban space is undergoing a drastic evolution. Understanding the spatial pattern of conversion and expansion of functional urban land, in the context of rapid urbanization, helps us to grasp the trajectories of urban spatial evolution in greater depth from a theoretical and practical level. Using the ESRI ArcGIS 9.3 software platform, methods, such as overlay analysis, transition matrix, and kernel density estimation, were used in order to analyze the spatiotemporal characteristics of different types of functional urban land conversion and expansion in the central city of Changchun. The results show that different types of functional urban land were often expanded and replaced, and the urban spatial structure was constantly evolving. The conversion and expansion of functional urban land show similar characteristics to concentric zone and sector modes and show dynamic changes in different concentric circles and directions at different periods. Our method can accurately identify the different types of functional urban land, and also explore the evolutionary trajectory of urban spatial structure. This study will help to coordinate the development of different functional urban spaces and to optimize the urban spatial structure in the future.

**Keywords:** functional urban land; urban space; urban land use/cover change; urbanization; Changchun

#### **1. Introduction**

As the urban population continues to increase, cities have always been the focus of scholars and policy makers [1]. More than 50% of the world's population is currently concentrated in urban areas [2] and it is estimated that this will reach 70% by the year 2050 [3]. With the rapid increase in urban population, urban space is constantly expanding [4]. Urban–rural population migration [5], cross-regional population mobility [6], and bidirectional urban flows [7], are considered to be the main driving forces of the increase in the urban population. Changes in the urban land use/cover and their environmental consequences, such as the sustainable development of urban space [8] and the spatial mode of urban land change [9], have become important topics in the field of urban land research. In the context of rapid urbanization, the phenomenon of conversion between different types of functional urban land often occurs in urban built-up areas [10]. However, scholars have not paid enough attention to the phenomenon. Although big data, such as point-of-interest (POI) and vehicle trajectories, have been used to identify functional urban areas, it is difficult to accurately classify the different types of urban land [11–13]. Relying on the urban land use maps of the past few years and other auxiliary data, we have obtained a detailed classification of urban construction land based on different uses, which is helpful for analyzing spatial patterns of functional urban land conversion and expansion.

**Citation:** Zhou, G.; Zhang, J.; Li, C.; Liu, Y. Spatial Pattern of Functional Urban Land Conversion and Expansion under Rapid Urbanization: A Case Study of Changchun, China. *Land* **2022**, *11*, 119. https://doi.org/ 10.3390/land11010119

Academic Editors: Víctor Hugo González-Jaramillo, Antonio Novelli and Piyush Tiwari

Received: 24 November 2021 Accepted: 10 January 2022 Published: 12 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Unlike urban land conversion, the expansion of urban land has become a relatively common phenomenon worldwide and has attracted widespread research attention [14]. Remote sensing and geographic information systems (GIS) facilitate the study of urban land expansion [15,16]. Xiao et al. (2006) used Landsat images and the annual urban growth rate to analyze land use change and urban spatial expansion in Shijiazhuang [17]. López et al. (2001) analyzed and predicted land use change on the fringes of a Mexican city [18]. Polimeni (2005) analyzed the replacement of agricultural land by residential land caused by the expansion of urban land in the Hudson River Valley [19]. The expansion of urban land occupies a large part of non-urban land [20], especially in the metropolitan areas of Southeast Asia [21]. For example, a large area of high-quality farmland has been replaced by urban land in Beijing [22]. In the above studies, the scholars have analyzed the expansion of urban land from a regional perspective. They paid attention to the overall expansion of urban land and did not subdivide the types of urban land. Therefore, our perspective is different from theirs as we are concerned about different types of urban land in the city.

Regarding the study of urban land within a city, scholars have focused on the identification of functional urban areas or land. Big data are often used to identify functional urban space and to analyze the organization of functional urban spaces [23–26]. Tian et al. (2010) used POI data to divide the functional urban space of Beijing and pointed out that Beijing belongs to an urban spatial structure similar to that of the concentric zone mode [24]. Chen et al. (2020) compared the similarities and differences in the spatial organization of the urban functions in 25 Chinese cities with help of POI data [25]. In contrast to Tian et al.'s conclusions, Liu et al. (2021) used carpooling big data to identify the regional centers and highlighted the polycentric spatial structure of Beijing [26]. In addition, the redevelopment of a certain type of urban land occasionally appears in the existing literature. Under the influence of industrial decentralization, industrial land is constantly being replaced by other urban construction lands [27–31]. Charney (2015) analyzed the conversion between office land and residential land in the central urban area of Tel Aviv [32]. There are many other types of land use in cities, but these have not received enough attention.

The opening to the world in 1978 was the starting point for the rapid development of Chinese urbanization. After 40 years of development, the urbanization rate increased to 59.78% in 2018, which is an increase of 41.68% from 1978 [33]. This rapid urbanization has caused a fast expansion of urban space in most Chinese cities, especially in provincial capitals and cities in east China. From 2003 to 2015, urban development land increased by 22,612.2 km2 [34]. However, the rapid expansion of urban space has not automatically led to high-quality urban construction or stable use of urban land. The conversion of urban land use continued to occur in urban built-up areas [35,36]. According to William Alonso's bid-rent theory, the price of commercial land is the highest, the price of industrial land is the lowest, and the price of residential land is between the two [35,37]. Land with a lower price will be gradually replaced by land with a higher price. The excessive aspirations of local governments for rapid urban development inevitably led to a lack of vision in urban planning, which intensified the land use conversion in urban built-up areas [35]. Due to the urban land conversion and expansion, the transformation of the urban spatial structure is also happening at the same time. In terms of urban diffusion, Sargolini (2015) pointed out that three different possible scenarios can be profiled [38]. In this work, we consider two scenarios of urban land conversion and expansion. The first scenario is the continuous expansion of urban construction land to replace non-urban construction land. The second scenario is the conversion of one kind of functional urban land to another kind of functional urban land, such as the conversion of industrial land to residential land. This phenomenon is usually caused by the unreasonable arrangement of different kinds of urban construction land in the urban built-up areas due to blind and rapid urban development. In Chinese cities, especially in developing cities, these two scenarios are very common, while the second scenario occurs less frequently in Western cities [35,36].

In China, there is a large gap in the level of development between cities. Cities, such as Beijing and Shanghai, have a high level of development. The intra-urban land use has hardly changed, and the phenomenon of functional urban land replacement is not significant [35]. Land use in developing cities that are undergoing rapid urbanization is changing [10]. The expansion and replacement of urban land coexist. Changchun is a developing city experiencing rapid urbanization and was chosen as the study area here. Therefore, our research methods and conclusions are generally applicable to developing cities, especially in Asia and Africa. Cities in Asia and Africa are still undergoing land use transformation due to rapid urbanization. The rapid urbanization of Asian and African cities will continue, and the urban population will further increase. Therefore, it is necessary to explore the transformation of functional urban land use in the context of rapid urbanization in order to develop a scientific and reasonable urban spatial development strategy.

Based on this, this article will discuss the spatiotemporal evolution of the conversion and expansion of the different types of functional urban land in the central city of Changchun. Functional urban land expansion means transforming non-construction land into urban use at the edge of built-up areas. Functional urban land conversion refers to intra-urban replacement in land use and includes the following two modes: "convert in" and "convert out". The mode of "convert in" refers to the conversion of other kinds of functional urban land into a certain type. The mode of "convert out" refers to the conversion of a certain kind of functional urban land into other types. For example, the "convert in" of residential land means the conversion of other kinds of functional urban land into residential land; the "convert out" of residential land means the conversion of residential land into other kinds of functional urban land. In other words, the conversion of functional urban land is the sum of "convert in" and "convert out" modes. Once we have understood the concept of functional urban land expansion and conversion, then can the analytical method of combining functional urban land expansion and conversion trace the process of urban spatial evolution? Furthermore, is the spatiotemporal process of expansion and conversion consistent across the different types of functional urban land? As we all know, with the expansion of urban space, the importance of the sustainable development of urban space is self-evident [8]. Our research will help to explore the trajectory of urban space evolution, promote the coordinated development of different urban functions, and thus contribute to the sustainable development of urban space.

The paper is organized into the following sections: In Section 1, the research background and the purpose of the research are discussed. In Section 2, we first introduce the study area, then explain the data source, and finally elaborate the research method in detail. In Section 3, we analyze the conversion and expansion of the four types of functional urban land. In Section 4, we discuss the static and dynamic changes of functional urban land conversion and expansion from the perspective of concentric circles and sector modes. The last section is the conclusion of this study.

#### **2. Materials and Methods**

#### *2.1. Study Area*

Changchun is a developing city, located in northeast China (Figure 1). The functional urban land is undergoing drastic conversion and expansion, which supports this research. Changchun is the capital of Jilin Province. In 2015, the population of Changchun reached 7.54 million, with a total area of 20,594 km2. The central city, the study area of this article, is the core component of the city and the center of regional economic development, covering an area of 612 km2. The urban construction land of Changchun is mainly distributed in the central city, with an area of 348.6 km2 as of 2015.

Roads are the "skeleton" of the urban space, especially the ring road, which helps to analyze the concentric circle pattern of the evolution of functional urban land. The four ring roads, namely the first, second, third, and fourth ring road, divide the central city of Changchun into five zones (Figure 2). The areas of the five zones (Z-1, Z-2, Z-3, Z-4, and Z-5) are 19 km2, 52 km2, 99 km2, 128 km2, and 314 km2, respectively. Based on the five sub-zones, we can analyze the concentric features of the conversion and expansion of the functional urban land. In addition to the concentric features, this article also focuses on the sector features of the replacement and expansion of functional urban land, similar to Hoyt's sector mode. With the People's Square as the center, the central city can be divided into eight sectors according to different directions (E, W, S, N, NE, SE, SW, and NW).

**Figure 1.** Location of study area.

#### *2.2. Data Sources and Processing*

#### 2.2.1. Data Sources

The data used in this article are urban land use maps from 2003, 2007, 2011, and 2015. The source of the urban land use maps is the Changchun city master plan, which can be accessed from the website of Changchun municipal planning and natural resources bureau (http://gzj.changchun.gov.cn, accessed on 25 December 2019). Urban land use maps were the basic data of the Changchun city master plan, and urban planners spent a lot of time creating them based on topographic maps with a scale of 1:100,000 and Google Maps with a spatial resolution of 10 m. Moreover, the urban planners conducted extensive field surveys to identify the dominant uses of parcels that carry multiple functions. Therefore, urban land use maps contain very accurate classification and layout information of functional urban lands. In the maps of urban land use, urban construction land consists of the following eight types of functional urban land: residential land (RL), public service land (PL), commercial land (CL), industrial land (IL), logistic and warehouse land (LWL), road and transport facility land (RTL), green space and square land (GSL), and municipal utility land (MUL). The rest of the Changchun city master plan is also an important support for this paper. In addition, land transaction data obtained from the website of the Changchun planning and natural resources bureau was used to assist in defining the projection of the urban land use maps.

**Figure 2.** Concentric circle division of the central city.

#### 2.2.2. Data Processing

The file format of the urban land use maps we obtained as JPG. Therefore, we needed to vectorize the urban land use maps. First, the urban land use maps in JPG format were imported into ESRI ArcGIS 9.3 to define the projection (ArcToolbox, data management tools, projections and transformations, define projection). Here, the land transaction data obtained from the website of the Changchun planning and natural resources bureau was used to help georeferenced correction. The land transaction data contains detailed location information and can be easily located in the urban land use maps. Moreover, the land transaction data also contains a coordinate system. Combined with the location and coordinates of the land transaction data, the coordinates of the urban land use maps can be accurately corrected. According to the coordinate system of land transaction data, we used the WGS-84 coordinate system to define the projection. We used the coordinate points of ten transaction plots to conduct spatial corrections of the urban land use maps in JPG format. Then, we vectorized the urban land use maps. In this work, the minimum detected unit is the plot enclosed by the urban road network. After that, we used vectorized urban land use maps and the following methods to analyze the expansion and conversion of functional urban lands. In this paper, the top four types of functional urban land (i.e., RL, CL, PL, and IL) are mainly discussed. The other four types of functional urban land (i.e., LWL, RTL, GSL, and MUL) are collectively referred to as "other land" (OL). The urban land use in the central city for the four years (2003, 2007, 2011, and 2015) is shown in Figure 3.

**Figure 3.** Urban land use maps of the central city (2003, 2007, 2011, and 2015).

*2.3. Methods*

#### 2.3.1. Overlay Analysis

Overlay analysis was used to study the expansion of functional urban land. Using the ESRI ArcGIS 9.3 software platform, maps of urban land use of the four years (2003, 2007, 2011, and 2015) were overlapped, and expansions of RL, CL, PL, and IL were obtained through overlay analysis. The following equation was used to calculate the spatial distribution of the expansion of the four types of functional urban land and to analyze the spatial features of urban spatial expansion:

$$P\_{\bar{\imath}\bar{\jmath}} = \frac{AT\_{\bar{\imath}\bar{\jmath}} - At\_{\bar{\imath}\bar{\jmath}}}{AT\_{\bar{\jmath}} - At\_{\bar{\jmath}}} \tag{1}$$

where *Pij* represents the proportion of the growth area of *j*-type urban land in zone *i* to the total growth area of *j*-type urban land in the central city; *ATij* represents the area of *j*-type urban land in zone *i* in the year *T*; *Atij* represents the area of *j*-type urban land in zone *i* in the year *t*; *ATj* represents the area of *j*-type urban land in the central city in the year *T*; *Atj* represents the area of *j*-type urban land in the central city in the year *t*; and *t* represents the first temporal threshold, while *T* represents the next temporal threshold.

#### 2.3.2. Transition Matrix

The transition matrix method is an important method for the analysis of functional urban land conversion as it can be used to analyze the number and spatial distribution of conversions between different types of functional urban land [39]. The mathematical expression of the transition matrix is provided in the following equation:

$$S\_{ij} = \begin{bmatrix} S\_{11} & \cdots & S\_{1n} \\ \vdots & \ddots & \vdots \\ S\_{n1} & \cdots & S\_{nn} \end{bmatrix} \tag{2}$$

where *Sij* represents an area converted from urban land type *i* to urban land type *j*; and n represents the number of functional urban land types. ESRI ArcGIS 9.3 and Excel 2016 were used to obtain the transition matrix of functional urban land conversion. First, we used overlay analysis to obtain the intersection layer data (ArcMap, ArcToolbox, analysis tools, overlay, intersect), added a new field to the attribute table of the intersection layer data, calculated the area, and exported the attribute table to a file in DBF format. Then, we used the pivot table in Excel to obtain the transition matrix of urban land use. Thus, we can calculate "convert in" and "convert out" values of the central city and the five zones. The value of functional urban land conversion for the five zones was calculated using the following equation:

$$\mathcal{L}LR = \mathcal{C}\_{\text{in}} + \mathcal{C}\_{\text{out}} \tag{3}$$

where *ULR* is the area of functional urban land replacement; *Cin* represents the conversion of other kinds of functional urban land into this kind of functional urban land; and *Cout* is the opposite of *Cin*.

#### 2.3.3. Kernel Density Estimation

Through the previous overlay analysis and transfer matrix, we have obtained the spatial layout map of functional urban land expansion and conversion. We used the kernel density estimation method for the analysis of the spatial distribution density of functional urban land with the aim of discovering where the land use changes are spatially clustered. The kernel density was calculated as follows:

$$f(\mathbf{x}) = \frac{1}{nh^d} \sum\_{i=1}^{n} K(\frac{\mathbf{x} - \mathbf{x}\_i}{h}) \tag{4}$$

where *K*( *<sup>x</sup>*−*xi <sup>h</sup>* ) represents the nuclear density equation; *h* is the threshold; *n* is the number of points within the threshold; and *d* is the dimension of the data.

#### **3. Results**

#### *3.1. General Characteristics of Functional Urban Land Conversion and Expansion*

Judging by the composition of the functional urban land in the central city of Changchun, the sum of the proportions of RL, IL, PL, and CL exceeded almost two-thirds. From 2003 to 2015, the area of urban construction land in the central city increased by 111.24 km2, of which the four types of functional urban land increased by 81.39 km2, and their contribution rate reached 73.17%. RL and IL have always been the two largest land types in the central city. The proportion of RL exceeded 30%, while the proportion of IL remained at around 22%.

The area of RL, CL, PL and IL in the three stages increased by 31.13 km2, 8.66 km2, 6.88 km<sup>2</sup> and 44.48 km2, respectively, through functional urban land expansion. The expansion of the functional urban land gradually shifted from the inner ring to the outer ring (Figure 4), and mostly occurred outside of the second ring road from 2003 to 2007, and outside of the third ring road in the latter two stages. From 2003 to 2007, the expansion of the four kinds of functional urban land outside of the second ring road accounted for 98%, including 26% in Z-3, 20% in Z-4, and 52% in Z-5. From 2007 to 2011, the expansion of the four kinds of functional urban land outside of the third ring road accounted for 86%, including 27% in Z-4 and 59% in Z-5. From 2011 to 2015, the expansion of the four kinds of functional urban land outside of the third ring road accounted for 93%, including 25% in Z-4 and 68% in Z-5. The functional urban land was mainly expanded to the southwest and southeast in the period from 2003 to 2015. Furthermore, from 2003 to 2007, the four kinds of functional urban land also expanded to the northwest and north. In the latter two stages, in addition to the southwest and southeast, the northeast and east were also the main directions for the functional urban land expansion. According to the Changchun city master plan, the government's urban development strategy guided the direction of the urban spatial expansion.

**Figure 4.** Overall characteristics of the kernel density of the functional urban land expansion. This figure is the kernel density map of the total expanding area of the four types of functional urban land, so the kernel density value is greater than that of a certain type of functional urban land.

Similar to the expansion of functional urban land, the area of functional urban land conversion was also very large. RL, CL, and PL increased by replacing other types of functional urban land, while IL decreased (Table 1). RL increased by 1105.17 ha by "convert in" (2109.73 ha) and "convert out" (1004.56 ha) modes, while CL increased by 100.88 ha and PL increased by 220.08 ha. The "convert in" mode of industrial land was much lower than the "convert out" mode, which resulted in a reduction in IL by 2402.06 ha. The conversion of the functional urban land mainly occurred in Z-1 during the first period, and mainly occurred in Z-3 from 2007 to 2015, which showed a trend of expansion from the inner ring to the outer ring. However, the third ring road was the spatial boundary of the conversion of functional urban land (Figure 5).

**Figure 5.** Overall characteristics of the kernel density of functional urban land conversion. This figure is the kernel density map of the total conversion area of the four types of functional urban land, including the area of "convert in" and "convert out", so the kernel density value is greater than that of a certain type of functional urban land.


**Table 1.** Transition matrix of functional urban land replacement from 2003 to 2015. Unit: ha.

#### *3.2. Spatiotemporal Analysis of Conversion and Expansion of Residential Land*

With further development, RL increased by 1357.58 ha between 2003 and 2007, 1851.99 ha between 2007 and 2011, and 1478.21 ha between 2011 and 2015. The expansion of RL mainly occurred outside of the second ring road and continued towards the outer ring road. The proportions of RL expansion beyond the second ring road in the three phases were 98.86%, 95.27%, and 99.72%, respectively; the proportions of RL expansion beyond the third ring road were 53.50%, 73.97%, and 91.64%, respectively; the proportions of RL expansion beyond the fourth ring road were 36.50%, 33.36%, and 62.49%, respectively. The residential space was constantly expanding from the inner circle to the outer circle (Figure 6). The directions of residential land expansion were heterogenous, changing from southeast, northwest, and north in 2003–2007 to southeast and east in 2007–2011, and towards the south, north, and northeast in 2011–2015.

**Figure 6.** Kernel density map of residential land expansion in different periods.

The conversion of RL was intense from 2003 to 2015, and the "convert in" area of RL was significantly larger than the "convert out" area. Through replacement, residential land increased by 127.16 ha between 2003 and 2007, 569.11 ha between 2007 and 2011, and 408.90 ha between 2011 and 2015. The "convert in" mode of RL during the three periods mainly occurred in Z-1, Z-2, and Z-3. (Figure 7). Due to the suburbanization of industry, IL has mainly been replaced by RL in the inner city. In the east and northeast of the central city, a large area of IL has been converted into RL. The "convert out" of RL mainly occurred in Z-1, Z-2, and Z-3 from 2003 to 2011, and in Z-1, Z-2, Z-3, and Z-4 from 2011 to 2015. Overall, the central city's RL increased continuously from 2003 to 2015.

#### *3.3. Spatiotemporal Analysis of Conversion and Expansion of Commercial Land*

Judging by the average area of the expanding CL plots, the increase in commercial space was greatly influenced by the construction of large-scale flagship projects. Often, the development of flagship projects leads to concentric and directional characteristics of commercial land expansion (Figure 8). Combined with the field investigation, we found that the construction of the automobile trade city (on an area of 40 ha) led to the expansion of CL to the west during 2003–2007. The construction of a number of automobile sale and service shops (4S) led to the expansion of CL to the southwest during 2007–2011. The

construction of JingYue CBD led the expansion of CL to the southeast. Since these large commercial facilities cover a large area, they were built on the edge of the central city. Therefore, the commercial space of the central city gradually expanded towards the outer circle, and the proportion of CL expansion outside of the fourth ring road was 34.02%, 58.88%, and 92.00% in the three stages, respectively.

**Figure 7.** Kernel density map of residential land conversion in different periods.

**Figure 8.** Kernel density map of commercial land expansion in different periods.

In general, with the replacement of functional urban land, the area of CL continued to expand, increasing by 328.07 ha, 85.91 ha, and 38.42 ha in the three stages, respectively. The "convert in" mode of CL mainly occurred within the first ring road during 2003–2007 and was mostly within the third ring road during 2007–2015. Compared to the first period, the "convert in" mode of CL was expanded in the latter two stages. CL has largely replaced IL, indicating that the process of "suppress the second industry and develop the third industry" is continuing. The price of CL is generally higher than that of the other types of functional urban land. Therefore, it is generally rare to replace CL with other kinds of functional urban land. However, in this study, a large area of CL has been replaced by other kinds of functional urban land, especially RL and PL. Field research has shown that CL, replaced by other urban land types, was usually low-grade commercial facilities with chaotic layout and wasted land. The conversion of the functional urban land has thus promoted the integration of this inefficient land into the surrounding built-up environment. The "convert out" mode of CL mainly occurred in Z-1 and Z-2 from 2003 to 2007 and in Z-3 and Z-4 from 2007 to 2011 (Figure 9). In the third stage, the area of the "convert out" of CL decreased sharply and was concentrated only towards the northeast in Z-3.

**Figure 9.** Kernel density map of commercial land conversion in different periods.

#### *3.4. Spatiotemporal Analysis of Conversion and Expansion of Public Service Land*

The expansion of PL over the three phases showed that the public service space mainly expanded towards the outside of the southeast in Z-5 (Figure 10). According to the city master plan of Changchun, the ecological environment in the southeast of Changchun is beautiful, and the forest is widespread. The master plan pointed out that the Changchun municipal government has been committed to transforming the area into a "University Town" and a "New Livable Town" since the mid-1990s. Accordingly, universities and exhibition facilities gradually gathered in the area. The relocation of Jilin Jianzhu University and the expansion of Jilin Agricultural University led to this expansion of PL during 2003–2007. In 2007–2011 and 2011–2015, the construction of convention and exhibition facilities, as well as the expansion of universities led to an expansion of PL. The Changchun modern agricultural park, which covers an area of 106 ha, was formally established in 2009. A scientific research institution, the Chinese Academy of Agricultural Sciences, was relocated to the "University Town" in 2011. The educational and exhibition facilities have led to an expansion of PL.

The conversion of PL became less and less frequent from 2003 to 2015. During the period of 2003–2007, an area of 358.97 ha was converted from PL to other urban land types. The area that was converted from the other urban land types to PL was 563.85 ha, resulting in an increase of 204.88 ha in PL. The area of "converted in" was 168.52 ha, and the area of "converted out" was 153.84 ha, indicating an increase of 14.98 ha from 2007 to 2011. From 2011 to 2015, the area of "converted in" was 67.72 ha, and the area of "converted out" was 57.65 ha, indicating an increase of 10.07 ha. The "convert in" and "convert out" modes of PL mainly occurred in Z-1, Z-2, and Z-3 (Figure 11). Similar to the conversion of CL, the replacement of PL has shown that the urban service functions needed to be further improved, especially in Z-1, Z-2, and Z-3. With the continuous improvement of the urban

service functions in the inner city, the conversion of PL can be moved from the inner circle to the outer circle. The conversion between residential, commercial and public service lands occurred frequently. From 2003 to 2007, 53.56% of the PL, converted from other urban land types, consisted of RL and CL. However, 75.84% of the PL converted to other urban land types consisted of RL and CL. From 2007 to 2011, 35.25% of the PL, converted from other urban land types, consisted of RL and CL. However, 78.88% of the PL, converted into other urban land types, consisted of RL and CL. From 2011 to 2015, 50.25% of the PL, converted from other urban land types, consisted of RL and CL. However, 99.92% of the PL, converted to other urban land types, consisted of RL and CL.

**Figure 10.** Kernel density map of public service land expansion in different periods.

**Figure 11.** Kernel density map of public service land conversion in different periods.

#### *3.5. Spatiotemporal Analysis of Conversion and Expansion of Industrial Land*

The expansion of IL showed predictable spatial characteristics. From 2003 to 2015, there was an expansion of IL mainly towards the southwest in Z-5. The northeast expansion beyond the fourth ring road was added during the latter two stages (Figure 12). It is clear that the southwest and northeast were the main directions of the expansion of IL, while Z-5 was the main circle of IL expansion. Between 2003 and 2007, the proportion of IL expansion in Z-5 accounted for 66.40%. IL expansion in the southwest direction was 44.91%. During

the latter two stages, the proportions of IL expansion in Z-5 were 84.19% and 79.95%, respectively. From 2007 to 2011, the proportions of IL expansion in the southwest and northeast were 41.57% and 34.63%, respectively. From 2011 to 2015, the proportions of the expansion of IL in these two directions accounted for 38.68% and 36.80%, respectively. According to the city master plan of Changchun, the automobile industry is the pillar industry of Changchun, and the expansion of the automobile industry in the southwest is the main reason for the observed IL expansion.

**Figure 12.** Kernel density map of industrial land expansion in different periods.

From 2003 to 2015, the conversion of IL was very frequent, and IL was continuously replaced by other urban land types (Figure 13). It is clear that IL is constantly being replaced by other urban land types. During all of the three periods, the areas of IL that replaced other urban land types were 834.66 ha, 94.63 ha, and 1.38 ha, respectively; the areas of IL converted into other urban land types were 818.16 ha, 667.46 ha, and 348.03 ha, respectively. The IL "converted out" area exceeded the IL "converted in" area. By replacing IL, it increased by 16.5 ha between 2003 and 2007, decreased by 572.83 ha between 2007 and 2011, and decreased by 346.64 ha between 2011 and 2015. IL was constantly replaced by RL. In all of the three stages, the proportions of the area converted from IL to RL were 38.50%, 75.03%, and 85.52%, respectively. The "convert in" mode of IL mainly occurred in Z-2, Z-3, and Z-4 between 2003 and 2007. In the latter two stages, the "convert in" mode of IL was not obvious. The "convert out" mode of IL in the three periods mainly occurred in Z-1, Z-2, and Z-3. IL in the inner city was gradually moved to the periphery of the city (industrial suburbanization).

**Figure 13.** Kernel density map of industrial land conversion in different periods.

#### **4. Discussion**

The existing literature on urban space mainly uses urban land expansion data to study urban spatial evolution by comparing the structure of urban space in different years [16–22]. In this study, we explored the dynamic changes of urban space from two perspectives of functional urban land expansion and conversion. Our research demonstrated that functional urban land conversion, similar to functional urban land expansion, is also very common in developing cities. Urban land use maps have facilitated this work. In fact, scholars have been able to identify the different kinds of functional urban spaces or areas through a variety of data [11–13,26]. We suggest that scholars can use the above identified functional urban space or area data to further analyze the replacement between different functional urban land types in order to explore the inherent characteristics of urban spatial evolution. Our work showed that the combination of urban land expansion and conversion can effectively reveal changes in intra-urban space. The structure of urban land use changed slightly from 2003 to 2015. The proportions of RL and CL increased from 27.82% and 3.71% in 2003 to 31.04% and 5.30% in 2015, respectively. The proportions of PL and IL decreased from 11.48% and 23.02% in 2003 to 10.42% and 21.54% in 2015, respectively. The change in urban land use structure is the joint effect of expansion and conversion, but we often ignore the conversion between different types of functional urban land. In addition, using urban land use data at four-year intervals (2003, 2007, 2011, and 2015), our findings prove that functional urban land expansion and conversion in developing cities are very frequent under the process of rapid urbanization. Usually, the topic of urban spatial evolution is studied over a long time period, such as five or ten years [17,18,20,22]. However, our research suggests that for developing cities, the time interval can be shortened, such as the four-year period in this article.

On the basis of the research conducted by Tian et al. (2010) [24], we not only analyzed the concentric characteristics of the functional urban space, but also analyzed the concentration along spatial directions. The concentric characteristics of the expansion of different types of functional urban land are similar, while the directional characteristics are different. RL, CL, PL, and PL have all expanded into the outer circle. RL was mostly expanded into Z-3, Z-4, and Z-5, while CL, PL and IL were mainly expanded into Z-5. RL gradually

expanded from the inner circle to the outer circle; the expansion of IL to the periphery was the result of industrial suburbanization, and the construction of large flagship projects has directed the expansion of CL and PL to the periphery. Influenced by urban development strategies and urban planning, the directions of expansion of the four types of functional urban land have shown different characteristics. Southeast, southwest, and south were the main directions for expanding RL; CL expanded in many directions, but there was only one main direction (southeast) for the expansion of PL, and IL has largely expanded to the southwest and northeast.

The conversion of the functional urban land in Z-1, Z-2, and Z-3 showed the characteristics of the transition from the inner circle to the outer circle. The conversion of functional urban land beyond the third ring road was rare. Only the "convert out" mode of RL and CL and the "convert in" mode of IL in a certain period could happen in Z-4. RL, CL, and PL were constantly replacing other types of functional urban land, while IL was constantly being transformed into other types of functional urban land. The conversion area of the four types of functional urban land in Z-1 was large in the first phase and small in the latter two phases. The Z-2 and Z-3 were zones where the area of functional urban land conversion was large in the latter two phases. The "convert in" mode of RL, CL, and PL was mainly located in Z-1, Z-2, and Z-3, which was a process of continuous improvement of the urban functions of the inner city. In the first period, the "convert in" mode of IL was larger than the other three types of functional urban land, and mostly concentrated in Z-2, Z-3, and Z-4. However, in the latter two periods, the "convert in" mode of IL was rare. The "convert out" mode of RL and CL changed from being concentrated in Z-1, Z-2, and Z-3 to being concentrated in Z-1, Z-2, Z-3, and Z-4. The redistribution of public service facilities and the suburbanization of industries have led to the "convert out" mode of PL and IL mainly in Z-1, Z-2, and Z-3.

Our research found that the conversion and expansion of the functional urban land in the central city of Changchun show similar characteristics as the concentric circle and the sector modes. The conversion and expansion of the functional urban land show the obvious characteristics of concentric circles. Influenced by urban planning and urban development policies, the conversion and expansion of functional urban land show predictable fanshaped (directional) characteristics. The Burgess concentric zone theory and Hoyt's sector theory are static modes of urban spatial structure [24,40,41]. However, the concentric zone mode and the sector mode are not static for a developing city but are dynamic. In other words, the static layout of urban space exhibits concentric circles and fan-shaped characteristics, and the evolutionary process of urban space also exhibits concentric circles and fan-shaped characteristics for developing cities. The conversion and expansion of the different functional urban lands show different concentric circles and fan-shaped characteristics at different periods. The expansion of RL, CL, PL, and IL shifts from the second ring road to the third ring road. The direction of expansion of RL, CL, PL, and IL in the three stages is different. The conversion of RL, CL, PL, and IL also show dynamic changes in the concentric circles and in the fan-shaped characteristics.

Our research will provide theoretical support for the coordinated development of different types of functional urban space in different areas of a city. The coordinated development of urban space will further contribute to the optimization of urban spatial structure and will promote the sustainable development of urban space. However, this study has one main shortcoming. We only studied urban space from the perspective of functional urban land. Urban space is complex, and apart from land, it also covers other elements such as population, architecture, transportation, society, and culture. In general, it is biased to explore the evolution of urban space using a single criterion rather than multiple criteria. We therefore should consider additional development elements to overcome this shortcoming in future research, rather than only using a single criterion. Moreover, it should also be noted that urban land use is one of the most important indicators of urban space and can effectively decode the geospatial pattern of cities. Our study proves that the phenomenon of intra-urban land transformation brought by rapid urbanization is common and needs to draw the attention of urban policy makers and decision makers to scientifically adjust the strategy of urban spatial development. We hope that this research will be open to other researchers who study urban space from different perspectives in order to obtain inspiration for urban research. In addition, the quantitative methods we used are retrospective tools. We will introduce predictive tools to study the temporal and spatial changes in the replacement and expansion of functional urban spaces in the future in order to better optimize the urban spatial structure.

#### **5. Conclusions**

With the help of urban land use maps and other auxiliary data, we analyzed the temporal and spatial characteristics of the replacement and expansion of different types of functional urban land using overlay analysis, transition matrix, and kernel density estimation methods. The results prove that the analysis of the replacement and expansion of functional urban land can reflect the dynamic evolution of urban space in the context of rapid urbanization. This research showed that changing the structure of land use in cities and expanding and replacing the functional urban land can be used to analyze whether the urban spatial structure has stabilized. The conversion and expansion of functional urban land in the central city of Changchun is still very active, which indicates that urban spaces have not yet stabilized and that they still require continuous improvement. When one kind of functional urban land replaces other kinds of functional urban land, there is further replacement by other kinds of functional urban land. This study demonstrated changes between RL, CL, PL, and IL. Combined with changes in the structure of the land use in cities, suggestions can be made for the future layout of urban land, and the urban spatial structure can be adjusted in order to promote sustainable urban spatial development.

The temporal and spatial characteristics of urban spatial evolution are a classic research theme involving multiple disciplines of urban planning, geography, and economics. The spatiotemporal changes in the functional urban land evolution and the path of optimization of urban spatial structure were explored in order to support the reasonable development of urban space. As we mentioned in the Introduction section, our research methods and conclusions are generally applicable to developing cities around the world in order to help urban policy makers propose sound strategies for the scientific development of urban spaces.

**Author Contributions:** Conceptualization, G.Z. and C.L.; methodology, G.Z. and Y.L.; software, G.Z.; formal analysis, J.Z.; investigation, G.Z. and J.Z.; writing—original draft preparation, G.Z. and J.Z.; writing—review and editing, G.Z., C.L. and Y.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China, grant number 42171191, 41771172, 41871158; the China Postdoctoral Science Foundation, grant number 2018M641760; the Education Department of Jilin Province, grant number JJKH20201173KJ.

**Acknowledgments:** The authors acknowledge the journal editors and guest editors of the Special Issue "Land Use Change from Non-urban to Urban Areas: Problems, Challenges and Opportunities" for their support and the granting of funding covering the fees related to the publishing of this paper.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Spatiotemporal Dynamics of Soil Impermeability and Its Impact on the Hydrology of An Urban Basin**

**Fernando Oñate-Valdivieso \*, Arianna Oñate-Paladines and Milton Collaguazo**

Departamento de Ingeniería Civil, Universidad Técnica Particular de Loja, C/. Marcelino Champagnat S/N, Loja 1101608, Ecuador; ayonate@outlook.es (A.O.-P.); mecollaguazo1@utpl.edu.ec (M.C.)

**\*** Correspondence: fronate@utpl.edu.ec

**Abstract:** The presence of impervious surfaces in catchments interferes with the natural process of infiltration, which has a marked influence on the hydrological cycle, affecting the base flow in rivers and increasing the surface runoff and the magnitude of flood flows. Like many Latin American cities, Loja (located in southern Ecuador) has experienced significant rates of urban growth in recent years, increasing the impervious surfaces in the catchment where it belongs. The aim of this study is to analyze the spatiotemporal dynamics of imperviousness in the study area for the period 1989–2020, using the Normalized Difference Impervious Surface Index (NDISI) and the supervised classification of Landsat images. The effect on flood flows was studied for each timestep using HEC-HMS hydrological model. Additionally, a future scenario of impervious surfaces was generated considering the observed spatiotemporal variability, possible explanatory variables, and logistic regression models. Between 1989 and 2020, there was an increase of 144.12% in impervious surfaces, which corresponds to the population growth of 282.56% that occurred in the same period. The period between 2001 and 2013 was the one that presented the most significant increase (1.06 km2/year). A direct relationship between the increase in impervious surfaces and the increase in flood flows was observed, reaching a significant variation towards the horizon year that could affect the population, for which measures to manage the surface runoff is necessary.

**Keywords:** urban hydrology; impervious surfaces; land use scenarios; urban surface growth; hydrological model; flood flows

#### **1. Introduction**

Features such as climate, topography, vegetation, and coverage of a natural watershed produce a natural water cycle and a given hydrological response. Different factors such as the impervious surfaces can affect this unique natural hydrological process and cause adverse effects to the catchment [1].

The impervious surface is usually defined as the collection of anthropogenic landforms that water cannot directly infiltrate into, including rooftops, roads, and parking lots [2,3]. The urbanization process has significant impacts on the hydrology of a basin; as urban areas expand, permeable and moisture-holding lands transform into impervious surfaces such as concrete and asphalt, causing a decrease in infiltration and base flow, as well as an increase in flood flows and runoff volumes. The storm drainage systems simplify the natural drainage systems, altering the response of the basin to precipitation events since shorter concentration and recession times occur [1,4,5]. On the other hand, the dynamics of impervious surfaces impact urban regional climate by altering the thermal environment and water quality [3,6].

Several studies have analyzed the effect of urbanization [7–9], land-use changes [10–14], or impervious cover change [15,16] on Hydrology.

Remote sensing has been extensively utilized for the detection of impervious surfaces [17,18]. The approaches have been diverse: index-based methods, classification-based

**Citation:** Oñate-Valdivieso, F.; Oñate-Paladines, A.; Collaguazo, M. Spatiotemporal Dynamics of Soil Impermeability and Its Impact on the Hydrology of An Urban Basin. *Land* **2022**, *11*, 250. https://doi.org/ 10.3390/land11020250

Academic Editors: Víctor Hugo González-Jaramillo and Antonio Novelli

Received: 13 December 2021 Accepted: 28 January 2022 Published: 8 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

methods, and mixture analysis. The index-based methods include indices such as normalized difference buildup index (NDBI) [19], normalized difference impervious surface index (NDISI) [20], modified NDISI (MNDISI) [21], biophysical composition index (BCI) [22], and perpendicular impervious surface index (PISI) [23]. The classification and regression approaches include maximum likelihood classifier [24], support vector machine [25], artificial neural networks [26], random forest [27], and object-oriented methods [28]. For the mixture analysis, spectral mixture analysis (SMA) [29] and temporal mixture analysis (TMA) [30] are applied.

Urbanization is a worldwide trend. Currently, more than 50% of the world's population lives in urban centers, and more than 500 cities in the world have a population above 1 million inhabitants [5]. The reasons for urban growth are diverse; in Latin American cities, we could highlight the natural demographic growth, migration from the countryside to the city in search of better living conditions, changes in the location patterns of economic activities, and housing, among others.

Several cities in Ecuador have experienced rapid growth, which is evidenced in a notable increase in the urbanized area in recent years. One of those cities is Loja, capital of the province of the same name, located south of Ecuador and bordering Peru. This study analyzes the influence of urban growth on the hydrology of the basin where the city is located and on the extreme flow events that occur in it. For this, using aerial photographs and satellite images, a Normalized Difference Impervious Surface Index (NDISI) was calculated, and a multitemporal analysis of the urban surface variation was carried out. Then, flood flows for various coverage scenarios were generated using precipitation data and applying a hydrological model to finally evaluate the effect of these flows over the areas surrounding the riverbanks in various points of interest. The study of the spatiotemporal variation of the impervious surfaces using a spectral index and supervised classification of images, combined with statistical techniques and artificial intelligence to define future scenarios of impervious surfaces, and the evaluation of their possible impacts through hydrological modeling, are the newest aspects of the present work.

#### **2. Materials and Methods**

#### *2.1. Study Area*

The Zamora River (A = 227 km2) is a tributary of the Santiago River and part of the hydrological system of the Amazon River. The Zamora River basin is located in the southern Andes of Ecuador, has an average height of 2400 m above sea level, an average slope of 30%, and an average slope of the main channel of 8.3% [31]. The basin is covered by vegetation in good condition, mainly composed of grasslands, scrublands, and forests [32]. Its climate is subhumid equatorial temperate, with a mean annual precipitation depth of 909.1 mm. The Zamora River presents dry periods between May and November and can present important flows during the rainy season (from December to April) [33]. The Zamora River, up to its pour point (79◦13 28" W, 3◦55 17" S) has six main tributaries that make up a network of 102.70 km in length, with the main channel of 22.89 km, which has stream order three, according to the Horton—Strahler Laws. The city of Loja occupies the middle and lower portion of the basin. The city has about 200,000 inhabitants and an area of 43 km2, being the only existing urban area in the Zamora River basin. The growth of the city in the last 30 years, as well as the construction and improvement of the road network, has created impervious zones.

The location of the study area is presented in Figure 1.

#### *2.2. Data Collection*

Three image sets, acquired from Landsat 5 Thematic Mapper (TM), Landsat 7 Enhanced Thematic Thematic Mapper Plus (ETM+), and Landsat 8 Operational Land Imager (OLI)-Thermal Infrared Sensor (TIRS) were collected in the study area [34]. Their acquisition dates, spectral bands, and spatial resolutions are listed in Table 1. Atmospheric correction was performed to each image using the Atmospheric/Topographic Correction

for Mountainous Terrain (ATCOR) software developed by the German Aerospace Center, Wessling, Germany [35]. The Landsat images archived in the U.S. Geological Survey (USGS) data clearinghouse have been georectified [36]. All images in Table 1 are geometrically matched to each other.

**Figure 1.** Location of the study zone.

Further, historical information on the type and land use of the study area was collected [37], which was used in the hydrological modeling component of this study.

#### *2.3. Analysis of Impervious Surfaces*

The Normalized Difference Impervious Surface Index (*NDISI*) [20,38] is used to enhance impervious surfaces and suppress land covers such as soil, sand, and water bodies.

$$NDISI = \frac{T\_b - (MNDWI + NIR + SWIR1)/3}{T\_b + (MNDWI + NIR + SWIR1)/3} \tag{1}$$

*Tb* refers to the brightness temperature of the TIRS1 thermal band. *MNDWI* represents the Modified Normalized Difference Water Index (Equation (2)), *NIR* refers to the pixel values extracted from the near-infrared band. *SWIR1* refers to the pixel values extracted from the first shortwave infrared band.

$$MNDWI = \frac{G - SWIR1}{G + SWIR1} \tag{2}$$

*G* represents the pixel values extracted from the green band.


**Table 1.** Satellite images from the three Landsat sensors (TM, ETM+, OLI-TIRS) used in this study.

Applying Equation (1), NDISI images were generated for each of the collected images (Table 1). A manually adjusted threshold was used to extract impervious surface features from the NDISI images generated. The pixels with values greater than the threshold are impervious surfaces and were assigned a value of 1, while the pixels with values equal to or less than the threshold are nonimpervious surfaces and were assigned a value of 0. Thus, the resultant image is a binary image, only showing the extracted impervious surfaces.

Additionally, the supervised classification of the collected images was carried out using the maximum likelihood method [39] in order to obtain the urban area (impervious surface) and its temporal variation and maps of the impervious and nonimpervious surface for the study area.

The performance of the NDISI and the supervised classification for the detection of impervious surfaces was evaluated by visual comparison with the images included in Table 1. Using the results of the technique that offered the most reliable results, the spatiotemporal analysis was performed, as well as the generation of scenarios towards the year 2030 and hydrological modeling in order to study the impact of the variation of impervious surfaces on the hydrology of the basin under study.

#### *2.4. Spatiotemporal Analysis of Impervious Surfaces and Scenario Generation*

Once the maps of the impervious and nonimpervious surface were obtained, the changes that occurred between 1989 and 2001 were analyzed, relating them to the possible explanatory variables to obtain a predictive model that could be validated by comparison with the coverage obtained for 2013.

The changes that occurred were studied by applying the methodology proposed by [40], which allows determining the persistence, gain, loss, and exchanges between the thematic categories considered in each land occupation map through the analysis of a cross-tabulation, identifying the transitions that occurred between 1989 and 2001. The relationships between the observed transitions and their possible explanatory variables are called transition submodels. The number of transition submodels will be equal to the number of transitions that occur in the study area; it is possible to group several transitions into a single model when it is considered that these are the product of the same causes. Each transition model includes a certain number of explanatory variables, which can be selected based on their explanatory potential, calculated by Cramer's V coefficient, or by testing various combinations of explanatory variables until the optimal fit between transitions and explanatory variables is obtained. Cramer's V values greater than 0.4 are acceptable [41]. Three explanatory variables were considered: Elevation (using a digital elevation model—DEM), which influences the presence of different types of vegetation; the slope, which limits urban growth; and the distance to streets and roads, which motivates and facilitates urban growth.

The transition submodels were calculated by logistic regression and by means of a multilayer perceptron neural network (MLP), obtaining the probability of occurrence of each transition according to the selected explanatory variables. Logistic regression [42] allows establishing a relationship between a binary dependent variable (transitions) and the explanatory variables considered, modeling their probability of occurrence according to the latter.

Neural networks of multilayer perceptrons are formed by a set of simple elements (neurons or perceptrons) distributed in layers and are connected to the intermediate layer or layers by means of activation functions. These functions are defined from a series of weights or weighting factors that are calculated interactively in the learning process of the network. The objective of this learning is to estimate known results (observed transitions) from some input data (explanatory variables); to later calculate unknown results from the rest of the input data. Learning is carried out from all the units that make up the network, varying the set of weights in successive interactions [39].

The land cover change modeling towards the horizon year (2013) was carried out applying Markov chains; using the land cover map of the end date (2001) along with the transition probability matrix previously calculated, to determine the zones are that will undergo a transition from the end date to the prediction date (2013).

The future land cover map was modeled using a multiobjective land-use allocation procedure (MOLA) [41,43]. Considering all transitions and using the selected explanatory variables, a list of host classes (which would lose some area) and a list of demanding classes (which would gain some area) are created. Loss or gain areas are determined by Markov chains and through the multiobjective allocation procedure, in which the explanatory variables determine the most suitable places for each change in occupation. Land from all host classes is allocated to all demanding classes. The results of each land occupation reallocation are overlaid to produce the final result [41]

Two maps were generated to predict land cover for the year 2013 based on the modeling of the relationships between the observed changes and the explanatory variables. These relationships were modeled with logistic regression and neural networks. For the validation, the map extracted from the 2013 image was considered as a reference, and, through confusion matrices, the correspondence between the reference map and those obtained through neural networks and logistic regression was studied. Forecast errors of land cover were determined for each model proposed, as well as omission and commission errors that may have occurred.

From the confusion matrix, the global reliability of the classification was calculated as the relationship between the number of pixels correctly assigned and the total number of pixels in the image [39]. Complementarily, the fit between the reference map and the maps generated was calculated using the Kappa index [40]. After analyzing the adjustment, we proceeded to generate a land cover map towards the year 2030, considering the land cover maps of 2013 and 2020, the explanatory variables selected for each transition, and applying the model that presents the best capacities.

The spatiotemporal analysis of impervious surfaces and scenario generation described was carried out by applying the land change modeler module of TerrSet 2020 [41].

#### *2.5. Hydrological Modeling*

The HEC-HMS model was developed to study the response of the Zamora River basin to extreme precipitation events, considering the different stages of urban growth in the city of Loja. The basin topology was developed based on a digital elevation model generated using a contour map at 1:50,000 scale [44]. This topological model included contributing sub-basins, junction points in which the contributions of the sub-basins are added, sections of the river network in which the hydrologic routing of the hydrographs is carried out, and the outlet point of the basin in which the flow resulting from the rain-runoff simulation is obtained. The topological model is presented in Figure 2.

**Figure 2.** Topological model of the Zamora River basin.

Synthetic storms were generated for return periods of 10, 25, 50, and 100 years, using intensity equations determined for the city of Loja [45].

$$I\_{TR} = 92.854Id\_{TR}t^{-0.4083} \tag{3}$$

$$I\_{TR} = 480.74Id\_{TR}t^{-0.8489} \tag{4}$$

where *IdTR* is the maximum intensity for a given return period, *t* is the duration of the storm in minutes, *ITR* is the intensity in mmh<sup>−</sup>1. Equation (1) is valid for durations between 5 and 43 min, Equation (2) is valid for durations between 43 min and 1440 min.

Abstractions were quantified using the curve number (*CN*) methodology of the U.S. Soil Conservation Service (USSCS) [4,46] for normal conditions, calculating the *CN* for each hydrological response unit obtained according to the intersection of type and land use for each date considered. The transformation of surface runoff into flow was carried out by applying the USSCS Unit Hydrograph. For the hydrologic flow routing, the Muskingum-Cunge method was applied. The concentration and delay times of each sub-basins were determined using the Kirpich formula [4,46].

#### **3. Results and Analysis**

#### *3.1. Analysis of Impervious Surfaces Using NDISI*

Figure 3 shows the temporal variation of the NDISI index. A visual comparison with the collected images allowed us to determine that the consolidated areas of the city center are identified in an acceptable way through the NDISI index. The areas surrounding the city center were consolidated as urban areas over time, and in the process, a transition is observed from the heterogeneous mixture of impervious and green areas to consolidated urban areas. The impervious surfaces of the southwestern portion of the city were underestimated in all analyzed images. Land surface emissivity (ε) varies with land cover on the ground surface. In urban environments, surfaces with vegetation have a higher thermal retention capacity and, therefore, have greater cooling effects than areas without vegetation [47], which is reflected in the temporal variation of the NDISI.

**Figure 3.** Spatiotemporal dynamics of impervious surfaces using NDISI: Period 1989–2020.

In suburban–rural areas, impervious areas have been identified to the west of the city. These areas do not correspond to urban areas but to surfaces with bare soil due to fallow agricultural areas or small areas under construction that have just been cleared. On the other hand, the selected images were taken between August and November (which are part of the dry season) to ensure less cloud cover. Therefore, during that period, the vegetation cover is less vigorous and frequently leaves the soil exposed. In Landsat images, bare soil is often mistaken for impervious surfaces due to their similar spectral characteristics, resulting in noisy salt and pepper appearances in supervised image classification [47]. Furthermore, the thermal response of the soil is quite similar to that of the impermeable surface, which causes spectral confusion between impermeable areas and bare soil when classifying it [20,48].

Despite its acceptable performance, it was considered that the NDISI and its temporal analysis were not completely adequate to study the evolution of the impervious surfaces in the study area.

#### *3.2. Analysis of Impervious Surfaces by Supervised Classification*

Figure 4 shows the impervious surfaces in the study area determined by supervised classification applying the maximum likelihood criterion. A visual comparison with the collected images shows an adequate representation of the impervious surfaces and their spatiotemporal variation, which is why they are selected for the following phases of the study.

**Figure 4.** Spatiotemporal dynamics of impervious surfaces using supervised classification for 1989–2020.

A summary of the area occupied by the city (impervious surfaces) and its respective population at each year considered is included in Table 2.

**Table 2.** Variation of urban area and population in the city of Loja. Period 1989–2020.


\* Reference year: 1989.

Table 2 includes a summary of the area occupied by the city (impervious surface), which was determined through the supervised classification, as well as its respective population at each year considered. For the period between 1989 and 2020, there is an increase of 144.12% in impervious surfaces, which corresponds to the population growth of 282.56% that occurred in the same period. Furthermore, there is a very significant increase in the annual variation of the impervious surfaces in the period between 2001 and 2013 (1.06 km2/year), which is linked to the receipt of money remittances sent by a large number of Ecuadorian citizens who emigrated overseas as of 1999. The population density shows a significant growth between 1989 and 2001, but in 2013 it reduced probably due to the mentioned migratory process that Ecuador experienced during the first decade of this century. By 2020, the population density recovers an increasing trend.

Figure 5 presents the variability of the impervious surfaces in the periods 1989–2001, 2001–2013, and 2013–2020. In the period 1989–2001, growth was observed based on the consolidation of the areas adjacent to the downtown, with the areas located to the southeast and north of the city achieving further development. The greatest increase in impervious surfaces occurred between 2001 and 2013, with the highest incidence in the southwest of the city, which, at that time, already had basic infrastructure which facilitated urban development. Something similar was observed in the east of the city, although on a smaller scale. For its part, in the 2013–2020 period, urban development was directed towards the west of the city, which, due to its better topographic conditions, has become the ideal place for the growth of the city.

**Figure 5.** Variability of impervious surfaces by period analyzed.

#### *3.3. Change Detection*

Table 3 presents a summary of the cross-tabulation of data for the period between 1989 and 2001. As shown in the table, there is a predominance of persistence in all covers. There are 224.97 km2 of stable areas, equivalent to 98.90% of the total study area, and 2.51 km2 of zones that have undergone changes, corresponding to 1.10% of the total area. There is an increase in urban areas and a consequent decrease in rural areas that were occupied before the urban expansion occurred.

**Table 3.** Cross-tabulation of land cover between 1989 (columns) and 2001 (rows).


#### *3.4. Explanatory Variables*

Table 4 shows the measure of association between the explanatory variables and the land covers present in the study area. Cramer's V value fluctuates between 0.1 and 0.3. The slope is the variable that has the greatest association with the existing land cover categories (Table 4); this is because the slope affects urban expansion, as well as land use in rural areas, such as crops or the presence of natural forests. Another important explanatory variable is the elevation (DEM), which conditions urban expansion.

**Table 4.** Creamer's V values: measure of association between quantitative explanatory variables and land covers studied.


The distance to roads has an acceptable Cramer's V, corroborating the initial assumption that the presence of roads encourages urban expansion. The values of Cramer's V for distances to rivers are <0.1, probably because there is no strict regulation of urban expansion in areas near rivers.

#### *3.5. Transitional Submodels*

Table 5 shows the transition submodels, their respective variables, and the results of the calculated logistic regression. The coefficients that affect each explanatory variable are included in the logistic regression equation and the correlation between variables and transitions (ROC).

**Table 5.** Logistic regression results: Modeled transition (transition submodel), correlation (ROC), explanatory variables, and coefficients of each explanatory variable in the regression equation.


Table 6 shows the inverse relationship between the transition from nonimpervious to impervious surfaces and all the variables considered in the transition model. Chances of urban expansion are reduced when there is higher elevation, steeper terrain, and longer distances to roads. The degree of correlation between the transition studied and the explanatory variables is high, around 95%. Table 6 shows the results of the neural network application. The learning rate is low, about 1/1000, with a training and validation error of about 2/10, which is well above the acceptable error (RMS). This demonstrates the limited performance of neural networks in the present case, even though the accuracy rate is greater than 90%.

The transition probabilities for the land cover considered are included in Table 7. It can be seen that the probability of maintaining the same land use predominate, reaching almost the value of 1 in the case of nonimpervious surface and with values >1 in the case of impervious surface. As expected, the impervious surface is not likely to change to a nonimpervious surface, whereas the impervious surface is always the same.


**Table 6.** Results of the application of neural networks.

**Table 7.** Probability of transition between land uses.


Table 8 shows the confusion matrix between the map extracted from the 2013 image and the map generated by neural networks (MLP). Table 9 shows the confusion matrix between the map extracted from the 2013 image and the map generated by logistic regression (LogReg). In both tables, the comparison of the maps shows a predominance in the number of pixels that have the same thematic class. The largest errors occur when the nonimpervious surface has been modeled as an impervious surface (1497 and 1493 pixels). The errors in which the impervious surface was modeled as a nonimpervious surface are lower (289 and 285 pixels). Similarly, commission errors vary between 0.61% and 3.68%, and the maximum value corresponds to the impervious surface in both tables. The errors of omission vary between 0.12% and 16.51%, having the highest error by the commission in the transition to impervious surface

**Table 8.** Confusion matrix between the map extracted from the 2013 image and the map created through neural networks (MLP).


Table 10 shows the values of the general reliability calculated from the confusion matrices included in Tables 7 and 8, as well as the Kappa index and the correlation coefficient between the reference map of 2013 and the maps generated with logistic regression and neural networks. The map generated by logistic regression has a total reliability of 99.30%, a Kappa index of 0.8913, and a correlation coefficient of 0.8938. These values are higher than those obtained using neural networks by a very narrow margin.


**Table 9.** Confusion matrix between the map extracted from the 2013 image and the map created by logistic regression (LogReg).

**Table 10.** Validation parameters between the map extracted from the 2013 image and the maps created using logistic regression (LogReg) and neural networks (MLP).


#### *3.6. Scenario of Impervious Surfaces to 2030*

The scenario calculated for 2030 and its comparison with the existing impervious areas in 2020 is presented in Figure 6. It can be seen that, in the horizon year, important areas to the west of the city will be consolidated; this growth will be facilitated by the existence of access roads, areas with relatively flat relief, as well as the existence of small urban centers. According to this scenario, the impervious surfaces have an area of 51.53 km2.

**Figure 6.** Urban growth scenario for 2030, compared to city size in 2020.

#### *3.7. Hydrological Modeling*

The morphological characteristics of the sub-basins are presented in Table 11. It can be observed that the infiltration parameters (CN) undergo an increase as the urban area increases in each sub-basin. This increase is relatively small since, for example, in the Malacatos sub-basin, which has one with the greatest variation in CN, the CN variation reaches a value of around 3%. The variation is small because the urban area in each of the sub-basins is relatively small when compared with their total area. The sub-basins with the largest urban area (Figure 3) present a higher CN value. The concentration time is relatively short since the maximum distance that runoff must travel is related to the slope of the main channel.


**Table 11.** Characteristics of the sub-basins of the study area.

The precipitation values for different durations and return periods are indicated in Table 12. As expected, the precipitation values increase as the return period and duration increase.

**Table 12.** Precipitation values associated with each return period.


The storms included in Table 12 applied individually according to the return period, and the state of urban area expansion of the city of Loja allowed obtaining the flows included in Table 13. There is a direct relationship between the return period and the flood flows, as well as between the growth of the urban area and the flood flows for the same return period.

**Table 13.** Flood flows (m3/s) for different urbanization states and return periods.


The relationship between flow, return periods, and urban growth is presented in Figure 7, in which a high correlation between the urban area extension and the magnitude of the flows is observed for all different return periods considered.

**Figure 7.** Relationship between flood flow, return periods, and impervious surfaces in the city of Loja.

During the study period, the urban area of the city of Loja experienced a considerable increase, going from 17.68 km<sup>2</sup> in 1989 to 43.15 km<sup>2</sup> in 2020, an increase of 144.12%. The total area of the Zamora River basin is 227.48 km2, thus in 2020, the city of Loja covered only 18.97% of the total basin, while grasslands, natural forests, and shrubs covered the remaining surface. These land covers can retain surface runoff as they support infiltration, causing an opposite effect to urbanization. This may explain why the increase in flow is moderate despite the significant growth of the city.

#### *3.8. Similarities*

The behavior observed in the city of Loja has certain similarities with other urban areas around the world that experienced accelerated growth of impervious areas. Such is the case of the Alto Atoyac river basin (Oaxaca, southern Mexico), which experienced an increase in impervious surfaces of the order of 135 km<sup>2</sup> in the period between 1979 and 2013 [49]. This affected the recharge areas causing a decrease of 2.65 × 106 m3 of water infiltration into the subsoil. A similar case was observed in Addis Abab (Ethiopia) in the period between 1986 and 2016 [50], in which the impervious surfaces increased by 27%, producing a variation of 4.5 ◦C in the average surface temperature of the soil. The Pearl River delta in China [51] also experienced a very significant increase in impervious surfaces, from 390 km2 in 1988 to 4837 km2 in 2013, with the 1994–1999 period being the one with the fastest growth.

In all the cases mentioned, urban growth is related to a significant increase in the population that extends from cities to suburban areas, affecting soils that were initially covered with grasslands, forests, and agricultural areas. Although each case is different, it is possible to perceive that the increase in impervious surfaces and its effects are present in urban watersheds around the world; therefore, the proposed methodology to generate future scenarios of impervious areas can become a valuable management tool.

#### **4. Conclusions**

The NDISI satisfactorily discriminated the impervious areas in the consolidated center of the city, but in the suburban areas, an overestimation of the impervious surfaces was observed, caused by spectral confusion between impervious surfaces and bare soil (product of fallow farms, not very vigorous vegetation, and small newly opened construction areas). On the other hand, the supervised classification of Landsat images presented better discrimination of impervious areas. Therefore, the latter was elected to carry out the study of the spatiotemporal dynamics of soil impermeability in the catchment under study.

A methodology has been proposed that allows modeling future growth scenarios of impervious zones by combining the observed spatiotemporal variability, possible explanatory variables, and logistic regression models.

Slope, elevation, and proximity to highways conditioned urban growth; therefore, there was the persistence of the different land covers in the study area. The best estimate

of the change in land cover was found by logistic regression; however, neural networks performed similarly.

There is a direct relationship between the increase of impervious surfaces and the magnitude of flood flows produced by an extreme precipitation event. The basins that experience the greatest growth in impervious surfaces are those that present a greater increase in their flood flows, observing a linear relationship. If the percentage of area covered by impervious surface use is reduced compared with the areas occupied by vegetation in good condition, the increase in flood flows will be moderate.

The urbanization process directly influences the hydrological cycle, increasing impervious surfaces, reducing the infiltration capacity, and increasing the magnitude of flood flows. This must be considered in urban planning.

The increase in impervious surfaces and their effects are present in urban watersheds around the world; therefore, the proposed methodology to generate future scenarios of impervious areas can become a valuable management tool.

**Author Contributions:** Conceptualization, F.O.-V. and A.O.-P.; methodology, F.O.-V. and A.O.-P.; software, F.O.-V., A.O.-P. and M.C.; validation, A.O.-P. and M.C; writing—original draft preparation, F.O.-V. and A.O.-P.; writing—review and editing, F.O.-V. and A.O.-P.; visualization, A.O.-P. and M.C.; supervision, F.O.-V.; project administration, F.O.-V. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


#### *Article*
