Spatially Explicit Subpixel-Based Study On The Expansion of Articial Impervious Surfaces And Its’ Impacts On Soil Organic Carbon

: Precise spatiotemporal datasets of artificial impervious surfaces (AISs) are essential for 9 evaluating urbanization processes and associated soil organic carbon (SOC) dynamics. However, s patially 10 explicit studies on SOC stocks based on high-quality AIS data remain deficient, which affects the accuracy 11 of urban SOC budgets. In this study, we used 30-m Landsat images and a subpixel-based model to 12 accurately evaluate and quantify the annual AIS of Kaifeng, an ancient city in China that experienced 13 intensive urbanization from 2000 to 2020. Soil organic carbon (SOC) dynamics were further estimated and 14 spatially exhibited based on the SOC densities (SOCD) of different land covers observed in the field. Our 15 results demonstrate that Kaifeng experienced drastic AIS expansion from 2000 – 2020, both in total area (an 16 increase of ~154.35%) and density (described by mean AIS abundance, 0.56 vs. 0.72). Spatially, AIS mainly 17 sprawled to the west, and infilling was observed in the old town. Moreover, the expansion of AIS in Kaifeng 18 has resulted in a total of 0.08 Tg of SOC loss over the past 20 years, and the study area has acted as a clear 19 carbon source. The greatest SOC losses occurred during 2010 — 2015, mainly in the west — with >30% 20 (~0.024 Tg) of the total loss occurring between 2010 and 2015. This study provides new insights into urban 21 growth through the mapping of growth patterns in terms of both outward sprawl and infill. We also provide a novel means of presenting the spatial patterns of urbanization-induced SOC dynamics using subpixel AIS maps.


Introduction 27
Urbanization has become the main theme of global land change and is the foremost factor affecting the 28 carbon cycle at multiple scales (Hutyra et al., 2014; Zhu et al., 2019). The most direct evidence of 29 urbanization is the conversion of agriculture and/or natural lands into artificial environments and the sealing 30 soils with artificial impervious surfaces (AISs). Generally, AIS is composed of materials that prevent the 31 natural infiltration of water into soils and include building roofs, cement squares, and road surfaces (Zhu et  scattering. To facilitate the selection of endmembers, the minimum noise fraction (MNF) method was used to 118 determine the intrinsic noise of each image, and to ensure that the primary information was concentrated in 119 the first three or four bands. 120 Endmember collection is a key step in the subpixel approach employed here. An endmember, which is 121 distinguished from the mixed pixels, represents pixels that contain only the spectral information of one land 122 cover type. Ideally, endmembers of a certain land cover are distributed at the top of the triangle generated by 123 the different MNF bands. With the support of high-resolution datasets (collected from Google Earth v. 7.3.3, 124 Google LLC, USA), the endmembers from four groups, which represented the four land-cover types of high-125 albedo objects, low-albedo objects, green vegetation, and bare soil, were collected from the two-dimensional 126 scatter plots generated by the first three bands (Fig. S1). 127 A linear spectral mixture analysis (LSMA) model was employed to generate urban fractional land-cover 128 maps. This is one of the main subpixel-based methods and is commonly used for extracting AIS from 129 medium-resolution remote sensing data (Wang and Li, 2019). The LSMA model assumes that the spectrum of 130 a single-pixel captured by a sensor is a linear combination of all components within that pixel (Equation 1, 131 where i is the number of bands used, k is the number of endmembers, such that k = 1, 2, …, n, R i is the 133 reflectance of band i, which may contain more than one endmember, f k is the abundance of endmember k 134 within a pixel, which represents the proportion of AIS within a single pixel and indicates the AIS density, R ik 135 is the spectral reflectance of endmember k in a single pixel on band i, and is the error for band i. A fully 136 constrained least-squares solution was then applied to unmix the remote sensing data into four fractional 137 maps. The spatial AIS data were generated by taking the sum of fractional maps of the high-albedo and low-138 albedo objects, according to the methods of Lu and Weng (Lu and Weng, 2004 and 2020, based on the aggregated AIS data (Fig. S3). For each year, these samples were geographically 146 linked to the corresponding high-resolution images and the reference AIS was digitized. The percentage of 6 AIS in each plot then was calculated (i.e., AIS/14400 m 2 , e.g., Fig. S3f). We used the root mean square error 148 (RMSE, Equation 2) and Pearson's correlation coefficient (R, Equation 3) to evaluate the accuracy of the data. 149 where ., is the visually interpreted AIS abundance in plot i, ., is the corresponding estimated AIS 151 value in plot I, .
̅̅̅̅̅̅ ̅̅̅̅̅̅̅ is the average value of visually interpreted AIS abundance and the corresponding 152 estimated AIS, and n is the number of sampling plots (n = 30). 153 In this study, we focused only on urban expansion, while excluding the renewal process. Thus, AIS 154 expansion was assessed based on the assumption that urban growth was irreversible. For the same pixel, if the 155 value in the later periods was lower than the former, then the value of the former was given to the latter. The 156 net gain in AIS was then used to reveal the expansion of AIS, which can be defined by Equation 4 as follows: 157 where and are the AIS areas at the beginning and end of the study period, respectively. The 158 intensity of AIS expansion was then quantified as follows (Equation 5): 159 where represents the annual growth rate of the AIS area, and are the same as in Equation (4),  160 and n represents the time period. The rate of change in AIS was estimated by the slope, K (Equation 6). 161 Where AIS represents the total AIS area in a certain year and i is the number of years. Finally, we randomly 162 extracted the values of 5000 points from the annual AIS map to analyze and thereby evaluate the changes in 163 the urban form of Kaifeng from 2000 to 2020. 164 7 We used field SOCD data collected from different land covers across Kaifeng to quantify SOC stocks in 166 each year and to evaluate SOC dynamics (Fig. 1) (Table S3) were 168 collected in 2009. Under the hypothesis that soil sealed by AIS is stable, the spatiotemporal pattern of SOC AIS 169 from 2000 to 2020 (at 5-year intervals) was visualized and quantified based on spatially explicit AIS data. It 170 must be noted that no change trajectory could be generated from the fractional AIS images; therefore, this 171 study is mainly focused on SOC dynamics due to the expanded AIS occupying other land covers. Here, the 172 measured SOCD OPEN was simplified as follows (Equations 7 and 8): 173 SOCD OPEN = n 1 SOCD 1 +n 2 SOCD 2 +⋯+n i SOCD k n 1 +n 2 +⋯+n i (7)

Assessing SOC AIS dynamics and locating carbon sources/sinks
where ̅ represents the weighted mean SOCD OPEN , represents the SOCD of land-cover type k, f is the 174 number of sites taken by this land cover, n is the sum of the sample points taken in the open soils, and 175 SOCD AIS is the mean value estimated based on the soil sampling points located along roads and at buildings 176 (Table S3): 177 where ̅ represents the mean SOCD AIS and denotes the observed SOCD AIS . 178 The SOC dynamics induced by AIS expansion were calculated similarly to the theory underlying LMSA, 179 wherein the total SOC stock of a certain pixel was composed of SOCD AIS and SOCD OPEN (i.e., urban green 180 space, bare lands, and croplands). The total SOC of this pixel could then be expressed as follows: 181 10 could be simplified as follows: 186 In this study, the reduced area of open soils was equal to the increase in AIS. Therefore, the SOC 187 dynamics could be simplified once more as: 188 . The SOC dynamics were then 189 spatially illustrated based on pixel-based results utilizing spatial analysis methods.
increasing from 0.56 in 2000 to 0.72 in 2010. After 2010, the rate of increase slowed and the AIS abundance 204 remained stable, though a linear fit to the mean values shows a declining trend since 2010 (Fig. 4b). 205 According to the abundance of AIS (Fig. S4) (Fig. 5a-d). Specifically, the net gain in the AIS area increased from 15.56 km 2 in 2005 208 (equal to 26.9% of AIS area in 2000) to 16.30 km 2 in 2010 (equal to 22.2% of AIS area in 2005) (Fig. 5f, g). 209 The most intensive expansion accrued between 2010 and 2015, during which the areas of AIS expanded by 210 ~27.62 km 2 in 2015, which was more than twice the AIS in 2000 (Fig. 5h). The spatial patterns of AIS 211 expansion could be drawn as infilled old towns (main body of metropolitan area in 2000 in the east of 212 Kaifeng), with most sprawl occurring to the west from 2000-2015 ( Fig. 5f-h). The intensity and extent of the 213 increase in AIS decreased after 2015 (Fig. 5e, i). The newly developed AIS was 15.03 km 2 , which was less 214 than the increase observed from 2000-2005 (Fig. 5f, i). Spatial expansion was predominantly characterized 215 by an infilling growth pattern based on the extent of expansion during the preceding period, which occurred 216 throughout the main urban area. Moreover, the newly developed AIS was scattered and spatially 217 discontinuous in the northwest, north, and southwest of the built-up area (Fig. 5j). 218

SOC Dynamics caused by the Extension of AIS 219
A total of 0.51 Tg (1 Tg = 10 12 g) of SOC was stored beneath the AIS in 2000, which had increased to 220 1.17 Tg by 2020 and had more than doubled since 2000 ( Fig. 6a-f). It should be noted that the period from 221 2010-2015 exhibited the largest growth in SOC AIS (~0.25 Tg), accounting for ~31% of that in 2010. Gains in 222 SOC AIS were spatially consistent with the overall expansion of AIS, and were mainly concentrated in the 223 western part of the metropolitan area Fig. 6i). However, the data indicated that SOCD AIS was lower than 224 SOCD OPEN (Table S2). Gains in the SOC AIS stock also indicated SOC was lost during AIS expansion. As 225 shown in Fig. 6 (Fig. 6f). These two periods 229 were dominated by slight losses in SOC that occurred throughout the study area, and small patches of strong 230 SOC sources in the northwest (Fig. 6a, b). From 2010-2015, severe SOC loss occurred over a large area ( Fig. 10 After 2015, the SOC loss driven by AIS installation was only 0.008 Tg (Fig. 6Error! Reference source not  234 found.f), which was the lowest loss among each period, and was mainly distributed in the northwest of the 235 city (Fig. 6d). 236

Changing characteristics of AIS in Kaifeng 238
As the main type of land cover, it is crucial to understand the magnitude and spatial distribution of 239 changes in AIS. The spatiotemporal patterns of AIS are closely related to urbanization-induced environmental 240 problems, particularly when evaluating soil-related ecological problems, such as SOC losses in our study. 241 Kaifeng has experienced dramatic growth since 2000 and has primarily sprawled to the west (Fig. 4). There 242 are many reasons for this spatial change. First, Kaifeng was the old provincial capital of Henan, and the main 243 part of the built-up area consisted of a high density of low-height buildings. This is also reflected in (Fig. 5a- Some remaining uncertainties in this work should be addressed in the future. The first stems from the 328 estimated SOCD AIS . Generally, many randomized soil sample plots are required to calculate a confident 329 SOCD value. However, it is difficult to collect soil samples beneath AIS because of the constraints of urban 330 management regimes, which is a common issue in most studies concerned with soil properties under AIS 331 (Table S1). We recognize that it is unconvincing to use limited data to assess the dynamics induced by AIS 332 expansion. However, the essence of AIS expansion is the conversion of soil with a higher SOCD into soil 333 with a lower organic carbon density beneath the AIS, which results in a considerable amount of SOC loss 334