Assessing the Impacts of Urbanization-Associated Land Use/Cover Change on Land Surface Temperature and Surface Moisture: A Case Study in the Midwestern United States

Abstract

Urbanization-associated land use and land cover (LULC) changes lead to modifications of surface microclimatic and hydrological conditions, including the formation of urban heat islands and changes in surface runoff pattern. The goal of the paper is to investigate the changes of biophysical variables due to urbanization induced LULC changes in Indianapolis, USA, from 2001 to 2006. The biophysical parameters analyzed included Land Surface Temperature (LST), fractional vegetation cover, Normalized Difference Water Index (NDWI), impervious fractions evaporative fraction, and soil moisture. Land cover classification and changes and impervious fractions were obtained from the National Land Cover Database of 2001 and 2006. The Temperature-Vegetation Index (TVX) space was created to analyze how these satellite-derived biophysical parameters change during urbanization. The results showed that the general trend of pixel migration in response to the LULC changes was from the areas of low temperature, dense vegetation cover, and high surface moisture conditions to the areas of high temperature, sparse vegetation cover, and low surface moisture condition in the TVX space. Analyses of the T-soil moisture and T-NDWI spaces revealed similar changed patterns. The rate of change in LST, vegetation cover, and moisture varied with LULC type and percent imperviousness. Compared to conversion from cultivated to residential land, the change from forest to commercial land altered LST and moisture more intensively. Compared to the area changed from cultivated to residential, the area changed from forest to commercial altered 48% more in fractional vegetation cover, 71% more in LST, and 15% more in soil moisture Soil moisture and NDWI were both tested as measures of surface moisture in the urban areas. NDWI was proven to be a useful measure of vegetation liquid water and was more sensitive to the land cover changes comparing to soil moisture. From a change forest to commercial land, the mean soil moisture changed 17%, while the mean NDWI changed 90%.

1. Introduction

Cities in America have experienced urbanization in a variety of types, shapes, and sizes since the 19th century. With the introduction of transportation techniques, public and private transportation revolutions and a newly developed cultural value of living reshaped the spatial distribution of cities [1]. The physical and socioeconomic distinctions between urban and suburban areas became blurry. The common scenario of urbanization was that the commercial land spread along major highways from the center of the cities to the suburbs, and the residential land replaced the farmland at the periphery [2]. During the urbanization process, natural land cover, such as vegetation, exposed soil, and standing water were replaced with anthropogenic materials such as concrete, metal, and asphalt. Along with the change in land use and land cover (LULC) were modifications in surface energy and water balance, which resulted in the urban heat island (UHI) phenomenon and the unique characteristics of urban runoff.

The relationship between UHI effect and LULC changes has been examined to understand the impact of LULC changes on surface thermal properties [3]. In Weng et al. [4], the relationship of Land Surface Temperature (LST) and vegetation abundance were investigated through various scales. The results showed a stronger relationship between LST and vegetation fraction than NDVI in different spatial resolutions and different land use types. The greatest negative correlation between surface temperature and the vegetation abundance indicator occurred at the resolution level of 120 m, which was roughly the length of a city block. A city block is a basic unit in an urban area, and the characteristics of the landscape and the level of energy exchange is relatively similar in a block. In the examination of the spatial pattern of surface temperature in transects, a higher spectral variability occurs when the proportion of different land cover types is distributed more evenly, lower spectral variability occurs when less land cover types were found in a transect or one land cover type occupied the majority of the surface. Similarly, Yuan and Bauer [5] investigated the relationship between LST and percent impervious surface areas, and found a strong linear relationship in all seasons. Therefore, the percent impervious surfaces provide a complementary metric to the Normalized Difference Vegetation Index (NDVI) for analyzing LST over all seasons. Weng and Lu [6] demonstrated that the sub-pixel technique is an effective approach to classify urban LULC and characterize the impact of changes of LULC on LST. Sub-pixel technique is also a solution for mixed pixel problems in 30m resolution images. Weng and Lu [7] further tested the sub-pixel technique and vegetation-impervious surface-soil models to characterize urban landscapes, and proved that the combination of the two is an alternative approach to quantify the spatial and temporal changes of urban landscape composition. Carlson [8] formed the Temperature-Vegetation Index (TVX) space by plotting LST and vegetation fractions and showed the variation of surface moisture availability in TVX space. The advantages of the TVX method are that the surface moisture and evapotranspiration can be generated easily, and no ancillary atmospheric or surface data, or land surface model is needed. In addition, since the pixel distribution itself is used to set the conditions, it is not sensitive to atmospheric correction and surface parameters. The disadvantage is that the determination of the warm edge is subjective. In Amiri et al. [9], the TVX space was constructed to demonstrate the pixel trajectory of pixels due to LULC changes. The results showed that the TVX method could be applied to monitor changes in biophysical parameters due to LULC changes. Carlson and Arthur [2] illustrated how LST, surface moisture availability, fraction impervious surface, and urban-induced surface runoff responded to urbanization in a TVX space. Jiang and Islam [10] used TVX space and an extension of the Priestley-Taylor equation to estimate surface evaporation over heterogeneous areas. The results showed that this approach is more reliable and easily applicable for evaporation estimation where ground based data are not available. Owen et al. [11] and Carlson and Arthur [2] calculated soil moisture by the Soil-Vegetation-Atmosphere Transfer (SVAT) model, relating the moisture availability to the TVX space to assess the impact of urbanization on microclimate and hydrology. Sun and Kafatos [12] further suggested that the TVX space may not hold true in cold seasons based on a study in the South Great Plains.

Despite the progress mentioned above, the understanding of surface moisture change caused by urban LULC changes is still limited. Urban runoff was found to be highly dependent on regional climate, season, and preceding moisture condition [13]. Thus, the applicability of using runoff to assess the surface moisture condition may be limited in arid and semi-arid areas, and in the seasons with little precipitation. High thermal conductivity and high thermal inertia of wet soil made soil moisture a significant component in the variation of LST [14]. The effectiveness of using soil moisture that derived from the surface energy balance models to evaluate the urbanization impact deserves further exploration, because various vegetated and soil areas were replaced by impervious surfaces during the urbanization process. Additionally, Gao [15] suggested that Normalized Difference Water Index (NDWI) may be used as a measure of vegetation liquid water. This finding provided an alternative way to assess surface moisture in urban areas.

In this study, the performance of the biophysical parameters such as vegetation fraction, soil moisture, NDWI, and LST to monitor the impact of urbanization on land surface characteristics were examined by employing three representative areas in Indianapolis, United States, between 2001 and 2006. Soil moisture refers to near surface soil moisture estimation and correlates well with 0-5cm depth soil moisture field data. It is volumetric soil moisture, which is the ratio of the volume of water and the total volume of soil, water, and air. Surface moisture does not equal soil moisture. It is the wetness availability over all types of surfaces. In this study, surface moisture was represented by soil moisture and NDWI. In the TVX space, when there is less vegetation, there is higher land surface temperature, and vice versa. At a fixed vegetation level, when there is high soil moisture, there is lower land surface temperature. Evaporation fraction is related to both soil moisture and surface temperature: higher moisture and higher temperature results in higher evaporation.

Parameters such as soil moisture and LST have the ability to characterize changes induced by urbanization-associated LULC changes. This type of impact assessment is significant because the focus of the existing literature has been exclusive of whether and how urbanization has a clear positive or negative effect on the environment [16]. Based on the land cover data from the National Land Cover Database (NLCD), three representative types of LULC changes in the Midwest were selected, including: (1) from cultivated to residential land; (2) from forest to commercial land; and (3) from open space to commercial land. The responses of LST, vegetation cover, and moisture availability to these changes were investigated by using the TVX method. First, surface temperature and soil moisture parameters were derived in the TVX space; then, the changes in surface temperature and moisture and their relationship between different land cover classes were examined in the TVX space, LST-soil moisture space, and LST-NDWI space, respectively. By performing these tasks, this study attempts to address the following research question: how LCLU changes within urban areas would alter surface temperatures, vegetation cover, and surface moisture conditions.

2. Data and Methodology

2.1. Study Area

The study area, Marion County, Indiana, was located in the Midwestern part of the USA (Figure 1). The county seat is Indianapolis, the capital and largest city in the state of Indiana. The distribution of fractional vegetation cover for the whole study area is lower vegetation cover in central urban areas, and higher vegetation cover outside of the circle of interstate highway in agricultural fields and forests. The city is located on a flat plain, which makes it possible for the city to expand in all directions. According to the NLCD 2006, the conversion to developed land from 2001 to 2006 mainly took place in the suburban areas, specifically, between the circle of Interstate Highway 465 and the county boundary, with the largest changes in the southern and eastern fringes. Sparse land cover change in terms of density took place in the urban areas within the circle of Highway 465. According to the NLCD 2001-2006 land cover change data, 3.65% of the total land cover was changed from 2001 to 2006 in Indianapolis, IN, among which, 51% was changed from cultivated land to developed land, 31% was changed from a lower level of developed land to a higher level of developed land, and 3.6% was changed from forest to developed. Three small study areas representing the three typical urbanization-associated LULC changes were selected. The resultant land cover type was “developed” in the three small areas in 2006. In terms of intensity, the land cover types were developed with high and medium intensity in Areas 2 and 3. For developed open space, low intensity and medium intensity development occurred in Area 1.

2.2. Methodology

The workflow is shown in Figure 2. First, the Landsat TM data was acquired and was pre-processed to correct the atmospheric effects. Second, NDVI and NDWI were calculated. Third, LST was retrieved by radiative transfer equation, and fractional vegetation cover was computed by Linear Spectral Mixture Analysis (LSMA). Fourth, the TVX space was formed by the scatter plot of LST and fractional vegetation cover, and soil moisture was calculated with air temperature. LST, fractional vegetation cover, soil moisture, NDWI, and percent of imperviousness was compared between the images acquired in 2001 and 2006, and the impacts of LULC change due to urbanization were analyzed.

Figure 1. Three sample areas for assessing the impact of urbanization in the City of Indianapolis, USA. The vegetation fraction was produced by the Landsat TM image that acquired on 16 June 2001. (a) The whole study area, (b) Area 1 is characterized by conversions from cultivated to residential lands, (c) Area 2 and Area 3 represent changes from forest to commercial, and from open area to commercial respectively.

Figure 1. Three sample areas for assessing the impact of urbanization in the City of Indianapolis, USA. The vegetation fraction was produced by the Landsat TM image that acquired on 16 June 2001. (a) The whole study area, (b) Area 1 is characterized by conversions from cultivated to residential lands, (c) Area 2 and Area 3 represent changes from forest to commercial, and from open area to commercial respectively.

2.2.1. Data Collection and Pre-processing

Two Landsat TM images were acquired on 17 June 2001 and 1 July 2006, respectively. Land cover, percent developed imperviousness in 2001 and 2006, and land cover change data were obtained from the NLCD. The land cover change data contained only the pixels identified as changed between NLCD2001 Land Cover Version 2.0 and NLCD2006. The weather conditions of two weeks before the image acquisition day were compared: the total precipitation for the two weeks were 84.1 mm and 83.5 mm, the average of daily maximum temperature was 24.16 ºC and 28.02 ºC , and the average daily minimum temperature was 13.01 ºC and 17.39 ºC, respectively. Thus it was assumed that antecedent soil moisture conditions on the two image acquisition days were similar. This assumption was made based not only on the weather conditions, but also considered the characteristics of urban surface and evaporation. Since our study area is in the urban setting, evaporation over impervious surfaces did not follow the evaporation characteristics that have been found in soil and vegetated areas. Ramamurthy and Bou-Zeid [17] suggested that substantial contribution from impervious surfaces to urban evaporation usually occurs in the first 48 hours after precipitation. Since the amount of precipitation was below 1mm within 48 hours prior to both image acquisition dates, the urban evaporation conditions of the two dates were assumed to be similar. Furthermore, since impervious surfaces increase urban surface runoff, most precipitation in the urban areas would become runoff instead of evaporation. Therefore, the impact of precipitation on urban evaporation cannot be as significant as that on natural surfaces.

Figure 2. The flowchart of the study. The blue polygons with white text are the resultant biophysical parameters that were compared between 2001 and 2006.

Figure 2. The flowchart of the study. The blue polygons with white text are the resultant biophysical parameters that were compared between 2001 and 2006.

Atmospheric correction was applied prior to image processing. The at-sensor radiance was converted from calibrated Digital Number (DN), according to the referenced values from Chander et al. [18]:

$$L_{at-sensor}=G_{rescale}\times Q_{cal}+B_{rescale}$$

(1)

where ${\text{L}}_{\text{at}-\text{sensor}}$ is the at-sensor radiance, ${\text{G}}_{\text{rescale}}$ is band-specific rescaling gain factor, ${\text{B}}_{\text{rescale}}$ is band-specific rescaling bias factor, and $Q_{cal}$ is the quantized calibrated pixel value (DN).

The method for atmospheric correction for bands 1-5 and 7 was adopted from Song et al. [19]. This simplified dark object subtraction method assumed no atmospheric transmittance loss and no diffuse downward radiation at the surface. The radiative transfer equation can be written as:

$$L_{at-sensor}= L_p+ \frac{\rho F_d{T}_v}{\pi \left(1-s\rho \right)}$$

(2)

where $L_p$ is the path radiance, $F_d$ is the irradiance received at the surface, $T_v$ is the atmospheric transmittance from the target toward the sensor, $s$ is the fraction of the upward radiation back-scattered by the atmosphere to the surface, and $\rho$ is the surface reflectance. The variable $s$ was neglected according to Song et al. [19]. The equation for calculating $\rho$ can be expressed as:

$$\rho = \frac{\pi \left(L_{at-sensor}- L_p\right)}{T_v(E_0(cos\left({\theta }_z\right)T_z+E_{down})}$$

(3)

where $E_{down}$ is the downwelling diffuse irradiance (W/m²/µm), $E_0$ is the exoatmospheric solar constant (W/m²/µm), $T_z$ the atmospheric transmittance in the illumination direction, and ${\theta }_z$ the solar zenith angle.

Assuming 1% surface reflectance for the dark objects, the path radiance was estimated as:

$$L_p=G_{rescale}\times D{N}_{min}+B_{rescale}-0.01\frac{\left[E_0 cos\left({\theta }_z\right)T_z+E_{down}\right]T_v}{\pi }$$

(4)

where $D{N}_{min}$ the minimum DN value for each scene, which was selected as the darkest DN with at least 1000 pixels for the entire image.

2.2.2. NDVI and Fractional Vegetation Cover Calculation

The Normalized Difference Vegetation Index (NDVI) was calculated using the atmospherically-corrected at-surface reflectance of Near IR band and Red band of TM data. Linear Spectral Mixture Analysis (LSMA) was applied to derive fractional vegetation cover. First, the Principle Component Analysis (PCA) was applied to all the spectral bands, and the two highest ranked components were selected and plotted. Second, the endmembers were selected. According to Johnson et al. [20], the potential endmembers lay at the vertices of these PCA band scatter plots. The three selected endmembers were vegetation, bare soil, and impervious surfaces. Details about the selection of end-members and estimation of the fraction were discussed in Weng et al. [21].

2.2.3. LST Computation

After converting the DN of the Thermal IR band in equation (1), the scene-specific atmospheric correction for the Thermal IR band was applied using the following equation from Coll et al. [22]:

$$L_{corrected}= \frac{L_{at-sensor}-L_{up}}{\epsilon \tau }- \frac{1-\epsilon }{\epsilon }{L}_{down}$$

(5)

where $L_{up}$ and $L_{down}$ are the upwelling radiance and downwelling radiance, respectively. $L_{at-sensor}$ is the radiance received by the satellite sensor. $\tau$ is transmittance, and $\epsilon$ is emissivity.

The transmittance, upwelling radiance, and downwelling radiance were calculated using the NASA atmospheric correction calculator [23] with a series of atmospheric profiles interpolated to particular dates as inputs for the radiative transfer model. Hence, the atmospherically-corrected radiance was converted to LST in Equation (6):

$$T= \frac{K_2}{In\left(\frac{K_1}{L_{corrected}}+1\right)}$$

(6)

where $T$ is LST in Kelvin, $L_{corrected}$ is the atmospherically corrected radiance. $K_1$ and $K_2$ are TM thermal band calibration constants, valued 607.76 $\text{W}/({\text{m}}^2\text{sr μm})$ and 1260.56 Kelvin, respectively. The LST was then converted to Celsius. The computation of emissivity values followed the procedures of [24,25].

2.2.4. Soil Moisture Computation

Soil moisture was retrieved based on the surface energy balance modeling [10,26]. First, the extension of the Priestley-Taylor parameter [10] was estimated in the TVX space; second, based on the results, the evaporative fraction was calculated; finally, the volumetric soil moisture was computed based on its relationship with evaporative fraction [26].

The detailed description of the TVX space can be found in [2]. In this case, φ is the extension of the Priestley-Taylor parameter α. In Figure 3, point E represents the pixel with the highest temperature and the lowest φ value, and point D is the lowest temperature and the highest φ value. Point A has the highest φ value, because it is on the wet edge. In the similar triangles ABC and ADE, the condition AC/AE = BC/DE can be obtained. Therefore, ${\text{φ}}_{\text{i}}$ can be written as:

$${\phi }_i= \frac{T_{max}-T_i}{T_{max}-T_{min}}\left({\phi }_{max}-{\phi }_{min}\right)+{\phi }_{min}$$

(7)

where ${\text{T}}_{\text{max}}$ is the maximum scaled LST, ${\text{T}}_{\text{min}}$ the minimum. ${\text{φ}}_{\text{max}}$ is the maximum $\text{φ}$ value, and ${\text{φ}}_{\text{min}}$ the minimum. ${\text{T}}_{\text{i}}$ and ${\text{φ}}_{\text{i}}$ are the LST and $\text{φ}$ values of pixel i. LST was transformed to scaled LST before the TVX space was formed. The equation for transforming LST to scaled LST was written as:

$$T^{*}=\frac{T-T_{min}}{T_{max}-T_{min}}$$

(8)

Then the evaporative fraction can be estimated by the following equation:

$$EvaporativeFraction=\phi \times \frac{\text{∆}}{\text{∆}+\gamma }$$

(9)

where $\text{φ}$ is the extension of the Priestley-Taylor parameter α, and$\text{ ∆}/(\text{∆}+\text{γ})$ is the air temperature control parameter [10]. In this study, the maximum value of $\text{φ}$ (under the wet surface condition) was assumed to be 1.26, and the minimum was 0 (under the dry surface condition). $\text{∆}/(\text{∆}+\text{γ})$ was estimated using a linear function with air temperature [27].

The method to calculate soil moisture from evaporation fraction was adapted from Lee and Pielke [28]:

$$\theta =\frac{{\theta }_{fc}\times arccos(1-\sqrt{4\times EF})}{\pi }$$

(10)

where θ is volumetric soil moisture, EF is evaporative fraction, and ${\theta }_{fc}$ is volumetric soil moisture at field capacity. Silt loam and clay loam were commonly found in our study area. Based on the soil-water characteristics of different soil types in Lee and Pielke [28], ${\theta }_{fc}$ was assumed to be 0.3.

Figure 3. Method to interpolate ${\text{φ}}_{\text{i}}$ for each pixel. ${\text{φ}}_{\text{i}}$ is the $\phi$ value for each random pixel in the TVX space. In the similar triangles ABC and ADE, AC/AE = BC/DE. ${\text{φ}}_{\text{i}}$ can be calculated using the values of $T_{min}, T_{max}, T_i, {\phi }_{max}$ and ${\phi }_{min}$.

Figure 3. Method to interpolate ${\text{φ}}_{\text{i}}$ for each pixel. ${\text{φ}}_{\text{i}}$ is the $\phi$ value for each random pixel in the TVX space. In the similar triangles ABC and ADE, AC/AE = BC/DE. ${\text{φ}}_{\text{i}}$ can be calculated using the values of $T_{min}, T_{max}, T_i, {\phi }_{max}$ and ${\phi }_{min}$.

2.2.5. NDWI Calculation

The NDWI was developed to depict open water present in remotely sensed images by using Near IR and Visible Green light to enhance water and to eliminate the presence of soil and terrestrial vegetation features [29]. Instead of focusing on the open water features, Gao [15] suggested that NDWI can also be an indicator of vegetation liquid water as a supplement for NDVI. Gao’s NDWI was defined in Equation (11), where ρ represents the reflectance. The lab results by Gao [15] showed that as leaf layer increased, NDWI increased, suggesting that NDWI was sensitive to the total amount of liquid water in stacked leaves. In the range from 0 to 0.15 of NDWI, NDVI was saturated while NDWI remained sensitive to liquid water in green vegetation. The spatial distribution of NDVI and NDWI were similar.

According to Gao [15], the calculation of NDWI was based on the Landsat TM data. Since the Landsat TM sensor did not cover the wavelength of 1.24 μm, it was replaced by the Landsat band 5 (1.55–1.75 µm). Thus $NDWI$ was calculated using Equation (12):

$$\text{NDWI}=\frac{\text{ρ}(0.86\text{μm})-\text{ ρ}(1.24\text{μm})}{\text{ρ}(0.86\text{μm})+\text{ ρ}(1.24\text{μm})})$$

(11)

$$NDWI= \frac{\rho 4-\rho 5}{\rho 4+\rho 5}$$

(12)

where $\rho$ 4 and $\rho$ 5 represent the reflectance for band 4 and band5.

3. Results

3.1. Impact of Urbanization-Associated LULC Changes in Three Selected Areas

The biophysical parameters in 2001 and 2006, and the percentage of changes from 2001 to 2006, using Area 1 as an example, are shown in Figure 4. As shown in Figure 4a,b, land cover types of Area 1 were changed from cultivated crops, deciduous forest, and open water to medium-density and low-density developed land in the area from 2001 to 2006. Figure 4c shows the new housing development in the southeast side of the city in 2001. Some houses and units had been built, and impervious surface was sparse and scattered. In Figure 4d, the impervious surface pixels were connected, which indicated that impervious surfaces had spread out and covered the majority of the area in 2006. The percentage of impervious surface was largely under 80%, with only a few pixels having a value higher than 80%. According to a high-resolution aerial photo, those pixels were located at either the turning of the street, or the intersection of streets and driveways. The spatial pattern of imperviousness showed the characteristics of a typical suburban residential neighborhood. Compared to the urban residential area, it usually had less dense single-family detached homes with trees and lawns, and curved roads.

Figure 4. Land cover, the percentage of impervious surface, scaled vegetation fraction, scaled Land Surface Temperature (LST), soil moisture, and Normalized Difference Water Index (NDWI) [15] in 2001 and 2006, and the percentage of changes of each parameter from 2001 to 2006. Land cover change and percent of imperviousness were acquired from National Land Cover Database (NLCD), while other parameters were generated from this study.

Figure 4. Land cover, the percentage of impervious surface, scaled vegetation fraction, scaled Land Surface Temperature (LST), soil moisture, and Normalized Difference Water Index (NDWI) [15] in 2001 and 2006, and the percentage of changes of each parameter from 2001 to 2006. Land cover change and percent of imperviousness were acquired from National Land Cover Database (NLCD), while other parameters were generated from this study.

Generally, the distribution of LST was consistent with the changes in land cover and the impervious surface fractions. In 2001 (Figure 4i), the cultivated land in the middle of Area 1 yielded relatively high temperatures because of the lower vegetation cover (Figure 4f). Forests and highly vegetated cultivated land in the northeast and the lake in the northwest yielded the lowest LST. Developed land at the southwest yielded the highest LST. In 2006 (Figure 4j), LSTs of forest patches and agricultural land remained low, while LSTs of the developed area became high.

Soil moisture decreased from 2001 to 2006. The result ranged from 18.2% to 22.0% with a mean of 20.1% in 2001 and from 17.6% to 21.5% with a mean of 19.6% in 2006 in Area 1. The areas with higher LST yielded lower soil moisture. Table 1 shows that the standard deviation increased from 0.048 to 0.058 for scaled LST, and from 0.006 to 0.008 for soil moisture. The increase of standard deviations indicated that urbanization contributed to the variations of LSTs, suggesting the possible increase in LST heterogeneity.

The distributions of NDWI shared a similar pattern to fractional vegetation cover. Dense vegetated areas yielded high NDWI values, while sparsely vegetated areas had low values. However, the mean NDWI increased from 0.081 to 0.089 from 2001 to 2006, which was contradictory to the results of vegetation cover and soil moisture. This was because the vegetation cover was not completely removed in the urbanization process: grass and small trees were planted around the houses and along the roads. As a measure of vegetation liquid water, NDWI largely depended on vegetation abundance and type, which determined the capability of vegetation for holding water, especially in leaves. Therefore, NDWI value of cultivated land was not necessarily higher than residential land. Relatively low value of NDWI in 2001 was consistent with the low vegetation fraction value of the same year.

Moreover, as Gao [15] pointed out, NDWI cannot remove the background soil reflectance completely. Wet soil with green vegetation yielded higher NDWI than dry soil with green vegetation in almost the entire range of vegetation fractions. An aerial photo taken on 26 July 2006 captured a few ponds in Area 1. Since these ponds can increase the soil moisture around the ponds, consequently higher NDWI was detected. Open water was removed from the study by Gao [15]; however, in this study, the existence of open water might be one possible reason for the relatively high NDWI values. Open water such as ponds in residential areas was included due to its potential use to track land use and land cover change and the associated impact on surface moisture conditions.

Similar observations were made for Area 2 (forest to commercial land) and Area 3 (open space to commercial land). Table 1 summarized the mean and standard deviation of each parameter, and the difference between the mean value of 2001 and 2006 for each parameter. Area 2 and Area 3 were part of the commercial land at the intersection of Highway 465 and West 86th Street. The shopping mall, restaurants, the hospital, the preschool, and large parking lots represented the major land use in this area. Area 2 and Area 3 showed sharper changes of all the biophysical parameters compared to Area 1. The mean value of vegetation fractions decreased more than 50%; the mean value of scaled LST increased more than 50%; the mean value of percent impervious surface increased around 700%; the mean value of soil moisture decreased more than 15%; and the mean value of NDWI decreased more than 90%. In 2006, the mean value of NDWI in Area 3 was negative. As a measure of vegetation liquid water and background soil reflectance, the negative NDWI value indicated that almost all of the land surfaces were modified to impervious surfaces, and that nearly no vegetation, wet soil, or pond remained.

Table 1. The mean and standard deviation of scaled vegetation fraction (Scaled Fr), scaed LST, soil moisture, NDWI, and percentage of impervious surface in 2001 and 2006.

	Area1 Cultivated to Residential (sample size 3264 pixels)					Area2 Forest to Commercial (sample size 192 pixels)					Area3 Open area to Commercial (sample size 225 pixels)
	2001		2006			2001		2006			2001		2006
	mean	std	mean	std	difference	mean	std	mean	std	difference	mean	std	mean	std	difference
Scaled Fr	0.508	0.111	0.465	0.105	− 0.043	0.715	0.093	0.351	0.114	− 0.364	0.549	0.083	0.291	0.052	− 0.258
Scaled LST	0.370	0.048	0.418	0.058	0.101	0.316	0.041	0.581	0.054	0.265	0.404	0.057	0.632	0.044	0.228
Soil moisture	0.201	0.006	0.196	0.008	− 0.005	0.209	0.005	0.173	0.008	− 0.036	0.197	0.008	0.166	0.006	− 0.031
NDWI	0.081	0.176	0.089	0.145	0.008	0.397	0.102	0.040	0.144	− 0.357	0.139	0.118	-0.038	0.075	− 0.177
Imperviousness	0.090	0.181	0.216	0.240	0.126	0.067	0.164	0.561	0.373	0.494	0.081	0.138	0.633	0.337	0.552

Note: The difference column was calculated by subtracting mean value of 2001 from that of 2006 for each parameter and each land cover conversion type.

Table 1. The mean and standard deviation of scaled vegetation fraction (Scaled Fr), scaed LST, soil moisture, NDWI, and percentage of impervious surface in 2001 and 2006.

	Area1 Cultivated to Residential (sample size 3264 pixels)					Area2 Forest to Commercial (sample size 192 pixels)					Area3 Open area to Commercial (sample size 225 pixels)
	2001		2006			2001		2006			2001		2006
	mean	std	mean	std	difference	mean	std	mean	std	difference	mean	std	mean	std	difference
Scaled Fr	0.508	0.111	0.465	0.105	− 0.043	0.715	0.093	0.351	0.114	− 0.364	0.549	0.083	0.291	0.052	− 0.258
Scaled LST	0.370	0.048	0.418	0.058	0.101	0.316	0.041	0.581	0.054	0.265	0.404	0.057	0.632	0.044	0.228
Soil moisture	0.201	0.006	0.196	0.008	− 0.005	0.209	0.005	0.173	0.008	− 0.036	0.197	0.008	0.166	0.006	− 0.031
NDWI	0.081	0.176	0.089	0.145	0.008	0.397	0.102	0.040	0.144	− 0.357	0.139	0.118	-0.038	0.075	− 0.177
Imperviousness	0.090	0.181	0.216	0.240	0.126	0.067	0.164	0.561	0.373	0.494	0.081	0.138	0.633	0.337	0.552

Note: The difference column was calculated by subtracting mean value of 2001 from that of 2006 for each parameter and each land cover conversion type.

3.2. Land Cover Types and Their Surface Characteristics

Each land cover type has a unique signature of a combination of biophysical parameters. To compare the differences of biophysical parameters in each land cover type, the alteration of land cover from green space to impervious surface resulted in a significant increase of LST and a moderate decrease of soil moisture (Figure 5). Soil moisture was not significantly related to urban land cover change, but small changes in soil moisture largely affect evapotranspiration and thus the surface energy budget [11]. This was according to the results of the three selected areas. A test covering the whole study area is also desirable.

Figure 5. The characteristics measured by scaled vegetation fraction (Scaled Fr), scaled LST, soil moisture, NDWI and percentage of imperviousness for different land types in three sample areas on 17 June 2001 and 1 July 2006.

Figure 5. The characteristics measured by scaled vegetation fraction (Scaled Fr), scaled LST, soil moisture, NDWI and percentage of imperviousness for different land types in three sample areas on 17 June 2001 and 1 July 2006.

Figure 6. Scatterplot of scaled LST (x axis) versus scaled fractional vegetation cover (y axis) for Landsat TM images that was acquired on 17 June 2001 and 1 July 2006. Compared to the shape of the scatter plot in 2001(upper), 2006 one (lower) became “shorter” and “wider”, which indicated the general trend of the surface condition changed to lower vegetation cover, lower moisture availability, and higher temperature. The scaled LST was transformed from LST by Equation (8) using maximum, minimum, and average LST. The scaled fractional vegetation cover was transformed by the same method.

Figure 6. Scatterplot of scaled LST (x axis) versus scaled fractional vegetation cover (y axis) for Landsat TM images that was acquired on 17 June 2001 and 1 July 2006. Compared to the shape of the scatter plot in 2001(upper), 2006 one (lower) became “shorter” and “wider”, which indicated the general trend of the surface condition changed to lower vegetation cover, lower moisture availability, and higher temperature. The scaled LST was transformed from LST by Equation (8) using maximum, minimum, and average LST. The scaled fractional vegetation cover was transformed by the same method.

The TVX space for the whole study area needs to be examined in detail to better understand the impact of urbanization on the selected areas (Figure 6). The most noticeable change in the TVX space was that the shape of the scatterplot appears “shorter” and “wider” in 2006 than in 2001. This change suggests that the dense vegetated area decreased, and temperature difference between wet edges and dry edges became larger. The pixels in the urbanization process migrated from the upper-left corner to the lower-right corner, which was consistent with our findings shown in the three selected areas in Figure 7.

Figure 7. (a) Pixel trajectories in the TVX space, (b) Temperature-soil moisture space, and (c) Temperature-NDWI space from 17 June 2001 to 1 July 2006. Cultivated to residential was represented by Area 1, Forest to commercial was represented by Area 1, and open area to commercial was represented by Area 3.

Figure 7. (a) Pixel trajectories in the TVX space, (b) Temperature-soil moisture space, and (c) Temperature-NDWI space from 17 June 2001 to 1 July 2006. Cultivated to residential was represented by Area 1, Forest to commercial was represented by Area 1, and open area to commercial was represented by Area 3.

To view individual panels in Figure 7 separately, the distance and slope of different land cover changes on certain biophysical parameters were compared. To view the three figures together, the impact of LULC changes on three biophysical parameters were compared. Figure 7a–c shows the pixel trajectories in the TVX space, temperature-soil moisture space, and temperature-NDWI space for the three selected areas, respectively. The starting and ending points of the vectors were located at the average value of selected parameters. The length of the trajectories shows the degree of alteration, and the slope of the vectors shows the change rate. Generally speaking, the pixels moved from densely vegetated, high moisture, and low temperature conditions, to sparsely vegetated, low moisture, and high temperature conditions. Forest to commercial land had the largest degree of alteration, and cultivated to residential land yielded the smallest alteration. The change rates in LST-soil moisture space were almost the same, which indicated that the change rate of soil moisture between different land cover types was constant. On the contrary, in the LST-NDWI space, the change rates of NDWI differ between the three types of land cover change. NDWI decreased faster from forest to commercial land than from open area to commercial land, while increasing slightly from cultivated to residential land.

4. Discussion

The major contribution of this study was to examine the impacts of LULC changes on the thermal variations as well as the patterns of soil moisture by utilizing two Landsat images. In addition to the investigation, the changes of surface temperature in TVX space, the LST-soil moisture space and LST-NDWI space were also explored to assess the impacts of LULC changes. The results pinpointed that the rate and path of temperature and moisture changes among land cover types were very different. Soil moisture and NDWI were both tested as measures of surface moisture in the urban areas. NDWI was proven to be a useful measure of vegetation liquid water. Compared to soil moisture, NDWI was more sensitive to the land cover changes. The impact of vegetation on LST was due mainly to its transpiration. To some degree, vegetation liquid water reflected the moisture condition of the environmental setting, and the liquid content in leaves may directly affect the amount of transpiration, and, consequently, LST variations. NDWI has its advantage in monitoring the moisture in both urban and suburban areas, because it can reflect both phenological cycles of vegetation and the background soil reflectance. Its relationship with LST needs to be tested in the future. In addition, other methods to measure the surface moisture condition of impervious surface are desirable.

Furthermore, soil moisture was found not to change significantly with the LULC changes. This was consistent with [11], which suggested that surface moisture availability was not significantly related to changes in urban land cover. Although the change of surface moisture was limited, its impact on evapotranspiration was substantial. Furthermore, evaporation fraction has been related to urban runoff [2]. Thus, evaporation fraction can be another important parameter to assess urban microclimate and hydrology. Its relationship with soil moisture and urban runoff can be utilized to further investigate the characteristics of impervious surfaces in the urban areas.

The evaporation and soil moisture were successfully captured using their relationship with LST in the TVX space. Based on the experiment with the three selected areas, the soil moisture showed reasonable results along with other parameters. Although it still needs to be validated with in situ data, this method provides a useful and simple way to derive evaporation and soil moisture directly from satellite images and apply it in urban environment studies. At this time, we do not have ground measured data of soil moisture. By comparing our estimation with a previous study on actual soil moisture [26], it is found that our results fell within the reasonable range. The methodology developed in this study can also be applied to thermal images from other sensors. As long as the land surface temperature data is available, the soil moisture can be estimated in the TVX space. In the Landsat Data Continuity Mission, the thermal data from Landsat 8 is useful for the continuity of surface temperature and moisture assessment. Landsat 8 collects image data by two Thermal IR bands, which allows atmospheric correction of the thermal data using a split-window algorithm [30,31].

It is assumed that air temperature was the same for the whole study area. The ground stations in two airports showed the records were quite similar: on June 16, 2001, the air temperature difference was 1 K; on 1 July 2006, the air temperature was the same. Therefore, the temperature records from one of the airports (Indianapolis International airport) was used because the airport and the three selected areas were close and had similar land cover types nearby. The temperature record from Eagle Creek Airpark Airport was not used because it is located close to a big water body: Eagle Creek Reservoir; therefore, the record may be affected by the amount of evaporation from the lake. A sensitivity analysis has been conducted. When there is a 5 °C difference in the air temperature, the difference between Evaporative Fraction is 0.04, and the difference between soil moisture is 0.01. Since the air temperature difference in the whole study area is less than 1 °C, its impact on soil moisture estimation can be ignored.

5. Conclusions

In this study, we assessed the impact of urbanization-associated LULC changes on surface temperature and moisture availability in three selected areas within Marion County, Indiana, USA, between 2001 and 2006. The significant achievement is that we applied the methods that developed in forest and agricultural areas to urban surface moisture estimation. The selected areas exemplified land cover changes from cultivated to residential land, from forest to commercial land, and from open area to commercial land. Fractional vegetation cover, LST, soil moisture, NDWI and the percentage of impervious surface were chosen as the parameters to monitor the urbanization associated changes. The results showed that, compared to the area with LULC change from cultivated to residential land, LULC change from forest to commercial land altered the surface temperature and moisture more intensively. For example, compared to the area changed from cultivated to residential land, the area changed from forest to commercial altered 48% more in fractional vegetation cover, 71% more in LST, and 15% more in soil moisture, compared to the initial state. According to the two images, the general change patterns of pixels in response to the LULC changes were from low temperature, dense vegetation cover, high surface moisture condition to high temperature, sparse vegetation cover, low surface moisture condition in the TVX space, T-soil moisture space, and T-NDWI space. These change patterns were found in the whole study area, as well as in the selected study areas. The study is unique because soil moisture and NDWI were both tested as measures of surface moisture in the urban areas. NDWI was proven to be a useful measure of vegetation liquid water. Compared to soil moisture, NDWI was more sensitive to the land cover changes. For example, for a change from forest to commercial land, the mean soil moisture changed 17%, while the mean NDWI changed 90% relative to the initial state. Future work can be directed by adding urban morphology into the estimation of urban surface moisture conditions.