Urban Land-Cover Change and Its Impact on the Ecosystem Carbon Storage in a Dryland City

Abstract

Lack of research into the complexity in urban land conversion, and paucity of observational data of soil organic carbon (SOC) beneath impervious surface area (ISA) limit our understanding of the urbanization effects on carbon (C) pools in dryland cities. Employing Landsat TM images acquired in 1990 and 2010, a hybrid classification method consisting of Linear Spectral Mixture Analysis and decision tree classification was applied to retrieve the land cover (water, ISA, greenspace, cropland, and remnant desert) of the largest dryland city in China—Urumqi. Based on vegetation carbon (VEGC) and SOC density data determined through field observations and literature reviews, we developed Urumqi’s C pool maps in 1990 and 2010, and assessed the urbanization impacts on ecosystem C. Our results showed that ISA tripled from 1990 to 2010 displacing remnant desert and cropland. The urban landscape, especially the greenspaces, became obviously fragmented. In 2010, more than 95% of the urban ecosystem C was SOC, 48% of which under the ISA. The city lost 19% of C stock from 1990 to 2010. About 82% of the ecosystem C loss was caused by the conversion of remnant desert and cropland into ISA, mainly in the northern city.

1. Introduction

Recently, processes of urbanization have become a major global driver of land use/land cover change [1]. At present, about 3%–5% of global land area has been converted into urban and developed land-use (hereafter referred to as urban) [2,3], 13%–17% of which were intensively developed [4]. With the current expansion rate, the global urban areas by 2030 would nearly triple those of the year 2000 [1]. A major finding of urban ecological research in the past decade is that urban ecosystems could account for a significant portion of terrestrial C storage at both local and regional scales [5,6,7,8,9,10,11]. Rapid urban expansion accompanied by intensive human disturbances could significantly alter ecosystem C pools [12,13]. Compared to other ecosystems, urban ecosystems are characterized by highly dynamic landscape structures and strong spatial heterogeneity in C storage. Some urban land types, like lawns or urban forests, were found to have higher C storages/sequestration capacities than their rural counterparts [5,6]. In contrast, the installation of impervious surfaces could result in significant C loss. Even in a single urban park, the soil organic C densities (SOCD hereafter) of different land cover types could differ by 10 times [14]. Therefore, an accurate assessment of the urbanization effect on C cycle requires detailed information regarding the types and locations of land-cover changes, as well as knowledge about the C densities in different pre-urban and urban land-cover types. Spatially explicit land-cover changes and C sinks/sources datasets could also help decision makers improve urban land planning and ecosystem management for the purpose of C management [15,16].

However, many previous studies tended to consider urban areas as a homogeneous land cover when assessing the urbanization effects on ecosystem C [4,17,18,19,20]. Furthermore, significant uncertainties still exist in the C density of urban land-cover types. Field observations of the SOCD in impervious surface areas (SOCD_ISA), which could account for more than 63% of the urban area [21], are particularly scarce due to the difficulties in sampling under pavement or buildings [22]. Previous studies that relied on untested assumptions of SOCD_ISA (=0 [2,23,24], 1 [25,26], or 3.3 kg·C·m⁻² [7,8]) were prone to underestimation of urban C storages and overestimation of urbanization-induced C losses [11,22]. Moreover, all of the few [11,22,27] field studies on SOCD_ISA focused in humid areas, while SOCD_ISA observations in dryland that covers more than 30% of the global land surface are still unavailable [28].

In this study, we measured the SOCD_ISA and C densities of other land-cover types in the largest dryland city of China—Urumqi city, and developed 30 m resolution multi-temporal land-cover datasets with Landsat TM images to track the land-cover changes between 1990 and 2010. Our objectives were to (1) access the amount and distribution of the ecosystem C pools of the dryland city; (2) reveal the spatial-temporal change of urban land cover from 1990 to 2010 and its impacts on the urban C dynamics; and (3) identify the major C sources/sinks and provide a scientific basis for effective C management of the city.

2. Experimental Section

2.1. Study Area

Urumqi city (43.80°N, 87.60°E), is located at the northern slope of the Tianshan Mountain Range and southern edge of Junggar Basin (Figure 1), in northwestern China. With a metropolitan area of 380 km², and a population of 2.63 million by the end of 2013, it is the largest city in the northwestern China [29]. The area has a typical arid climate, with an annual precipitation of 263.2 mm, and a mean temperature of 7.3 °C [30]. The urban soils are not well developed and the soil depth are usually <1 m. The typical soil types are solonetzs and castanozems [31].

Figure 1. Location of the study area and sampling sites with photos showing the soil profiles of the impervious surface (red box 1) and the adjacent pervious surface (red box 2) in a site. The map on the right is Landsat TM image in true color composite (bands 321(RGB)) in 2010.

Figure 1. Location of the study area and sampling sites with photos showing the soil profiles of the impervious surface (red box 1) and the adjacent pervious surface (red box 2) in a site. The map on the right is Landsat TM image in true color composite (bands 321(RGB)) in 2010.

2.2. Remote Sensing Data

Because high-resolution remote sensing images that cover the whole study area in 1990 is not available, we have to use the Mid-resolution Landsat TM images in this study. Using Landsat TM imagery also helps to avoid several problems often associated with high-resolution imagery used for urban landscape mapping: (1) shadows caused by tall objects; (2) spectral confusion between dark ISA and shadow, or between bright ISA and bare soils; and (3) high spectral variation of the same land cover due to different spectral signatures of various construction materials found in urban areas [32]. Cloud free and geometrically rectified Landsat TM 5 images (path/row: 143/30, in LT1 format, UTM coordinate system zone 45N with WGS84 datum) acquired in September 1990 and August 2010 were used in this research. Both images were composed of six reflective bands with 30 m spatial resolution and one thermal infrared band with 60 m. The six reflective bands were used to classify urban land cover, by applying a hybrid method consisting of Linear Spectral Mixture Analysis (LSMA) and decision tree classifier (DTC) (see the following Section 2.3).

Previous studies indicated that atmospheric calibration would have no effect on the fraction maps generated by endmember selection in LSMA [33,34]. In a preliminary study, we compared the fraction images developed with and without atmospheric calibration and found the two approaches did not result in obvious differences in the spatial distributions (supplementary Figures S2 and S3) or histograms (Figures S4 and S5) of the derived land-cover fractions. A paired t-test based on 1600 random pixels (0.1% of the study area) also indicated atmospheric calibration had no significant effects (p-value > 0.05) on the un-mixing results (supplementary Table S1). Therefore, atmospheric calibration was not applied in this study.

2.3. Urban Land Cover Mapping

This study aims to assess the effect of urbanization on the carbon stock of a city. Because different urban land-cover types could have very different carbon storage [35], it is important to quantify the detailed “from-to” change trajectory information in land-use change detection. To this end, our study used the post-classification approach, which has the advantage of producing complete change matrix [32,36]. Furthermore, direct use of multi-temporal spectral responses may not lead to satisfactory change detection results, since change detection based on spectral response or textures is sensitive to factors such as different atmospheric conditions, soil moisture conditions, sun elevation angle and plant phenology [37]. In order to avoid these impacts, we used the post-classification method by conducting image classification for each date of imagery separately before conducting change detection, like many previous land-use change studies [38,39,40,41].

A detailed description of the methodology was provided in Figure 2. The LSMA method was applied to generate the fraction maps of urban land covers. LSMA is based on the assumption that the spectrum of one single pixel obtained by the sensor is a linear combination of the spectral of all components within the pixel [34,42,43]. According to the Vegetation-Soil-Impervious surface (VIS) model [44], the urban landscape was treated as the combination of four components—vegetation, soil, high-albedo, and low-albedo objects [34,45]. The key to proper use of LMSA is to select the endmember for each component. First, minimum noise fraction transform (MNF) was used to reduce the six reflective bands from the Landsat images to three principle components. Then from the spectral space of the three MNF bands, the endmembers of the four land-cover components were visually selected. Finally, a full-constrained least squares solution was used to un-mix the Landsat TM images into four fractional maps.

The low albedo component included water bodies and dark impervious surfaces. The high albedo component included bright soils and bright impervious surfaces. These land-cover types were separated in two steps:

(1) The water bodies were identified using the normalized difference water index (NDWI) (Equation (1)).

$$NDWI=\left(T{M}_2-T{M}_5\right)/\left(T{M}_2+T{M}_5\right)$$

(1)

where $T{M}_2$ and $T{M}_5$ are the Landsat TM band 2 and band 5, respectively.

(2) The low-albedo and soil different index (LSDI, Equation (2)) was applied to separate bright soils from urban impervious surfaces [46].

$$LSDI=f_{Low}-f_{Soil}$$

(2)

where $f_{Low}$ and $f_{Soil}$ represent the fraction maps of the low-albedo object and soil, respectively.

Figure 2. Flow chart of urban land cover mapping. MNF—minimum noise fraction transform; V-I-S model—vegetation—impervious surface—soil model; LSMA approach—linear spectral mixture analysis approach; NDWI—normalized difference water index; GV—green vegetation; LSDI—low-albedo and soil different index.

Figure 2. Flow chart of urban land cover mapping. MNF—minimum noise fraction transform; V-I-S model—vegetation—impervious surface—soil model; LSMA approach—linear spectral mixture analysis approach; NDWI—normalized difference water index; GV—green vegetation; LSDI—low-albedo and soil different index.

Based on the fractional maps of the major land components, a decision tree classifier was applied to classify the urban land cover into four types: water, green vegetation, remnant desert/bare soil and ISA (Figure 2). We used the 1: 1,000,000 land use/cover map of China in 1990 and 2010 to separate the croplands from urban greenspaces [47].

A total of 200 sampling plots, each with a size that matched the corresponding Landsat TM pixel, were randomly selected in the study area for the purpose of accuracy assessment in 1990 and 2010 (Figure 3). Because different land-cover types had distinct land coverages in the study area (e.g., water only accounts for 0.63% of the study area, while the coverage of remnant desert/bare soil >47% in 1990), it was difficult to set a universal minimum sampling size for all land-cover types. Therefore, stratified random sampling was not suitable for this study [48]. With visual classification, the major land cover type of each sampling plot was determined based on 1-m resolution Worldview images (in 2010) and 0.5-m resolution aerial photos (in 1990) (Figure 3(a-2),(b-2)). By comparing the land-cover type of the sampling plots with that form the corresponding pixels in the land-cover maps, error matrices were generated, with both user’s and producer’s accuracies calculated [46,49,50,51] (Table 1). According to our assessment, the overall classification accuracies in 1990 and 2010 were 90% and 92%, respectively.

To reduce the bias in the land cover and land-cover change maps, we did area adjustments were conducted based on the confusion matrix using “direct” estimator (Equation (3)). Suppose the population consists of n pixels and there are k classes of land-cover types. P presents a (k × k) confusion matrix with P_ij denoting the area proportion of map class i and reference class j. A simple random sample of size n is usually used to estimate the real area instead of the whole matrix [52]. Let the (k × 1) column vector t = (t₁, t₂,…, t_k)′ denote the area proportion in each class based on the reference classification, and the (k × 1) column vector r = (r₁, r₂,…, r_k)′ denote the area proportion in each class based on the map classification. The unknown parameter of interest is t, whereas r is known. Then, the “direct” estimator is given by

$${\widehat{t}}_{dir}=(\widehat{P}{\widehat{R}}^{-1})$$

(3)

where $\widehat{P}$ was the area proportion of actual class i was misclassified into j, and $\widehat{R}=diag ({\widehat{P}}^{\text{'}}1)$ where 1 denotes a (k × 1) column vector of 1 and diag denotes a diagonal (k × k) matrix with the diagonal entries obtained from the (k × 1) column vector [53].

Figure 3. Location of the sampling plots (yellow dots) for accuracy assessments. (a) and (b) show the Landsat TM images in false color composite (bands 432 (RGB)) in 1990 and 2010, respectively; (a-1) and (b-1) are the land-cover maps developed in this study, for the year 1990 and 2010, respectively; (a-2) and (b-2) are the high-resolution reference images used in the accuracy assessments for the year 1990 and 2000, respectively. Note that our high-resolution aerial photo in 1990 only covered the center of the city.

Figure 3. Location of the sampling plots (yellow dots) for accuracy assessments. (a) and (b) show the Landsat TM images in false color composite (bands 432 (RGB)) in 1990 and 2010, respectively; (a-1) and (b-1) are the land-cover maps developed in this study, for the year 1990 and 2010, respectively; (a-2) and (b-2) are the high-resolution reference images used in the accuracy assessments for the year 1990 and 2000, respectively. Note that our high-resolution aerial photo in 1990 only covered the center of the city.

Table 1. Accuracy assessment of urban land cover classification results of Urumqi in 1990 and 2010.

Year	Classified Data	Reference Data					Reference Totals	Classified Totals	Number Correct	Producer’s Accuracy	User’s Accuracy
		Water	ISA	Urban Greenspace	Remnant Desert/Bare Soil	Cropland
1990	Water	3	0	0	0	0	3	3	3	100.0%	100.0%
	ISA	0	73	0	9	0	83	82	73	88%	89.0%
	Urban greenspace	0	0	41	0	0	42	41	41	97.6%	100.0%
	Remnant desert/bare soil	0	10	1	60	0	69	71	60	87%	84.5%
	Cropland	0	0	0	0	3	3	3	3	100.0%	100.0%
	Overall Classification Accuracy = 90% (i.e., 180/200), Overall Kappa Statistics = 0.849
2010	Water	2	0	0	0	0	2	2	2	100.0%	100.0%
	ISA	0	114	0	6	0	124	120	114	91.9%	95.0%
	Urban greenspace	0	0	30	0	0	30	30	30	100.0%	100.0%
	Remnant desert/bare soil	0	10	0	34	0	40	44	34	85.0%	77.3%
	Cropland	0	0	0	0	4	4	4	4	100.0%	100.0%
	Overall Classification Accuracy = 92% (i.e., 184/200), Overall Kappa Statistics = 0.857

Abbreviation: ISA—impervious surface area.

Table 1. Accuracy assessment of urban land cover classification results of Urumqi in 1990 and 2010.

Year	Classified Data	Reference Data					Reference Totals	Classified Totals	Number Correct	Producer’s Accuracy	User’s Accuracy
		Water	ISA	Urban Greenspace	Remnant Desert/Bare Soil	Cropland
1990	Water	3	0	0	0	0	3	3	3	100.0%	100.0%
	ISA	0	73	0	9	0	83	82	73	88%	89.0%
	Urban greenspace	0	0	41	0	0	42	41	41	97.6%	100.0%
	Remnant desert/bare soil	0	10	1	60	0	69	71	60	87%	84.5%
	Cropland	0	0	0	0	3	3	3	3	100.0%	100.0%
	Overall Classification Accuracy = 90% (i.e., 180/200), Overall Kappa Statistics = 0.849
2010	Water	2	0	0	0	0	2	2	2	100.0%	100.0%
	ISA	0	114	0	6	0	124	120	114	91.9%	95.0%
	Urban greenspace	0	0	30	0	0	30	30	30	100.0%	100.0%
	Remnant desert/bare soil	0	10	0	34	0	40	44	34	85.0%	77.3%
	Cropland	0	0	0	0	4	4	4	4	100.0%	100.0%
	Overall Classification Accuracy = 92% (i.e., 184/200), Overall Kappa Statistics = 0.857

Abbreviation: ISA—impervious surface area.

2.4. Landscape Pattern Analysis

Five landscape indexes, which were generated by the Fragstats (version 4.2), were used to quantify the spatial characteristics of land-cover change, including number of patches (NP), largest patch index (LPI), landscape shape index (LSI), contagion index (CONTAG), and Shannon diversity index (SHDI). Detailed information about the 5 indexes could be found in Table 2.

Table 2. List of landscape indexes [54].

Landscape Indexes	Equation	Description
Number of patches (NP)	$NP=n_i$	It is a simple measure of the extent of subdivision or fragmentation of the patch type.
	(ni = number of patches in the landscape of patch type (class) i.
Largest patch index (LPI)	$LPI=\frac{Max\left(a_{1,}\dots , a_n\right)}{A} \times 100$	It is a simple measure of dominance.
	(ai = area (m2) of patch i, A = total landscape area (m2).)
Landscape shape index (LSI)	$LSI=\frac{0.25E}{\sqrt{A}}$	It provides a standardized measure of total edge or edge density that adjusts for the size of the landscape. LSI increases without limit as landscape shape becomes more irregular and/or as the length of edge within the landscape increases.
	(E = total length (m) of edge in landscape, A = total landscape area (m2).)
Contagion index (CONTAG)	$CONT=\left[1+{\displaystyle {\displaystyle \sum }_{i=1}^n}{\displaystyle {\displaystyle \sum }_{j=1}^n}\frac{P_{ij}\text{ln}\left(P_{ij}\right)}{2\text{ln}\left(n\right)}\right]$	It measures the degree of clumping of different patches in a landscape.
	(n = number of the patch types, PiJ. = the proportions of different attribute adjacency types i and j.)
Shannon’s diversity index (SHDI)	$SHDI= -{\displaystyle {\displaystyle \sum }_{i=1}^m}{P}_i\text{ln}\left(P_i\right)$	It is a popular measure of diversity in community ecology, applied here to landscapes.
	(Pi = proportion of the landscape occupied by patch type (class) i.)

Table 2. List of landscape indexes [54].

Landscape Indexes	Equation	Description
Number of patches (NP)	$NP=n_i$	It is a simple measure of the extent of subdivision or fragmentation of the patch type.
	(ni = number of patches in the landscape of patch type (class) i.
Largest patch index (LPI)	$LPI=\frac{Max\left(a_{1,}\dots , a_n\right)}{A} \times 100$	It is a simple measure of dominance.
	(ai = area (m2) of patch i, A = total landscape area (m2).)
Landscape shape index (LSI)	$LSI=\frac{0.25E}{\sqrt{A}}$	It provides a standardized measure of total edge or edge density that adjusts for the size of the landscape. LSI increases without limit as landscape shape becomes more irregular and/or as the length of edge within the landscape increases.
	(E = total length (m) of edge in landscape, A = total landscape area (m2).)
Contagion index (CONTAG)	$CONT=\left[1+{\displaystyle {\displaystyle \sum }_{i=1}^n}{\displaystyle {\displaystyle \sum }_{j=1}^n}\frac{P_{ij}\text{ln}\left(P_{ij}\right)}{2\text{ln}\left(n\right)}\right]$	It measures the degree of clumping of different patches in a landscape.
	(n = number of the patch types, PiJ. = the proportions of different attribute adjacency types i and j.)
Shannon’s diversity index (SHDI)	$SHDI= -{\displaystyle {\displaystyle \sum }_{i=1}^m}{P}_i\text{ln}\left(P_i\right)$	It is a popular measure of diversity in community ecology, applied here to landscapes.
	(Pi = proportion of the landscape occupied by patch type (class) i.)

2.5. Assessing the Urbanization Effects on Ecosystem C Stock in Urumqi

The spatial maps of ecosystem C pools in the Urumqi metropolitan area were developed based on the land-cover maps in 1990 and 2010 and the VEGC and SOC density data of different land-cover types. The VEGC and SOC densities of the four major land-cover type (ISA, urban greenspace, cropland, remnant deserts/bare soil) were determined by field observation and literature review.

Using the soil pit sampling method, we collected soil samples of 0–80 cm depth from 11 ISA plots and 11 urban greenspace plots during an urban construction project, which began in March 2013. To ensure their representativeness, the 11 sampling sites were carefully selected based on location, soil types, land-use history, and land-cover stabilities (Table 3). The sampling sites scattered across the city (Figure 1), representing four major urban land-use types (industrial, commercial, transportation, and residential areas) and four common impervious surface types (street, highway, parking lot, and paved yard) (Table 3). Each of the urban greenspace plots was sampled adjacent to its paired ISA plot. The urban greenspace plots also represented four major kinds of urban green land—urban lawn, street tree, urban woodland, and residential greenspaces. Furthermore, the 11 sampling sites were evenly distributed among the two soil types (four for Solonetzs soils and four for Castanozems soils) in the study area (Table 3).

Table 3. Description of the sampling sites, and the observed BD and SOC in the ISA and PSA at 0–80 cm depth.

Site ID	Longitude (DD)	Latitude (DD)	Elevation (m)	Soil Type	Land-Use Type	ISA Type	PSA Type	ISA	PSA
								SOCDISA (kg·C·m−2)	SOCDPSA (kg·C·m−2)
1	87.65	43.84	826	Solonetzs	Industrial	Street	Street tree	5.75	9.90
2	87.48	43.87	725	Solonetzs	Industrial	Street	Street tree	3.66	2.40
3	87.56	43.86	769	Solonetzs	Institutional	Street	Street tree	2.15	10.00
4	87.53	43.86	797	Castanozems	Commercial	Parking lot	Urban woodland	6.10	8.06
5	87.54	43.84	816	Castanozems	Transportation	Highway	Urban lawn	8.03	10.30
6	87.56	43.82	853	Castanozems	Industrial	Street	Street tree	4.44	8.12
7	87.56	43.80	879	Castanozems	Transportation	Highway	Urban lawn	5.44	8.08
8	87.60	43.79	887	Unavailable	Commercial	Street	Street tree	5.28	4.03
9	87.57	43.87	765	Solonetzs	Commercial	Street	Street tree	7.64	8.29
10	87.64	43.97	591	Solonetzs	Residential	Paved backyard	Residential green space	4.35	11.72
11	87.66	43.96	613	Solonetzs	Residential	paved backyard	Residential green space	6.16	7.97

Abbreviations: ISA—impervious surface area, PSA—pervious surface area.

Table 3. Description of the sampling sites, and the observed BD and SOC in the ISA and PSA at 0–80 cm depth.

Site ID	Longitude (DD)	Latitude (DD)	Elevation (m)	Soil Type	Land-Use Type	ISA Type	PSA Type	ISA	PSA
								SOCDISA (kg·C·m−2)	SOCDPSA (kg·C·m−2)
1	87.65	43.84	826	Solonetzs	Industrial	Street	Street tree	5.75	9.90
2	87.48	43.87	725	Solonetzs	Industrial	Street	Street tree	3.66	2.40
3	87.56	43.86	769	Solonetzs	Institutional	Street	Street tree	2.15	10.00
4	87.53	43.86	797	Castanozems	Commercial	Parking lot	Urban woodland	6.10	8.06
5	87.54	43.84	816	Castanozems	Transportation	Highway	Urban lawn	8.03	10.30
6	87.56	43.82	853	Castanozems	Industrial	Street	Street tree	4.44	8.12
7	87.56	43.80	879	Castanozems	Transportation	Highway	Urban lawn	5.44	8.08
8	87.60	43.79	887	Unavailable	Commercial	Street	Street tree	5.28	4.03
9	87.57	43.87	765	Solonetzs	Commercial	Street	Street tree	7.64	8.29
10	87.64	43.97	591	Solonetzs	Residential	Paved backyard	Residential green space	4.35	11.72
11	87.66	43.96	613	Solonetzs	Residential	paved backyard	Residential green space	6.16	7.97

Abbreviations: ISA—impervious surface area, PSA—pervious surface area.

The soil samples were collected with a 100-cm³ sample ring then dried and weighted to calculate the soil bulk density (BD) using the volumetric method. Then the Mebius method was used to measure the SOCC [55]. The SOCD was calculated based on the paired BD, SOCC, and soil depth h as:

$$SOCD=SOCC\times BD\times h$$

(3)

The SOCD of the cropland and remnant desert/bare soil in the study area was estimated based on China’s second national soil survey [56]. A total of 10 remnant desert/bare soil profiles and 14 cropland soil profiles were found within the administrative boundary of Urumqi. Information related to VEGCD for each land-cover type in Urumqi was collected based on literature review [57,58].

3. Results

3.1. The Areas and Distribution of Major Land-Cover Types

In 1990, remnant desert/bare soil was the dominant land cover type (48%), and was mainly distributed in the north and southwestern part of study region, followed by the ISA that covered about 24% of the study area (Table 4). Large patches of cropland were found in the northwestern region. Clusters of urban greenspaces were found at the southern and mid-eastern edges of the ISA.

Table 4. Land cover transfer matrix of study area from 1990 to 2010 (Area: km²).

	2010	Water Body	ISA	Urban Greenspace	Remnant Desert/Bare Soil	Cropland	1990 Total
1990
Water body		2.15	0	0.24	0	0	2.39
ISA		0	90.45	0	0	0	90.45
Urban greenspace		0.05	15.32	16.37	0	0	31.74
Remnant desert/bare soil		0	81.01	26.75	72.79	0	180.55
Cropland		0	33.05	7.74	26.84	6.96	74.59
2010 Total		2.20	219.83	51.10	99.63	6.96	379.72

Abbreviation: ISA—impervious surface area.

Table 4. Land cover transfer matrix of study area from 1990 to 2010 (Area: km²).

	2010	Water Body	ISA	Urban Greenspace	Remnant Desert/Bare Soil	Cropland	1990 Total
1990
Water body		2.15	0	0.24	0	0	2.39
ISA		0	90.45	0	0	0	90.45
Urban greenspace		0.05	15.32	16.37	0	0	31.74
Remnant desert/bare soil		0	81.01	26.75	72.79	0	180.55
Cropland		0	33.05	7.74	26.84	6.96	74.59
2010 Total		2.20	219.83	51.10	99.63	6.96	379.72

Abbreviation: ISA—impervious surface area.

From 1990 to 2010, the ISAs increased by 140% and occupied about 58% of the study area in 2010; while the cropland shrank by 90%, covering less than 2% of the study area by the year 2010. It is noteworthy that compared to their spatial pattern in 1990, the greenspaces became more fragmented and scattered across the ISA in 2010 (Figure 4).

Figure 4. The land-cover types of Urumqi in 1990 (a) and 2010 (b), and the conversion from 1990 to 2010 (c and d). The circled areas indicate the two large downtown parks that disappeared in 2010. ISA indicates the impervious surface area.

Figure 4. The land-cover types of Urumqi in 1990 (a) and 2010 (b), and the conversion from 1990 to 2010 (c and d). The circled areas indicate the two large downtown parks that disappeared in 2010. ISA indicates the impervious surface area.

3.2. Changes in the Landscape Pattern from 1990 to 2010

The urban landscape became more fragmented and complex from 1990 to 2010, particularly for the urban greenspaces and remnant deserts (Figure 5). The number of patches (NP) increased significantly for all land-cover types, except for the cropland (Figure 5b). The NP of the urban greenspace more than quadrupled, and that of the remnant desert increased by more than six times during the two decades. The shape of the land patches became more complex as indicated by the increased landscape shape index (LSI) except the cropland (Figure 5a,b). Again, the LSI of the greenspaces and remnant deserts increased more significantly than the other land-cover types (Figure 5c). While their patch number increased, the patch sizes of the urban greenspaces and the remnant deserts decreased sharply (Figure 5d). The largest patch index (LPI) of the urban greenspace in 2010 was less than 57% of that in 1990.

At the city scale, the landscape dynamic from 1990 to 2010 was dominated by the changes in ISA. In contrast to the more fragmented urban greenspaces and remnant deserts, the patch size and connectivity of the ISA increased significantly (Figure 5d), leading to increased LPI and contagion index (CONTAG) in Urumqi (Figure 5a). As the result of rapid ISA expansion (Table 4), urban landscape was gradually dominated by the impervious surface, which reduced the Shannon diversity index (SHDI) of the city (Figure 5a). Meanwhile, other land-cover types became smaller and more fragmented patches that are embedded in the matrix of urban impervious surface (Figure 4).

Figure 5. Change of number of patch (NP), largest patch index (LPI), landscape shape index (LSI), contagion index (CONTAG) of and Shannon diversity index (SHDI) values of all land—cover types (a); and NP (b); LSI (c); LPI (d) of each land cover types in Urumqi from 1990 to 2010.

Figure 5. Change of number of patch (NP), largest patch index (LPI), landscape shape index (LSI), contagion index (CONTAG) of and Shannon diversity index (SHDI) values of all land—cover types (a); and NP (b); LSI (c); LPI (d) of each land cover types in Urumqi from 1990 to 2010.

3.3. The Land Conversions during 1990–2010

Urbanization in Urumqi has been characterized by the dramatic expansion of ISA and fragmentation of pervious surfaces. The majority (~63%) of the newly developed ISA was converted from the remnant desert/bare soil, followed by the converted croplands (~26%) (Table 4). Meanwhile, a large amount of cropland, accounting for 36% of the cropland area in 1990, was abandoned and changed to remnant desert/bare soil.

The spatial pattern of land conversion showed that the ISA mostly sprawled northward at the cost of remnant desert/bare soil and cropland (Figure 4c,d). Urban infill was found in the south and middle part of the city. As the landscape became increasingly dominated by the ISA, several large patches of greenspaces in the old downtown areas disappeared, and many small patches of greenspaces emerged in the new built-up areas (Figure 4c,d). Although the total greenspace area only increased moderately from 8% to 13%, intensive land-conversion took place in the greenspace areas during the two decades. Approximately half of the original greenspaces were converted to ISA, while about 68% of the greenspaces in 2010 were newly converted from other land covers (Table 4).

3.4. The Ecosystem C Dynamic and Its Spatiotemporal Pattern

According to our field measurements, impervious surface has the lowest SOCD and total ecosystem C (TOTEC) density among all land-cover types. The SOCD_ISA = 5.36 ± 0.51 kg·C·m⁻², was even lower than the SOCD of the remnant desert/bare soil (5.55 ± 0.60 kg·C·m⁻²) (Table 5). Cropland had the highest TOTEC (10.58 kg·C·m⁻²), mainly because of its high SOCD (9.94 ± 0.56 kg·C·m⁻²), which was the highest among all land-cover types in the study area. The highest vegetation C (VEGC) density was found in the urban greenspaces (1.69 ± 0.65 kg·C·m⁻²), followed by the cropland (0.64 ± 0.21 kg·C·m⁻²) (Table 5, [56,57,58]).

A total of 2.34 Mt organic C was stored in the Urumqi metropolitan area in 1990. More than 95% of the C stock was located in the soil (Table 6). Over 70% of the SOC was stored in the cropland and remnant desert/bare soil, which was mainly distributed in the north (Table 6; Figure 6a). More than 92% of the VEGC stock in Urumqi was stored in the croplands and urban greenspaces, with urban greenspace having a slightly larger VEGC stock than the cropland (Table 6). Less than 7% of the VEGC was found in the remnant desert.

Table 5. Vegetation C (VEGC), soil organic C (SOC), and total ecosystem C (TOTEC) of different land-cover types in the study region.

Land-Cover Types	VEGC ± SE (kg·C·m−2)	SOCD ± SE (kg·C·m−2)	TOTEC (kg·C·m−2)	Source
Urban impervious surface	NA	5.36 ± 0.51 (N = 11)	5.36	This study
Urban greenspaces	1.69 ± 0.65 (N = 45)	8.08 ± 0.82 (N = 11)	9.77	[57] & This study
Cropland	0.64 ± 0.21 (N = 17)	9.94 ± 0.56 (N = 14)	10.58	[56,58]
Remnant desert/bare soil	0.04 ± 0.03 (N = 10)	5.55 ± 0.60 (N = 10)	5.59	[56,58]

Table 5. Vegetation C (VEGC), soil organic C (SOC), and total ecosystem C (TOTEC) of different land-cover types in the study region.

Land-Cover Types	VEGC ± SE (kg·C·m−2)	SOCD ± SE (kg·C·m−2)	TOTEC (kg·C·m−2)	Source
Urban impervious surface	NA	5.36 ± 0.51 (N = 11)	5.36	This study
Urban greenspaces	1.69 ± 0.65 (N = 45)	8.08 ± 0.82 (N = 11)	9.77	[57] & This study
Cropland	0.64 ± 0.21 (N = 17)	9.94 ± 0.56 (N = 14)	10.58	[56,58]
Remnant desert/bare soil	0.04 ± 0.03 (N = 10)	5.55 ± 0.60 (N = 10)	5.59	[56,58]

Table 6. Ecosystem C pools (Mt) of Urumqi in 1990 and 2010, and the C dynamics from 1990 to 2010.

Land-Cover Type	1990			2010			2010–1990
	SOC	VEGC	Sum	SOC	VEGC	Sum	SOC	VEGC	Sum
ISA	0.336	0	0.336	0.877	0	0.877	0.540	0	0.540
Urban greenspace	0.263	0.055	0.318	0.371	0.078	0.449	0.108	0.023	0.131
Remnant desert/bare soil	0.992	0.007	1.000	0.496	0.004	0.500	−0.496	−0.004	−0.500
Cropland	0.644	0.041	0.686	0.070	0.004	0.074	−0.575	−0.037	−0.612
Total	2.236	0.104	2.340	1.814	0.086	1.899	−0.423	−0.018	−0.441

Abbreviations: SOC—soil organic carbon; VEGC—vegetation carbon.

Table 6. Ecosystem C pools (Mt) of Urumqi in 1990 and 2010, and the C dynamics from 1990 to 2010.

Land-Cover Type	1990			2010			2010–1990
	SOC	VEGC	Sum	SOC	VEGC	Sum	SOC	VEGC	Sum
ISA	0.336	0	0.336	0.877	0	0.877	0.540	0	0.540
Urban greenspace	0.263	0.055	0.318	0.371	0.078	0.449	0.108	0.023	0.131
Remnant desert/bare soil	0.992	0.007	1.000	0.496	0.004	0.500	−0.496	−0.004	−0.500
Cropland	0.644	0.041	0.686	0.070	0.004	0.074	−0.575	−0.037	−0.612
Total	2.236	0.104	2.340	1.814	0.086	1.899	−0.423	−0.018	−0.441

Abbreviations: SOC—soil organic carbon; VEGC—vegetation carbon.

In 2010, the TOTEC in Urumqi declined to 1.90 Mt. Like in 1990, over 95% of the C stock was stored in soils. However, soils beneath urban ISA become the main SOC pool, accounting for the 48% of the SOC in the city, larger than the total SOC stored in the croplands and the remnant desert/bare soil (Table 6; Figure 6c). In 2010, about 90% of the VEGC in Urumqi was found in the urban greenspaces (Table 6, Figure 6d).

Land conversion resulted in 0.44 Mt TOTEC loss during the period of 1990–2010, primarily in the northern, middle, and southern parts of the city. Over 95% of the TOTEC loss was from the SOC pool. Nearly 19% of the SOC stock in 1990 has been lost during 1990–2010. ISA expansion was the major C release process, accounting for about 89% of TOTEC C loss (Figure 6). Small patches of carbon sinks were found in the newly developed urban greenspaces that accounted for 20% of TOTEC increase during 1990–2010 (Table 6, Figure 6g).

Figure 6. Spatial patterns of the ecosystem C pools in 1990 and 2010 (a–f) and the C dynamics in response to land-cover change from 1990 to 2010 in Urumqi (g).

Figure 6. Spatial patterns of the ecosystem C pools in 1990 and 2010 (a–f) and the C dynamics in response to land-cover change from 1990 to 2010 in Urumqi (g).

4. Discussion

Land cover change in Urumqi from 1990 to 2010. Since urban ecosystem C storage is related to pre-urban land-cover types as well as the current land cover, knowledge related to land-cover change patterns is important for our understanding of urbanization effects on the C cycle and locating of C sinks/sources in urbanized areas [59]. The land-cover change in Urumqi from 1990 to 2010 has been characterized by dramatic ISA expansion at the cost of cropland and remnant desert/bare soil (Table 4). Although the urban greenspace increased slightly, it, like all other pervious land-cover types, experienced dramatic fragmentation (Figure 5). As small patches of pervious surfaces were being separated by the ever-expanding impervious surfaces, most native plants and animals that lived in these “urban green-islands” would be isolated from each other. Meanwhile, the increased contagion index, the shrinking pervious patch sizes, and the more complex patch shapes (Figure 5) indicated intensified interactions between the ecosystems in urban and the paved areas (buildings, roads, etc.) where humans live and work. In other words, the environments (climate, atmosphere, ground water, etc.) of the urban ecosystems became more susceptible to human disturbances [60,61]. The intensified human disturbances together with the habitat loss and isolation threatened the stability of the urban ecosystem functions (e.g., C sequestration) as well as the ecosystem structure (e.g., native species abundances) in Urumqi [59]. Because the city is located in an arid climate, and because the ecosystem is relatively fragile compared to other eco-zones [62], rapid urbanization in this area may have more serious environmental/ecological consequences than in other areas.

The effect of urbanization on ecosystem C storage. Our study showed that urbanization from 1990 to 2010 resulted in about 1.23 kg·C·m⁻² loss in Urumqi. The urbanization effects in Urumqi were weaker than those found in the southern US (−2.6 kg·C·m⁻²) [63], United Kingdom and the Northern Ireland (−8 ± 4 kg·C·m⁻²) [25]. This is because the dryland ecosystem had extremely low VEGC storage (0.04 kg·C·m⁻²), and thus suffered less VEGC loss (compared to other regions) when the biomass was removed due to ISA expansion (Table 4 indicates about 63% of the newly developed ISA was converted from remnant deserts). More importantly, we found that the SOCD of the impervious surfaces (5.36 ± 0.51 kg·C·m⁻²) was only slightly lower than that of the remnant deserts (5.55 ± 0.60 kg·C·m⁻²), meaning the SOC loss due to remnant deserts to ISA conversion was also small. Previous studies that assumed the ISA conversion would remove 100% [25] or 50% [63] of the SOC might overestimate the negative effect of urbanization on ecosystem C balance [8,22]. Furthermore, we found that about 15% of the disappeared remnant deserts have been converted to urban greenspaces, which had much higher VEGC (1.69 kg·C·m⁻² vs. 0.04 kg·C·m⁻²) and SOC (8.08 kg·C·m⁻² vs. 5.55 kg·C·m⁻²) than the native deserts (Table 5). Management activities such as irrigation and fertilization enable urban greenspace to sequestrate more C than natural ecosystems in dryland [6,63]. The carbon sinks from the urban greenspaces could partially offset the C loss during urbanization. However, it should be noted that any gains in C by urban greenspace management were at the cost of increasing water, fertilizer, and fossil fuel usages.

Locating the C sinks and sources in Urumqi. In the face of rapid global urbanization in the 21st century, the potential for C management in urban and developed areas has drawn attention from both ecologists and decision makers [64]. Effective C management requires detailed information about the locations and magnitudes of the C sinks/sources. Our results showed large patches of intensive C sources in the northern and northwestern outskirts of Urumqi, where many croplands were converted to ISA as the built-up areas expanding northward. Although the northeastern expansion of the ISA seemed to be more intensive than the northwestern expansion (Figure 4d), fewer strong C sources were found in the northeast, because most ISA in this area were converted from remnant deserts that had low C density (Figure 4c, Table 4). Furthermore, the newly developed built-up areas in the northeast had relatively high urban greenspace coverage (Figure 4b), which acted as C sinks (Figure 6g). Meanwhile, small patches of C sources emerged across the built-up areas, mainly due to the urban infill that converted the greenspaces to ISA (Figure 4d). The conversions of two large greenspaces in the mid-eastern and southeastern downtown areas to ISA from 1990 to 2010 (Figure 4c,d) especially created two strong C sources (Figure 6g).

Based on these results, we recommend the municipal government to (1) limit northwestern urban expansion and direct the expansion to the desert areas in the northeast, trying to avoid conversion of cropland to ISA; (2) protect urban greenspaces and recover the large parks in the downtown area. Although the local government did a good job in maintaining the overall greenspace coverage in Urumqi, we found that the landscape pattern of the urban greenspaces has become highly fragmented during 1990–2010. The environmental services (e.g., climate regulation) and ecological functions (e.g., wildlife habitats) provided by the few large urban parks in downtown could not be fully compensated by many small patches of greenspaces that scattered across the metropolitan area [65].

Uncertainties. When using the post-classification method, errors in the classification maps could influence the accuracy of land-cover change detection, and affect the carbon assessment results in this study. According to the overall classification accuracies in the 1990 (90%) and 2010 (92%), the change detection error was 17%, or 64 km² in the study area (Table 1). Our results showed the land-cover change area from 1990 to 2010 was 191 km² (Table 4). Therefore, the change detection error might have influenced up to one-third of the actual changed area. We used the direct estimation adjustment to reduce the impacts of land-cover misclassification on our estimation of the carbon stock. Comparing the pre-adjustment and the post-adjustment results, we found the urban carbon stock could be overestimated by 2.7% in 1990 and be underestimated by 0.3% in 2010, due to the land-cover misclassifications. The relatively low impacts from the land-cover misclassifications on the carbon stock estimations were due to the similar carbon densities between the ISA (5.36 kg·C·m⁻²) and remnant desert (5.59 kg·C·m⁻²) lands (Table 5), which were the major sources of the overall classification error (Table 1). When applied in other regions with greater differences in the carbon densities among the converted land-cover types, the change detection errors from post-classification method may severely affect the carbon dynamic estimations. In such cases, some other multi-temporal classification change detection (i.e., pre-classification methodology) may be more suitable, especially if there is no need to quantify the detailed “from-to” change trajectory information.

Another major source of uncertainty is related to the estimated carbon density values for the urban ecosystems [66]. Ideally, a randomized sampling with large number of sampling plots is required for a confident field investigation on carbon density. In reality, however, our sampling of ISA was limited to the 11 excavation sites during an urban construction project. Although the location of these sampling sites generally represented the major ISA types, land-use types, and soil types of the study area, the approach was not a true random sampling approach. Other SOC_ISA studies faced similar problems, due to the inaccessibility of impervious-covered soil. Moreover, the SOCD of the cropland and the remnant deserts in the Urumqi metropolitan area were estimated based on the soil data from a national soil survey that focused on rural area. The underlying assumption was that urbanization-induced environmental changes might not significantly alter the SOCD of the croplands and remnant deserts that were mainly distributed in the city outskirts. Although this was a popular assumption in many previous urban carbon studies [7,8,25], its validity in the dryland city is yet to be tested. In general, more observation data from the SOC and VEGC of urban ecosystems were needed to improve statistical confidence and reduce uncertainties.

5. Conclusions

Land conversion related to urbanization has profound effects on regional ecosystems not only because of the occupation of the natural ecosystem but also due to the complex landscape patterns and intensive human disturbances. Our study investigated the effects of urbanization from 1990 to 2010 on the landscape characteristics and carbon stocks of the largest dryland city of China (i.e., Urumqi city). Urbanization in Urumqi has been characterized by the dramatic expansion of ISA, at the cost of remnant desert and cropland. It is noteworthy that compared to their spatial pattern in 1990, the greenspaces became more fragmented in 2010. These changes in land-cover composition and landscape structure threatened the native plants and animals that mainly live in the “urban green-islands” in the metropolitan area. Meanwhile, as large areas of croplands that had high carbon density were converted to ISA with relatively low carbon density, the metropolitan area acted as a carbon source, esp. at the northwestern outskirts. The carbon sources from urban areas contribute to global atmospheric CO₂ enrichment and may affect the global climate system.