Assessing the Effect of Urban Growth on Surface Ecological Status Using Multi-Temporal Satellite Imagery: A Multi-City Analysis

Abstract

Quantification of Surface Ecological Status (SES) changes is of great importance for understanding human exposure and adaptability to the environment. This study aims to assess the effect of urban growth on spatial and temporal changes of SES over a set of neighboring Iranian cities, Amol, Babol, Qaemshahr, and Sari, which are located in moderate and humid climate conditions. Firstly, the built-up footprint was prepared using Landsat images based on the Automatic Built-up Extraction Index (ABEI). Then, the surface biophysical characteristics were calculated. Secondly, the SES was modeled using the Remotely Sensed Ecological Index (RSEI), and the spatio-temporal changes of the SES were evaluated. The results revealed that the average RSEI for these cities increased from 0.48, 0.51, 0.53, and 0.55 in 1986 to 0.69, 0.77, 0.75, and 0.78 in 2022, respectively. The proportion of the poor ecological condition class in these cities rose from 10%, 3%, 5%, and 1% to 74%, 64%, 54%, and 41% during the 1986–2022 period. Our findings indicate that the SES of these cities significantly decreased while they experienced large physical growth. The findings and the methodical approach of the study provide a data-driven approach for monitoring SES in fast growing regions, which is required for studying the impact of climate change on society.

1. Introduction

We live on a dynamic planet in which the Earth surface is continuously changing, especially within and on the outskirts of cities in developing countries [1,2,3,4]. These transformations have a significant impact on energy balance, environmental elements, weather patterns, and climatic conditions [2,5,6,7,8,9]. Therefore, understanding the drivers behind these changes and managing or mitigating their adverse effects is of utmost importance. These drivers and sources of extensive change can be broadly categorized into two main groups: controllable and uncontrollable. It is important to note that the majority of factors influencing changes in land surfaces are of anthropogenic origin [10,11,12]. Population growth, urban expansion, and the conversion of natural landscapes into developed areas are among the most significant controllable factors originating from human activities [13].

Urban expansion is a pivotal factor driving extensive alterations in natural landscapes, land use/cover, ecosystems, Surface Ecological Status (SES), urban hydrology, amplified urban heat islands, and climate transformation [14]. One of the most significant adverse outcomes of urban expansion is the transformation of SES, which encompasses the detrimental consequences of human actions and environmental shifts on the ecosystem [15]. The SES serves as an evaluation of the health and functionality of the ecosystem [16]. Considering the rapid population growth and the escalating adverse effects it brings, further shifts in SES are anticipated [17]. Hence, it is crucial to identify and continually monitor the influence of urban expansion on SES and to model these shifts with an appropriate spatial and temporal resolution [18].

One of the methods for analyzing SES is field observation, or in-field surveying. However, it has certain drawbacks, including insufficient coverage, insufficient spatial accuracy, and inefficiency for inaccessible areas [19]. In contrast, the use of remote sensing (RS) data to monitor SES has become popular in recent years. The advantages of RS data, such as ease of access, extensive and continuous coverage at different scales, near-real-time availability, and extraction of surface biophysical characteristics, are effective tools to improve SES management [20]. Previous research has extensively employed spectral indices derived from satellite imagery to model Surface Ecological Status (SES) [17,21,22,23,24,25,26,27,28,29]. These indices encompass the Normalized Difference Built-up Index (NDBI), Normalized Differential Build-up and Bare Soil Index (NDBSI), Normalized Difference Vegetation Index (NDVI), Soil-Adjusted Vegetation Index (SAVI), Leaf Area Index (LAI), Land Surface Moisture (LSM), Normalized Difference Water Index (NDWI), Normalized Difference Soil Index (NDSI), and Land Surface Temperature (LST). Moreover, the combination of reflective and thermal data obtained from remote-sensing sensors has a high ability to model the parameters related to the Pressure–State–Response (PSR) framework [25,27,30]. Recently, many studies have shown that the combination of several remote sensing indicators can model the SES with appropriate accuracy [16,31]. Some of these composite models include the Ecological Index (EI) [32,33], Remote Sensing-based Ecological Index (RSEI) [34], Comprehensive Ecological Evaluation Index (CEEI) [16], Arid Remote Sensing Ecological Index (ARSEI) [35], Moving Window-based RSEI (MW-RSEI) [36], Remotely Sensed Urban Surface Ecological Index (RSUSEI) [37], Land Surface Ecological Status Composition Index (LSESCI) [25], and Remote Sensing Ecological Index considering Full Elements (RSEIFE) [38].

Previous studies have concentrated on assessing ecological conditions in various settings, including urban areas [27,30,39,40,41], forests [26], natural areas [42], agricultural regions [43,44], watersheds [45], wetlands [27], mining sites [25,46], etc. Many studies have focused on assessing the impact of land cover changes resulting from urban growth on the SES of cities. For example, Hu and Xu [31] analyzed the SES variations in the city of Fuzhou from 2000 to 2016. They introduced the Remotely Sensed Ecological Index (RSEI) to identify SES variations based on the integration of surface vegetation cover, impervious cover, LST, and wetness information. They applied Principle Component Analysis (PCA) based on the framework of the PSR for the integration of surface information in SES modeling. Their results showed that urban physical expansion led to the deterioration of surface ecological conditions, while urban afforestation contributed to an enhanced ecological quality within the urban core. Afterward, RSEI was used in several studies to evaluate the spatial and temporal changes of SES [17,21,34]. Xu et al. [34] analyzed the effect of population growth and land cover characteristics on regional ecological conditions in northern China. Their study results demonstrated that the expansion of urban infrastructure and population growth resulted in an increase in impervious surface cover (ISC), which had a notably adverse impact on regional surface ecological conditions. Zhang et al. [47] studied the effect of land cover change on SES in Chinese coastal cities using Landsat satellite imagery. Liu et al. [48] evaluated the SES of Xiamen and Kinmen Islands in southeastern China at different urban densities using RSEI.

Climate change, a pressing threat to human civilization, has its roots in escalating anthropogenic activities, with urban growth playing a pivotal role in this phenomenon [40]. The expansion of urban areas, particularly in regions characterized by humid climates, leads to the conversion of agricultural lands and green spaces into impermeable surfaces. This transformation has far-reaching consequences, including the exacerbation of climate change and the droughts that have become increasingly prevalent worldwide. In the context of climate change and drought, cities situated in humid climates hold greater significance compared to their counterparts in other climatic zones for two key reasons: (1) The escalating population and physical expansion of these humid climate cities exert more detrimental impacts on surface biophysical characteristics, as supported by previous research [30,34,49,50,51,52,53,54]. (2) Cities in humid climates are inherently more susceptible to shifts in climatic conditions and drought events when compared to cities in other climatic zones [34]. Consequently, the adverse consequences of urban growth in these cities on global climatic and environmental conditions are more pronounced.

The surface biophysical properties in the northern cities of Iran and the Caspian Sea basin characterized by moderate and humid climatic conditions have undergone significant changes primarily driven by population growth and the physical expansion of urban areas. These transformations are particularly notable due to the moderate and humid climatic conditions of the region, as well as the strategic positioning of green spaces and agricultural lands within and around these cities. These cities are in close proximity to each other, and the rural areas with substantial populations have witnessed expansion, largely facilitated by favorable weather conditions and fertile soil. Collectively, these factors have endowed this region with paramount importance in terms of surface characteristics and their impact on urban growth. The substitution of agricultural lands and green spaces with built-up areas induces alterations in surface biophysical properties. These changes involve a reduction in surface greenness and wetness, along with an increase in LST and surface imperviousness [17,34,55,56]. Assessing and quantifying the variations in each of these biophysical characteristics resulting from the physical expansion of cities have been the focus of numerous studies in the past [57]. Therefore, the assessment and quantification of changes in SES, a product of the conditions of four surface biophysical characteristics (wetness, greenness, LST, and imperviousness) stemming from urban growth, are of paramount significance, as underscored in this study. Therefore, this study aims to evaluate the influence of urban growth on the spatial and temporal changes in SES within northern cities of Iran characterized by moderate and humid climatic conditions. In this study, for the first time, the effect of urban growth on SES is evaluated for cities located in moderate and humid climates in Iran.

2. Study Area

The climatic and geographical attributes of cities situated in the northern plains of Iran closely resemble those of many cities worldwide, particularly those in close proximity to the sea. Consequently, this region has been chosen as a test bed, representing similar regions across the globe, to illuminate and substantiate the detrimental influence of population growth and the physical expansion of cities on SES. The study area encompasses four cities: Amol, Babol, Qaemshahr, and Sari, along with their surrounding suburbs, all located within the province of Mazandaran, Iran (as depicted in Figure 1). These cities share similar climatic and weather conditions, characterized by a moderate and humid climate, owing to their proximity and comparable geographical and topographical settings. Notably, they represent the largest cities in Iran within this specific climate category. The respective average elevations above sea level for Amol, Babol, Qaemshahr, and Sari are approximately 76 m, −2 m (below sea level), 51.2 m, and 54 m. In terms of land area, these cities span 38.5 km², 34.5 km², 33.1 km², and 44.3 km², respectively. Amol is one of the oldest cities in Iran and serves as the industrial and export hub of the Mazandaran Province. It is also known as the rice capital of Iran. According to the latest population census, the city is home to 237,528 people. The population density in Amol is 6169 people per km². Babol, on the other hand, is one of the key cities in the north of the country, known for its significance in the fields of medicine, agriculture, education, politics, culture, and commerce. It is also the largest county in the Mazandaran Province in terms of area. The main sources of income for the people of this city are agriculture, particularly rice cultivation, and the buying and selling of consumer goods. The people of Babol actively participate in government and local employment opportunities. Additionally, urban construction and development have been significant, largely driven by investment. Over the past three decades, Babol has witnessed substantial expansion in its built environment and, among the cities in the Mazandaran Province, it has experienced the most significant physical development after Sari. Furthermore, Qaemshahr is an industrial and export center in Mazandaran and is a significant city in terms of transportation, agriculture, tourism, mining, and industry in Iran. According to the latest population census, the city has a population of 204,953 people. The population density in this city is 6192 people per km². Sari is the most populous and largest city in the Mazandaran Province. Due to its strategic location on the transportation route from eastern to central Iran, Sari holds particular significance. The region’s economy is primarily based on agriculture, animal husbandry, handicrafts, and fishing. According to the latest population census, the city of Sari has a population of 309,820 people, with a population density of 6993 people per km².

3. Data and Methods

3.1. Data Collection

Numerous studies have used remote sensing data to model and monitor SES [58]. Similarly, this study used a set of remote sensing and climatic data. Remote sensing data and products include Landsat and Moderate Resolution Imaging Spectroradiometer (MODIS) images. Landsat 5 and 8 images were used for extracting built-up lands and surface biophysical characteristics including the Automatic Built-up Extraction Index (ABEI), LST, Normalized Difference Vegetation Index (NDVI), and wetness. In addition, the daily atmospheric profile containing product of water vapor (MOD07) at the Landsat satellite overpass time was used to estimate the LST from Landsat thermal data based on the Single-Channel (SC) algorithm.

Climatic data were employed to calculate atmospheric water vapor content. Considering the influence of climatic conditions and seasonal fluctuations on SES, the criteria for satellite image selection included (1) clear skies during satellite acquisition and (2) a lack of rainfall in the two days leading up to the satellite overpass. The specifics of the data utilized are outlined in

Table 1. Details of the research data.

Data Type	Satellite (Sensor)	Acquisition Date	Spatial Resolution of Thermal Bands (m)	Spatial Resolution of Reflective Bands (m)
Satellite Image	Landsat 5 (TM)	1986	120 (resampled to 30 m)	30
		1993
		2001
		2008
	Landsat 8 (OLI/TIRS)	2015	100 (resampled to 30 m)
		2022
	Terra (MODIS)	2001–2022	-	5000
Air Temperature	Meteorological stations	1986–2022	-	-
Relative Humidity			-	-

.

3.2. Methods

The conceptual model depicted in Figure 2 was employed to model and assess changes in the ecological conditions of the cities under investigation. Initially, Landsat images underwent pre-processing, from which maps of surface biophysical properties such as wetness, heat, greenness, and imperviousness were generated using spectral indices and various methods at different time points. Subsequently, the trends in surface biophysical property changes within the studied cities were evaluated and compared. In the subsequent step, by amalgamating information regarding surface biophysical characteristics via Principal Component Analysis (PCA), the SES was modeled utilizing the RSEI index. Finally, the alterations in the SES of the cities under investigation were computed and compared for analysis.

3.2.1. Extraction of Built-Up Areas

The ABEI index was employed for mapping built-up lands [59]. This index utilizes a combination of six reflectance bands (Equation (1)) for Landsat 5 and a combination of seven reflectance bands (Equation (2)) for Landsat 8.

$$\mathrm{L}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{s}\mathrm{a}\mathrm{t} 5=0.825\times \mathrm{B}\mathrm{l}\mathrm{u}\mathrm{e}-0.086\times \mathrm{G}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{n}-0.441\times \mathrm{R}\mathrm{e}\mathrm{d}+0.052\times \mathrm{N}\mathrm{I}\mathrm{R}-0.198\times \mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1+0.278\times \mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1$$

(1)

$$\mathrm{L}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{s}\mathrm{a}\mathrm{t} 8=0.312\times \mathrm{U}\mathrm{l}\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{B}\mathrm{l}\mathrm{u}\mathrm{e}+0.513\times \mathrm{B}\mathrm{l}\mathrm{u}\mathrm{e}-0.086\times \mathrm{G}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{n}-0.441\times \mathrm{R}\mathrm{e}\mathrm{d}+0.052\times \mathrm{N}\mathrm{I}\mathrm{R}-0.198\times \mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1+0.278\times \mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1$$

(2)

Based on collected samples derived from Google Earth images and color composites of Landsat imagery for both the built-up and non-built-up categories, the accuracy of the generated built-up maps was assessed using producer′s, user′s, and overall accuracy metrics. In cases where Google Earth images were unavailable for the study area in 1986, 1993, and 2001, only the color composites of Landsat images were utilized. The test samples for each year were strategically chosen to represent various spatial locations within the study area, encompassing cities, villages, agricultural zones, barren land, green areas, and water bodies. A minimum of 500 test samples were collected for both the built-up and non-built-up categories each year. Subsequently, changes in the built-up class area across different years were calculated and compared.

3.2.2. Spectral Indicators

The spatial and temporal pattern of physical and biophysical features on the Earth’s surface plays a significant role in the SES. To explore the relationship between changes in SES and the spatio-temporal patterns of surface biophysical characteristics, four spectral indicators were utilized, including greenness (NDVI), moisture (wetness), impervious surface (ABEI), and heat (LST), as outlined in

Table 2. Indicators and spectral methods for modeling surface properties.

Surface Biophysical Characteristics	Parameter	Equation	References
Greenness	NDVI	$\frac{\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{R}\mathrm{e}\mathrm{d}}{\mathrm{N}\mathrm{I}\mathrm{R}+\mathrm{R}\mathrm{e}\mathrm{d}}$	[61]
Impervious surface	ABEI	Equations (1) and (2)	[59]
Moisture	Wetness	$\mathrm{L}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{s}\mathrm{a}\mathrm{t} 5=0.14\times \mathrm{B}\mathrm{l}\mathrm{u}\mathrm{e}+0.17\times \mathrm{G}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{n}+0.32\times \mathrm{R}\mathrm{e}\mathrm{d}+0.33\times \mathrm{N}\mathrm{I}\mathrm{R}-0.62\times \mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1-0.41\times \mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}2$	[62]
		$\mathrm{L}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{s}\mathrm{a}\mathrm{t} 8=0.15\times \mathrm{B}\mathrm{l}\mathrm{u}\mathrm{e}+0.19\times \mathrm{G}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{n}+0.32\times \mathrm{R}\mathrm{e}\mathrm{d}+0.34\times \mathrm{N}\mathrm{I}\mathrm{R}-0.71\times \mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1-0.45\times \mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}2$
Heat	LST	$\mathrm{L}\mathrm{S}\mathrm{T}=\mathsf{\gamma }\left(\frac{1}{\mathrm{L}\mathrm{S}\mathrm{E}}\left({\mathsf{\psi }}_1{\mathrm{L}}_{\mathrm{s}\mathrm{e}\mathrm{n}}+{\mathsf{\psi }}_2\right)+{\mathsf{\psi }}_3\right)+\mathsf{\delta }$	[63]

. NDVI serves as an indicator of vegetation density, with its values ranging from −1 (indicating water) to +1 (indicating dense vegetation). The NDVI thresholding method was employed to estimate the land surface emissivity coefficient for LST estimation (refer to [60] for further details). In this study, the wetness component derived from Tasseled Cap Transformation (TCT) was applied to surface moisture estimation. Additionally, the ABEI index was utilized to estimate the extent of impervious surfaces.

The values of surface biophysical characteristics, such as NDVI, wetness, ABEI, and LST, for various dates were normalized within the range of 0 to 1. Subsequently, the mean and standard deviation of each standardized surface biophysical characteristic were calculated both for the overall study area and for each of the cities. These values were then compared. Finally, the trends in surface biophysical characteristics for different cities during the period from 1986 to 2022 were investigated and compared.

3.2.3. SES Modeling

In this study, we employed the PSR-framework-based RSEI to model SES. This framework comprises three key components: environmental states, climate responses, and anthropogenic pressures. The level of human activity within a region exhibits a robust correlation with impervious cover. Consequently, we utilized the imperviousness index to account for pressure intensity. The greenness index was applied to assess environmental states, serving as an indicator of vegetation cover and quality. Additionally, heat and moisture indices were employed to address climate responses. It is worth noting that pressure intensity on the environment resulting from human activities can also exhibit a strong correlation with heat.

In this study, the NDVI component serves as the greenness index, while the wetness component, derived through TCT, represents moisture levels. LST is utilized as an indicator of heat in response to climatic factors. ABEI, on the other hand, signifies impervious surfaces, reflecting the environmental pressures resulting from human activities. The RSEI was formulated exclusively based on remote sensing data, combining greenness, wetness, imperviousness, and heat indices through PCA. It has demonstrated a robust correlation with SES [17,31,34]. Consequently, the term “RSEI” can be expressed as Equation (3):

$$\mathrm{R}\mathrm{S}\mathrm{E}\mathrm{I}=\mathrm{P}\mathrm{C}1\left(\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I},\mathrm{W}\mathrm{e}\mathrm{t}\mathrm{n}\mathrm{e}\mathrm{s}\mathrm{s},\mathrm{L}\mathrm{S}\mathrm{T},\mathrm{A}\mathrm{B}\mathrm{E}\mathrm{I}\right)$$

(3)

To facilitate the study and comparison of spatio-temporal variations in SES derived from spectral indices, the RSEI values were normalized within the range of 0 to 1 for each image. Here, values of 0 and 1 correspond to the best and worst SES conditions, respectively. In the subsequent phase, the SES of the study area was categorized into five distinct groups: excellent (0–0.2), good (0.2–0.4), moderate (0.4–0.6), fair (0.6–0.8), and poor (0.8–1) (as outlined in

Table 3. Method of classification of SES [16].

Level	Class	Distance to Each Class	Ecological Quality Conditions
1	excellent	0–0.2	Ecosystem functions are healthy. Vegetation is very dense and no human disturbances occur.
2	good	0.2–0.4	Ecosystem functions are healthy. Vegetation is high and human disturbances are very rare.
3	moderate	0.4–0.6	Ecosystem performance is not affected. The vegetation level has reached the average level (about 50%). Human disturbances often occur, and the effect of an urban heat island is not obvious.
4	fair	0.6–0.8	Ecosystem functions have been significantly lost but are recoverable. Vegetation is poor. Human activities and the effects of the heat island are also evident in these areas.
5	poor	0.8–1	Ecosystem and vegetation functions have been completely destroyed and are very difficult to recover. Intensification of human activities, drought, and heat islands is also evident in these areas.

).

To quantify spatial and temporal changes in the SES of the study area, the following steps were taken:

• Mean values and standard deviations of RSEI for both the entire study area and individual cities were calculated and compared over the period from 1986 to 2022.
• The areas corresponding to different SES classes were computed and compared for both the overall study area and each of the cities throughout the period spanning from 1986 to 2022.

4. Results

The overall accuracy of land cover maps, generated based on ABEI for the years 1986, 1993, 2001, 2008, 2015, and 2022, stood at 92.0%, 91.1%, 93.3%, 90.3%, 92.2%, and 90.1%, respectively. Furthermore, the user′s (producer′s) accuracy for the built-up land class during these years was 91.8% (92.5%), 90.7% (91.9%), 94.2% (93.2%), 91.5% (90.5%), 92.5% (91.7%), and 90.4% (89.9%), respectively (as indicated in

Table 4. Confusion matrix and producer′s, user′s, and overall accuracy (%) of selected classes within the 1986, 1993, 2001, 2008, 2015, and 2022 land-cover maps.

1986		References Data
		Built-up lands	Non-built-up lands	Total	Producer′s accuracy (%)
Classified data	Built-up lands	620	50	670	92.5
	Non-built-up lands	55	590	645	91.4
	Total	675	640	1315
	User′s accuracy (%)	91.8	92.1		92.0
1993		References data
		Built-up lands	Non-built-up lands	Total	Producer′s accuracy (%)
Classified data	Built-up lands	582	51	633	91.9
	Non-built-up lands	59	547	606	90.2
	Total	641	598	1239
	User′s accuracy (%)	90.7	91.4		91.1
2001		References data
		Built-up lands	Non-built-up lands	Total	Producer′s accuracy (%)
Classified data	Built-up lands	626	45	671	93.2
	Non-built-up lands	38	540	578	93.4
	Total	664	585	1249
	User′s accuracy (%)	94.2	92.3		93.3
2008		References data
		Built-up lands	Non-built-up lands	Total	Producer′s accuracy (%)
Classified data	Built-up lands	674	74	748	90.1
	Non-built-up lands	62	599	661	90.6
	Total	736	673	1409
	User′s accuracy (%)	91.5	89.0		90.3
2015		References data
		Built-up lands	Non-built-up lands	Total	Producer′s accuracy (%)
Classified data	Built-up lands	720	65	785	91.7
	Non-built-up lands	58	739	797	92.7
	Total	778	804	1582
	User′s accuracy (%)	92.5	91.9		92.2
2022		References data
		Built-up lands	Non-built-up lands	Total	Producer′s accuracy (%)
Classified data	Built-up lands	701	78	779	89.9
	Non-built-up lands	74	689	763	90.3
	Total	775	767	1542
	User′s accuracy (%)	90.4	89.8		90.1

).

The spatial distribution of the built-up class has undergone significant changes throughout the study period, as illustrated in Figure 3. The physical extents of the studied cities have expanded in various geographical directions.

The area of the built-up class in the study area has increased from 80 km² in 1986 to 219 km² in 2022 (

Table 5. Area of built-up lands in the different cities and the study area for the period 1986–2022 (km²).

Area/Year	1986	1993	2001	2008	2015	2022
Amol	16.02	18.48	21.88	24.58	26.81	30.61
Babol	14.31	17.25	20.91	22.33	23.77	26.60
Sari	13.61	16.53	19.38	21.00	22.44	24.96
Qaemshahr	13.61	18.71	24.42	28.06	31.77	33.61
Whole study area	80.83	101.12	128.99	149.24	180.42	219.61

). Over the period from 1986 to 2022, the built-up class area in the study area and each of the cities, including Amol, Babol, Qaemshahr, and Sari, has increased by 171%, 91%, 85%, 83%, and 146%, respectively. In these cities, the area of the built-up class has increased by 14 km², 12 km², 11 km², and 20 km², respectively. The increasing area of built-up lands has reduced the area of agricultural lands and gardens in the study area. Among the studied cities, Sari and Qaemshahr had the highest and lowest increase in built-up areas, respectively. In 2022, the most and least built-up land areas were associated with Sari and Qaemshahr, respectively.

The maps of the standardized biophysical properties in different years show different spatial distributions (Figure 4). An examination of these surface characteristic maps reveals that in urban regions, LST and ABEI are higher compared to non-urban areas. Conversely, NDVI and wetness in urban regions are lower than in non-urban areas. In simpler terms, as urban areas expand, the regions within the urban environment characterized by high temperatures, impervious surfaces, and lower levels of greenness and wetness are expanding as well. Furthermore, due to the prevalence of bare lands surrounding the cities of Sari and Qaemshahr, the LST in the suburbs of these cities is higher than that of Babol and Amol. The mean ABEI for the study area has increased from 0.39 in 1986 to 0.49 in 2022, primarily due to the expansion of built-up areas. The standard deviation and mean of standardized values for NDVI, LST, and wetness exhibit variations corresponding to changes in land cover across different years.

The mean values of standardized biophysical properties for each city in the 1986–2022 period show an increasing and decreasing trend (Figure 5). Between 1986 and 2022, the mean standardized LST and ABEI increased for all four cities. In contrast, the mean NDVI and wetness for these cities have experienced a significant decrease during this timeframe. Specifically, the cities of Sari and Amol saw the highest and lowest increases in ABEI, respectively. The mean standardized LST for Amol, Babol, Qaemshahr, and Sari from 1986 to 2022 increased by 112%, 90%, 57%, and 35%, respectively. On the other hand, the mean standardized NDVI for these cities decreased from 0.55, 0.54, 0.57, and 0.55 in 1986 to 0.32, 0.34, 0.36, and 0.39 in 2022. The mean standardized wetness in the urban areas of Amol, Babol, Qaemshahr, and Sari during this period decreased by 47%, 37%, 36%, and 29%, respectively.

The spatial distribution of RSEI values, which represent SES, has undergone significant changes during the study period, as illustrated in Figure 6. A value approaching 1 indicates worse ecological conditions, while a value less than 1, closer to 0, signifies better ecological conditions in the area. Between 1986 and 2022, areas with RSEI values approaching 1 (indicating poor ecological conditions) have expanded within the urban environment. Agricultural areas with high wetness exhibit the lowest RSEI values in the study area, reflecting excellent ecological conditions. Similarly, orchard lands characterized by abundant greenery also boast low RSEI values, although these values are higher than those found in agriculturally wet areas. In contrast, urban areas with a high proportion of built-up cover demonstrate the highest RSEI values in the study area. This is primarily due to the limited presence of greenery and wetness, an abundance of impervious surfaces, and susceptibility to drought and heat. Given the substantial differences in surface biophysical characteristics between urban built-up areas and the agricultural lands and green spaces surrounding the cities, the ecological conditions within the urban environment starkly contrast with those in the suburbs.

The mean (standard deviation) of RSEI for the study area in 1986, 1993, 2001, 2008, 2015, and 2022 was 0.44 (0.17), 0.42 (0.19), 0.45 (0.19), 0.47 (0.18), 0.47 (0.17), and 0.44 (0.17), respectively, as shown in Figure 7. The mean RSEI values for Amol, Babol, Qaemshahr, and Sari increased from 0.48, 0.51, 0.53, and 0.55 in 1986 to 0.69, 0.77, 0.75, and 0.78 in 2022, respectively. Among these cities, Babol experienced the highest increase in mean RSEI, while Amol had the lowest increase.

The spatial distribution of SES varied heterogeneously during the 1986–2022 period, as depicted in Figure 8. Visual examination reveals a gradual shift of areas with poor and fair ecological conditions from the suburbs to urban areas. In 1986, only a negligible portion of urban areas exhibited poor ecological conditions, but over the subsequent years, this percentage experienced a substantial increase. Regions with excellent ecological conditions were predominantly found in lands characterized by very high wetness or greenness, such as forests, wetlands, and extremely wet agricultural areas. Alterations in the wetlands’ area in the study region in different years emerged as one of the primary factors contributing to spatial changes in areas with excellent ecological conditions throughout the study period.

Numerous regions within the study area exhibited favorable ecological conditions in various years, comprising agricultural lands and green spaces. However, due to the transformation of these green areas and agricultural lands surrounding cities into urbanized areas, the ecological conditions in these regions have shifted from good and moderate to fair and poor. Over the course of the study period, the extent of different ecological condition classes in the region has experienced fluctuations.

The area covered by the poor ecological conditions class in 1986, 1993, 2001, 2008, 2015, and 2022 was 39 km², 26 km², 64 km², 99 km², 72 km², and 97 km², respectively. Across all these dates, the largest area was associated with the moderate ecological conditions class. Over time, there has been a decreasing trend in the area of regions with good ecological conditions and an increasing trend in regions with poor ecological conditions. Additionally, the areas with excellent, moderate, and fair ecological conditions have shown similar and closely related trends throughout these years.

A significant portion of the SES classes has shifted from excellent, good, and moderate to fair and poor (

Table 6. Area percentage of SES classes for different cities from 1986 to 2022.

Area/Year		1986	1993	2001	2008	2015	2022
Amol	Excellent	2.763	0.744	0.072	0.011	0.012	0.009
	Good	31.876	16.463	11.605	5.979	7.880	4.880
	Moderate	37.298	4.556	25.372	20.496	27.927	27.931
	Fair	27.384	40.335	43.929	39.774	41.268	35.900
	Poor	0.676	1.869	19.019	33.766	35.909	41.300
Babol	Excellent	1.294	0.062	0.203	0.002	0.020	0.020
	Good	28.167	11.740	11.605	4.915	1.138	4.149
	Moderate	35.953	30.460	21.027	18.021	15.937	16.296
	Fair	29.308	32.648	38.674	34.987	34.805	25.174
	Poor	5.275	25.088	28.516	42.073	48.098	54.358
Sari	Excellent	0.129	0.119	0.085	0.000	0.201	0.134
	Good	10.515	7.744	11.215	2.305	1.021	0.639
	Moderate	34.248	39.115	23.892	18.434	12.546	5.390
	Fair	44.576	41.786	43.880	40.591	41.522	19.734
	Poor	10.529	11.234	20.926	38.68	44.708	74.101
Qaemshahr	Excellent	0.225	0.008	0.070	0.000	0.000	0.000
	Good	12.691	5.702	8.405	2.288	0.643	0.464
	Moderate	37.189	35.392	23.387	18.187	16.295	10.794
	Fair	46.221	48.862	43.309	48.138	41.176	24.280
	Poor	3.672	10.034	24.826	31.386	39.884	64.460

). Notably, in Sari, Qaemshahr, Babol, and Amol, the area classified as poor increased from 10%, 3%, 5%, and 1% in 1986 to 74%, 64%, 54%, and 4% in 2022, respectively.

Areas characterized by excellent ecological conditions consistently represent the smallest proportion of urban areas across all years, in contrast to the other categories. Comparatively, Babol and Amol consistently exhibit the highest and lowest proportions of areas classified as having poor ecological conditions, respectively. In 2022, the largest percentage of ecological condition classes in all four cities belongs to the poor category.

5. Discussion

Throughout the past decades, human activities such as deforestation, overgrazing, and urbanization have significantly contributed to the degradation of terrestrial ecosystems. This degradation is evident through phenomena such as soil erosion, desertification, and the loss of biodiversity. Moreover, some of these actions have led to the increased emission of human-generated greenhouse gases, which further exacerbate climate change by trapping heat in the atmosphere. The results of this study highlight that the primary drivers behind the substantial spatial and temporal changes in SES are urban physical expansion and the transformation of green spaces and agricultural lands into built-up areas. The findings reveal that the built-up land in the study area has increased by 170% over the past few decades. This significant increase has led to substantial alterations in land surface characteristics, including imperviousness, vegetation, wetness, and LST (as observed in Figure 4 and Figure 5, and detailed in Table 5). Considering the historical context and topographical features of the study area, changes in imperviousness are likely to contribute to heightened flood risks, given the region’s annual precipitation volume. Furthermore, an elevated LST, stemming from diminished vegetation, reduced wetness, and increased imperviousness, can have adverse implications for public health and the suitability of climatic conditions for certain tree species and agricultural activities [64]. Xu et al. [34] showed that the expansion of impervious surfaces leads to decreased moisture, reduced vegetation density, elevated surface heat, and increased dryness, ultimately exacerbating surface ecological deterioration. In addition, Chen et al. [42] revealed a negative correlation between the quality of surface ecological conditions and urbanization indices in their research.

The population in areas characterized by fair ecological conditions within the study area has experienced a substantial increase over the past decades, highlighting one of the most significant adverse impacts of urban expansion.

The various surface characteristics employed in the RSEI model exhibit varying degrees of importance in SES modeling, which can be quantified and compared based on the eigenvector values in PCA [25]. Within the study area, it is evident that the significance of LST in SES modeling surpasses that of other surface characteristics [34]. Firozjaei et al. [37] showed that the significance degree of LST in SES modeling is higher than those of the characteristics of vegetation, wetness, and imperviousness. Additionally, Shan et al. [32] identified positive effects on ecological condition quality through LSM and NDVI, whereas LST and NDBSI had negative effects. Notably, their results underscored the pronounced significance of NDBSI compared to other indices.

This study adopted RSEI as a suitable indicator for SES modeling that provides the following advantages: (a) comparable at different geographical scales; (b) scalable; (c) customizable to error minimization due to the weight definitions; (d) visualizable and able to separate the ecological states of different land covers [31,34]. The efficiency of this method in modeling SES has been confirmed by previous research [31,34]. Nonetheless, employing CEEI and RSUSEI for SES modeling and subsequently comparing their outcomes with those of RSEI can prove beneficial. The RSEI is one of the most widely used indices for assessing SES. This index is entirely based on remote sensing data and can quickly monitor SES aspects for a large area. However, the use of RSEI has its limitations. In summary, its application is not universally uniform, the results obtained are somewhat random, and its calculation process cannot account for all ecological elements [25,37,65,66]. Sun et al. [67] demonstrated in their study that RSEI has the least utility in highly ecological areas such as deserts, and scenarios prone to land degradation, where the results are unstable. CHENG et al. [40] showed that the contribution of PC1 (the first principal component) in the standard RSEI calculation process is not constant, typically varying from 60% to 90% in different regions. Therefore, relying solely on PC1 as an SES assessment index may lead to the loss of precise information [68,69].

Due to its potential to integrate data obtained from both reflective and thermal bands, Landsat satellite imagery can efficiently capture surface properties and simulate changes in ecological conditions and biophysical characteristics [16,34]. Data obtained from MODIS can also be valuable for assessing SES. Although MODIS data have lower spatial resolution compared to Landsat data, they can still offer higher efficiency in evaluating temporal changes. Qureshi et al. [27] demonstrated that simultaneous use of MODIS and Landsat data can be beneficial for assessing both spatial and temporal changes in SES. Furthermore, the utilization of high spatial resolution, radar images, and hyper-spectral imagery data (such as synthetic aperture radar (SAR), Hyperion, AVIRIS, Sentinel 1 and 2, Worldview-3, etc.) could offer even greater possibilities.

The results of this study shed light into how the absence of state-of-the-art nature-based solutions in urban planning can lead to reducing SES quality. Urban expansion in the vertical direction while improving greenification and urban/climate-friendly building materials should be prioritized. Forecasting SES by simulating urban growth and land cover changes can provide valuable insights for future planning. In terms of future research, it is advisable to explore the seasonal variations’ impact on SES changes, although this was constrained in our study area due to substantial cloud cover during the colder seasons. One limitation of this study was the lack of access to historical data regarding the socio-economic characteristics of the cities under investigation. Furthermore, it is important to acknowledge that economic and social conditions can influence the extent of SES changes resulting from urban expansion. Hence, it is recommended to investigate these factors in future studies. Finally, it is important to note that the SES assessment process takes into account various natural and human factors. However, the actual composition of the SES is even more intricate. When assessing a vast region, topographic elements play a pivotal role in shaping SES patterns. Factors such as slope and digital elevation model exert some influence on the SES. This aspect should be given additional attention in future SES assessments.

6. Conclusions

The primary objective of this study was to model and assess changes in ecological conditions across selected cities using satellite imagery. The results revealed that, in the studied cities, which have experienced significant expansion, areas characterized by high humidity in agriculture exhibited the most favorable SES. Similarly, garden lands with lush greenery showed good ecological conditions, albeit lower than those of agriculturally humid regions. In contrast, built-up areas had the least favorable SES among all SES classes due to their low moisture and greenery levels, coupled with high imperviousness, aridity, and heat. Consequently, the transformation of suburban agricultural and garden lands into built-up areas led to reduced surface moisture and greenery, alongside increased imperviousness, dryness, and heat, ultimately diminishing SES quality. As a result, urban expansion in recent years has had a detrimental impact on SES. Presently, roughly 36% of the population resides in areas characterized by poor ecological conditions, marking a significant increase compared to 1986. The findings indicate that the changes in ecological conditions are mainly caused by the process of changing the landscape and the surface biophysical characteristics. These results can prove quite helpful and practical to local decision makers and urban planners for improving ecological conditions while planning urban landscapes. Our study can serve as a source of inspiration for rapidly developing cities in developing countries, encouraging them to adopt a systematic approach and leverage publicly available satellite data to monitor the influence of historical developments on SES. This approach can facilitate the assessment of urban structures in alignment with innovative urban planning concepts, such as urban greening, green roofs, sponge cities, and the attainment of UN sustainable development goals.