Assessing and predicting changes of ecosystem service values in response to land use/land cover dynamics in Ibb City, Yemen: a three-decade analysis and future outlook

ABSTRACT

Assessing ecosystem services values (ESV) within land use/land cover (LULC) changes is crucial for promoting human well-being and sustainable development of regional ecosystems. Yet, the spatial relationship between LULC changes and ecosystem services is still unclear in Yemen. This study aimed to assess the impacts of LULC changes on ESV in Ibb City, over three decades (1990–2020), and to predict the ESV changes in 2050. The hybrid land use classification technique for classifying Landsat images, the CA-Markov model for LULC prediction, and the benefit transfer method (BTM) for assessment ESV were employed. Our findings revealed that there was a continuous increase in built-up areas and barren land, with a decrease in cultivated land and grassland, which are predicted to continue for the next 30 years. Consequently, the total ESV has decreased from US$ 68.5 × 10⁶ in 1990 to US$ 65.2 × 10⁶ in 2020 and is expected to further reduce to US$ 61.2 × 10⁶ by 2050, reflecting the impact of urban expansion and socio-economic activities on ESV. The study provides insights for future ESV monitoring, which will contribute to the formulation of effective land-use strategies for more sustainable ecosystem services, particularly in rapidly urbanizing data-limited regions.

KEYWORDS:

• LULC Changes
• CA-Markov
• coefficient of elasticity
• ecosystem service values (ESV)
• ecosystem services functions

1. Introduction

Ecosystem services (ES) refer to the functions and activities of natural ecosystems and their organisms that provide human well-being and sustainable development; or are the benefits that humans obtain directly or indirectly through ecosystems (such as services of provisioning, regulatory, cultural, and supporting) (Costanza et al. 1997; Liu et al. 2019; Msofe et al. 2020). Consequently, the ES not only furnishes material services to humans such as food, water, and various raw materials but also provides non-material services such as climate and gas regulation, etc., which support system for sustaining life and establishing ecological conditions for human well-being (M. Das and Das 2019; Morshed et al. 2022). Globally, the values of Ecosystem services have deteriorated dramatically over time and across different geographical areas (M. Das and Das 2019). According to (Costanza et al. 1997; 2014), the total value of global ES during the period 2007–2011 was US$145–125 trillion per year. From this perspective, land use/land cover (LULC) change is the first and direct driving force of ES, and its change process plays a crucial role in changing the functions of ES (Jiang, Fu, and Lü 2020; Msofe et al. 2020). Thus, the study of the impact of LULC change on changing ecosystem service values (ESV) is of great significance for understanding regional environmental changes, maintaining ecological balance, and promoting the coordinated development of regional economy and environment (N. Li et al. 2021). In recent years, scholars have gradually begun to strengthen research in this area and achieved a series of research results, making the research on the value of ES develops from a static to a dynamic direction (Hasan et al. 2020). However, studies have shown that in the historical period, LULC changes were mostly caused by human socioeconomic activities, thus leading to changes in the value of ecosystem service (L. Gao et al. 2019). Therefore, as the global environment deteriorates and people's awareness of environmental protection increases, scholars and research organizations globally are increasingly aware of the impact of human activities on environmental changes (B. Shrestha, Ye, and Khadka 2019).

In the literature, diverse methodologies have been developed for the assessment of terrestrial ES. These approaches can be broadly classified into four groups including (a) revealed preference techniques; (b) stated preference methodologies; (c) cost-based methodologies; and (d) the benefits transfer method (BTM) (Chopra, Khuman, and Dhyani 2022). These methods differ in the source of information and resources (monetary and time) required. Among these methods, the BTM is the most widely employed approach on a global scale for valuing ES because of its feasibility and simplicity (M. Shrestha and Acharya 2021). Most notably, these researchers categorized the global biosphere was first classified into 16 sub-ecosystems and 17 ecosystem service functions, and the value of each ES was evaluated by BTM (Chatterjee et al. 2022). Although Costanza's estimates of ESV are considered to be rough approximations, they remain the most comprehensive and widely recognized set of estimates available for measuring changes in the value of services provided by various ecosystems (Belay, Melese, and Senamaw 2022) for ecosystem service assessment research. Following this groundbreaking estimation of global ESV, scholarly interest in the valuation of ES has witnessed marked attention on an international scale among researchers (Woldeyohannes et al. 2020). For instance, in Bangladesh, (Gashaw et al. 2018; Haque et al. 2023) evaluated the impact of LULC changes on ESV; in China, Liu et al. (2019)Assessed the Impacts of urbanization-induced land-use changes on ESV; In India, (Basu et al. 2023; Chatterjee et al. 2022; M. Das and Das 2019; M. Das, Das, and Pandey 2022) analysis of ESV Based on LULC change; and in Ethiopia, (Gashaw et al. 2018; Tolessa, Senbeta, and Kidane 2017) monitored the impact of LULC change on ESV. In addition, many scholars have designed a variety of models to predict the future impact of LULC and ESV at different scales. For example, land change modeler in TerraSet (known as IDRISI) (Taiwo et al. 2023), machine learning algorithms (A. Al Kafy, Saha, et al. 2022), the Cellular Automata (CA) algorithms (Rahaman, Kafy, Faisal, et al. 2022), Multi-Layer Perceptron Neural Network and Markov Chain (MLPNN-M) (Basu et al. 2023; A. Wang et al. 2023), model future land use simulation (FLUS) model (Hou et al. 2021), artificial neural networks (ANNs) (Morshed et al. 2022), The Cellular Automata-Artificial Neural Network (CA-ANN) (A. Wang et al. 2023; A. Al Kafy et al. 2021b), RapidEye images (A. Al Kafy et al. 2022), and Markov chain and cellular automata models (J. Li et al. 2019; F. Zhang, Yushanjiang, and Jing 2019), etc. These models provide a strong foundation for evaluating the impacts on EVS of current and future LULC, as the input data requirements and simulation objects are different, the application effects have their advantages and disadvantages (A. Al Kafy et al. 2021a). The Markov model mainly has advantages at the regional macro scale, and the cellular automata (CA) model has advantages at the spatial unit scale (Y. Zhang et al. 2021). As a result, the integration of CA with a Markov model, commonly denoted as CA-Markov, within a GIS framework, has significantly minimized simulation inaccuracies (Rahaman, Kafy, Faisal, et al. 2022). This is the rationale behind employing this model in the present study for simulating perspective LULC changes and their associated ESV.

The majority of studies adopt value coefficients to approximate the ESV and create maps of their distribution through a combination of BTM with remote sensing and GIS estimation techniques (Liu et al. 2019). This approach is feasible for estimating ESV and suggests management alternatives for areas with a lack of data (Hoque et al. 2022). Furthermore, the Economics of Ecosystems and Biodiversity (TEEB) database has been established to identify value coefficients for different biomes from various regions worldwide, serving as the foundation for developing these coefficients (Msofe et al. 2020; Morshed et al. 2022). The TEEB database offers valuable insights into the valuation of ES across diverse ecosystems (Kornatowska and Sienkiewicz 2018). Despite efforts to estimate ESV for regions with data scarcity, there are few studies on the evaluation of ESV related to the dynamics of LULC change in the Middle East (Bindajam et al. 2021). In Yemen, rapid population growth and economic development are a problem of concern, where LULC types are transforming into built-up areas and influencing degradation of ecology and environment, which hinders the ability to Maintain and manage ecosystem services (Dammag, Jian, et al. 2023). One of the underlying factors for the degradation of Yemen’s ecosystems is the lack of integration of ecosystems’ values in decision-making at national and local levels, lack of acknowledgment that the ecosystem is a closed system, and that ecosystem deprecation feedback into slow economic growth (Masike 2014). Consequently, this study has entailed two areas of focus in this field. The primary focus has been given to assessing the impacts of historical LULC changes on ESV than projected changes in ESV in the future based on the predicted LULC pattern under different scenarios.

One such urban center grappling with the complexities of urbanization is Ibb city, which is located in the southwest highlands of Yemen, is the tourist capital of Yemen, and is the third-largest city of the Republic in terms of population (Al-Darwish et al. 2018). In recent years, the city has witnessed an influx of internally displaced persons (IDPs) as a result of recent conflicts in Yemen. The continuous intensification of LULC change towards agricultural land and urbanization due to population pressure, government policies, and climate change have led to ecological problems in the region, such as food production, biological control, the decrease in grassland coverage, and the degree of desertification (Dammag et al. 2023). Changes in Ibb's ecosystem and its services are closely related to maintaining regional ecological security and sustainable economic and social development (Dammag et al. 2023a). Rapid urbanization and increasing populations have led to fundamental changes in land use, which have widespread impacts on ES in the region. Therefore, understanding the impacts of LULC changes on ESV in Ibb City is essential to give public awareness and indicate policy directions about the loss/gain in ecosystem services. Previously, there is no study carried out in the study region concerning the impacts of LULC changes on ESV. Therefore, this study is considered crucial to reveal the detection impacts of LULC changes on ESV over the past and next three, which provides a basis for sustainable development and management of the ecology and environment in the study area. To achieve this, (1) to examine the temporal and spatial dynamics patterns of LUCC in the study area during the whole study period; (2) to estimate the impacts of LULC changes on ESV in monetary values in the study area over 1990–2020 periods. in addition, the study will predict the ESV in the future (2035, 2050) using the CA-Markov model to predict the trend and spatial structure of LULC classes; (3) to explore the changes in individual ESV due to the effects of dynamics in each LULC class; and (4) to assess the impact of LULC change on ESV through determining the elasticity of ESV across the entire period. This study offers a distinctive approach to estimating local-level ESV by leveraging the knowledge of local experts instead of relying solely on global databases, which may yield inaccurate or exaggerated results. In addition, the discussion offers suggestions for economic development and balancing regional differences of ESV that can inform decision-making processes and be utilized in other data-scarce regions and facing similar conditions.

2. Data and methods

2.1. Study area

The study area, Ibb city, is located in the southwestern highlands of central Yemen. It is the tourist capital of Yemen and is about 193 km away from Sana’a (the capital of Yemen) (Dammag et al. 2023). It is located between 13°54′ and 14° 73′ N latitude and 43° 08′ and 45° 21′ E longitude, as depicted in Figure 1. Its importance is due to its role as a commercial and distribution center of regional importance (Al Sabahi et al. 2009). The area of Ibb city is about 5300 km² distributed over twenty districts, and its population is estimated at 3,243,580 people, or 10.8% of the total population of Yemen, at an annual rate of 2.47% (Dammag et al. 2023b). It is the third-largest city in terms of population and is inhabited by about 78% rural population. Furthermore, it is known as one of the most productive food warehouses in the country, with the city’s agricultural production accounting for (5.6%) of Yemen’s total production, making it significant not just ecologically, but also socio-economically (Varisco 2018). In recent years, the city has suffered unexpected eco-environment degradation as a result of rapid urbanization, population growth, and the influx of IDPs. Due to the mountainous areas surrounding the city, as well as the volcanic soil and heavy amounts of rain throughout the year, the city produces wheat, coffee, fruits, and vegetables, and is famous for its simple handicrafts (Ward 2014). The diversity of agricultural products is due to the heavy rainfall throughout the year, which makes the soil fertile (Dammag et al. 2023). Therefore, Ibb City holds significant importance as an agricultural hub and primary source of food for the inhabitants of the neighboring regions (Dammag et al. 2023a). The city's lands contain some minerals, the most important of which are clay minerals used in the manufacture of cement and refractory bricks, mineral (zeolite) used in the manufacture of detergents, basalt used in the manufacture of building stone, and the mineral used in the manufacture of ornamental stones (Berghof foundation 2020).

Figure 1. Location map of the study area.

Due to Ibb city location at an altitude of around 1,900 meters above sea level, it relishes a moderate climate throughout the year, and it is accompanied by heavy rains with coldness as a result of the monsoon rains saturated with water that blows from the southeast and southwest, as the annual precipitation rate is more than (1000 mm) on the western and southern mountainous heights at an altitude of (1500 m above sea level), with an annual mean temperature 18 ◦C (Https://weatherspark.com 2020). Geologically, it’s high in the southwest and low in the northeast, with undulating landforms, forming three types of terrain: middle mountains, low mountains, and river valleys (NIC 2015). The natural vegetation in the region depends on agricultural terraces, followed by marginal grasslands occasionally with trees (CEPF 2012). The area above 1500 m above sea level accounts for 49.40% of the administrative area. Most of these valley’s flow into the Tihama Plain in the west, while the valleys east of Ibb flow into the Gulf of Aden (Ward 2014).

2.2. Data acquisition and processing

The overall approach employed in this study (Figure 2) involved estimating ESV based on LULC change. Thus, LULC is the core data for this study, and the Landsat images were collected from the United States Geological Survey (USGS) official website (https://earthexplorer.usgs.gov)to elucidate LULC dynamics. Landsat images were comprised of Landsat TM (5) from 1990, ETM+ (7) from 2005, and Landsat OLI_TIRS (8) from 2020, which were interpreted for remote sensing images at a resolution of 30 m × 30 m. Track numbers 128/47 and 128/48 were obtained (Table 1). ERDAS's LULC change classification will use a 15-meter Google Earth connected to the software, allowing interpretation of the correct location of LULC for remote sensing blurry images. We selected Landsat data for the season of dry in each year because the cloud-free images were available during the season of dry. In addition, there was no rain during the season of dry, which ensures the actual presence of water bodies.

Figure 2. Flowchart of the methodology of this study.

Table 1. Detail of data sources used for analysis of LULC in the study.

Satellite	Sensor	Resolution (m)	Path/Row	Acquisition date	Season	Multispectral
Landsat 5	TM	30 × 30	128/47, 128/48	Feb 15, 1990	Dry	1–5 & 7
Landsat 7	ETM+	30 × 30	128/47, 128/48	Jan 17, 2005	Dry	1–7 & 8
Landsat 8	OLI-TIRS	30 × 30	128/47, 128/48	Feb 13, 2020	Dry	1–8, & 9

This study uses BTM to estimate ESV coefficients, which involves adjusting existing values or data from one location to estimate the ESV of other similar novel sites in the absence of site-specific evaluation information (Msofe et al. 2020). It is mainly useful when raw data cannot be collected due to financial or time restrictions (Aneseyee et al. 2022). The modification of ESV coefficients is performed based on BTM, which refers to the adaptation of existing values or data from one site to estimate ESV for other similar new sites in the absence of site-specific evaluation information (Msofe et al. 2020). This study considers the value of tropical areas with equivalent LULC classes for the geographical area studied to verify the significance of data extracted from the TEEB evaluation database.

ENVI 5.0 was used to perform atmospheric correction and geometric correction on the remote sensing images, and then supervised classification was used to interpret the remote sensing images. The socioeconomic data come from the ‘interviews’, and the precipitation data come from meteorological services. The process characteristics and reasons for LULC change from 1990 to 2050 at a 15-year interval were analyzed, including image processing methods, the creation of interpretation markers, and a remote sensing system for LULC classification. This study selected data from 1990 to 2005 as Phase I; Phase II data from 2005 to 2020; and Phase Final in 2035 and 2050 to analyze the impact of LULC change on the regional ESV in the study area.

2.3. LULC classification and simulation

Satellite images were pre-processed using ENVI 5.3 software for various corrections to ensure an accurate representation of biophysical phenomena. This pre-processing makes sure the imagery aligns correctly with the geographic coordinate system (Darem et al. 2023). Accurate land cover change detection requires meticulous radiometric and geometric corrections. Image classification involves sorting pixels in Landsat images into LULC types (Dammag et al. 2023). Specific bands from Landsat 5 TM and Landsat 8 OLI were used to produce false-color composite maps. The types of LULC were categorized into five classes, namely, cultivated land, grassland, built-up area, bare land, and waterbody by applying the maximum likelihood supervised classification method (Table 2). Then the process characteristics and reasons for LULC change during the study period 1990–2050 at a 15-year interval have been analyzed, and the spatiotemporal gain and loss for each type of land cover were calculated.

Table 2. Description of LULC classification classes.

LULC classes	Description
Built-up areas	Residential, industrial, and commercial areas, and infrastructure.
Barren lands	Exposed soils, un-vegetated lands, and bare land and rock outcrop.
Cultivated lands	Grow wooded crops, grain crops, and irrigated land.
Grassland	Herbaceous, scattered trees, and pastures.
Waterbody	Rivers, streams, lakes, dams, and various reservoirs.

The supervised maximum likelihood classifier (MLC) algorithm was applied post–pre-processing. This method is adaptable to various data and systems (Servadio and Convertino 2018). Some classifications faced challenges due to similar spectral characteristics, and post-classification refinement was employed for better accuracy. Mixed pixels, common in medium-resolution data, were corrected using additional data sources. For verification, 50 forms were created using public land use data and Google Earth maps. Then the accuracy of image classification maps was assessed to understand their reliability. An error/confusion matrix was employed a recommended method by many researchers (Jamal et al. 2022; Jodder et al. 2023), to compare categorized images with reference data. This matrix utilized reference points from multiple sources, including on-site observations, historical maps, and satellite imagery. For 2020, on-site observations with a 60CSx GPS map were used, encompassing an area of ∼100 m × 100 m (http://www.cso.gov.ye 2020). In total, 143 reference points formed the matrix, with specific allocations based on area coverage. This study used the CA-Markov model in GIS to simulate changes in LULC. The CA-Markov LULC matrix determined past trends and predicted future LULC configurations, transitions, and driving factors. This matrix calculated overall accuracy and the kappa coefficient (K) to measure the classification's accuracy, considering various types of errors. K values range from −1 to 1, indicating the type of land cover changes. Mathematically, the Kappa test measures the value of the coefficient in the point Equation (1) (1) $\mathrm{K}={\displaystyle \frac{\mathrm{S}{\sum }_{i=1}^r{\mathrm{A}}_{\mathrm{ii}}-{\sum }_{i=1}^r({\mathrm{A}}_{i+}\times {\mathrm{A}}_{+\mathrm{i}})}{{\mathcal{S}}^2-{\sum }_{i=1}^r({\mathrm{A}}_{i=1}\times {\mathrm{A}}_{+\mathrm{i}})}}$(1) Where K is the kappa coefficient, and S is the number of columns and land cover class in the error matrix. A_ii is the number of pixels in the first row and first column of the error matrix. A_i+ and A_+i are the total number of pixels in row i and column i, respectively. r is the number of rows in the matrix.

The results revealed that in 1990, the overall accuracy and kappa coefficient were 87.4% and 0.84, respectively. In 2005, these figures reached 86.4% and 0.82, while in 2020 they reached 89.2% and 0.86. The overall classification showed strong performance across all maps, reaching an overall accuracy of 87.6% and a resulting kappa coefficient of 0.84 (Table 3). Therefore, the image processing method applied in this study has been effective in producing time-matched LULC data.

Table 3. Accuracy assessment of LULC classified for 1990, 2005 and 2020.

LULC classes	1990		2005		2020
	PA (%)	UA	PA (%)	UA	PA (%)	UA
Cultivated land	84.3	90.7	93.4	88.7	82.7	93.3
Grassland	91.4	92.3	89.1	91.1	94.5	92.4
Built-up Area	88.3	80.0	85.7	90.4	92.5	91.3
Barren Land	85.5	86.4	80.1	85.7	85.5	86.9
waterbody	88.7	86.1	83.9	84.5	91.5	82.2
Overall accuracy %	87.4		86.4		89.2
Kappa Statistics	0.84		0.82		0.86

2.4. Assignment of ecosystem services values (ESV)

Most ESV evaluations use BTM based on universal or adjusted value coefficients developed by other researchers (Costanza et al. 1997; Pham and Lin 2023). This study has used the BTM because of its simplicity and parsimony in the approach and resources required for assessment, particularly in data-poor regions (Woldeyohannes et al. 2020). Although value proposition coefficients have been criticized due to the limitations and uncertainties in their local use, several researchers have used them in areas where data are scarce and improved upon initial values (Costanza et al. 1997) for ecosystem service assessment (Liu et al. 2019; Costanza et al. 2014). Furthermore, the availability of global datasets on LULC based on remote sensing and regional datasets based on validation and/or correction of the global data facilitates the application of BTM (Woldeyohannes et al. 2020). Therefore, we used this method in our study to evaluate ecosystem services for Ibb City. To calculate ESV associated with different types of LULC changes, we conducted a comparison between the LULC classifications obtained from the Landsat TM/ETM dataset and the biomes identified in the ES Assessment Model developed by Costanza et al. (2014). The most relevant biome for each category was assigned as the proxy for that LULC type. For example, cropland biome for ‘Cultivated land’, desert land for ‘Bare land’, urban/village area for ‘Built-up area’, grass/rangelands ‘Grassland’, and lakes/rivers ‘Waterbody’. Despite the LULC categories and equivalent biomes do not precisely match, their use has been possible in many other studies (Morshed et al. 2022; Msofe et al. 2020; Rahman and Szabó 2021). LULC types, equivalent biome, and corresponding value coefficient over historical years 1990, 2005, and 2020, and projected years 2035 and 2050 for the study area are presented in Table 4.

Table 4. Corresponding biomes for LULC classes and ES coefficients based on the adjusted estimates.

LULC types	Equivalent biome	Ecosystem service coefficient (US$ ha^−1 year^−1)
Cultivated land	Cropland	225.6
Grassland	Grass/rangelands	293.3
Built-up Area	Urban	0
Bare Land	Desert	0
Waterbody	Lakes/rivers	8103.5

Source: Costanza et al. (1997).

In order to assess changes in ESV, all value coefficients were converted to US dollars (US$) per hectare per year, based on the 1994 currency value. Moreover, the study incorporated adjustments to the value coefficients using the Consumer Price Index (CPI) and the Producer Price Index (PPI) to account for changes in inflation and price levels over time (available at http://data.imf.org/) to analyze the impact of the temporal evolution of the coefficients on the total estimate of ESV change. In Table 5, we present specific information regarding the annual ecosystem service value coefficients adjusted for each LULC change scenario. The study area maps were created using other datasets in a GIS. Therefore, these unit values were more reliable for calculating ESV at the regional level.

Table 5. Details of coefficients for annual conservative adjustments and ecosystem service functions for each LULC type in 1994 US$ ha−1 year−1.

Ecosystem services	LULC change types/Biome
	Cropland	Grassland	Urban	Waterbody
Provisioning services	181.2	119.6	0.0	2581.2
Food protection	181.2	94.1		157.8
Water supply				2423.4
Raw material		24.3
Genetic resources
Medical services		1.2
Regulating services	27.2	133.1	0.0	5452.7
Climate regulation		58.1
Disturbance regulation
Water regulation				5145
Waste treatment				307.7
Erosion control		74.3
Gas regulation
Biological control	27.2	0.7
Supporting services	17.2	29.6	0.0	0.0
soil formation		10.2
Nutrient cycle
Pollination	17.2	19.4
Habitat/refugia
Cultural services	0.0	11	0.0	69.6
Recreation		5.3		69.6
Cultural		5.7
Sum economic ESVs	225.6	293.3	0.0	8103.5

2.5. Calculation of ESV

Considering the LULC class, we made adjustments to the coefficient of ESV based on the fundamental equivalent table of ESV per unit area, as proposed by Yang et al. (2020), and Costanza et al. (1997; 2014). The total value of ecosystem services in a particular period, denoted as ESVt was determined by the following equation (2) (2) $\mathrm{ES}{V}_t=\sum _{k=1}^n(A_{\mathrm{kt}}\times V{C}_k)$(2) Where ESVt is the total value of ecosystem services in an estimated period, Akt represents the area of LULC type ‘k’ (ha), and VCk is the ecological value coefficient (US$ ha−1 year−1) associated with LULC type ‘k’. System service value coefficient, ‘n’ is the number of LULC types (n = 5).

To calculate the change rate of ecosystem service value (ESVcr) was determined using the following equations Equation (3) (3) $\mathrm{ES}{V}_{\mathrm{cr}}={\displaystyle \frac{\mathrm{ES}{V}_{t2}-\mathrm{ES}{V}_{t1}}{\mathrm{ES}{V}_{t1}}\times 100\text{\%}}$(3) Where, ‘ESVt1’ represents the total value of ecosystem services in the initial year ‘t1’ and ‘ESVt2’ in the recent year ‘t2’, respectively. A positive ESVcr indicates an increase in the amount of ecosystem service value, while a negative value suggests a decrease.

To estimate the secondary ecosystem service values (ESVft) for a specific item ‘f’ within a given period ‘t’, was calculated using the equation below Equation (4) (4) $\mathrm{ES}{V}_{\mathrm{ft}}=\sum _{k=1}^n(V{C}_{\mathrm{fk}}\times A_{\mathrm{ft}})$(4) Where VC_fk represents the value per unit area of ecosystem service ‘f’ associated with LULC type ‘k’, ‘m’ represents the number of ecosystem service types, and ‘Aft’ represents the area (ha) of LULC type ‘k’ contributing to the ecosystem service ‘f’.

2.6. Calculating elasticity of ESV in response LULC changes

Describes elasticity as a measure that quantifies the responsiveness of one variable to changes in another variable (Yuan et al. 2019). We will use this approach to understand how the total ESV changes with changes in LULC by calculating the elasticity of the ESV generated by measuring the percentage change in the ESV in response to the percentage change in LULC (F. Zhang, Yushanjiang, and Jing 2019), using the following Equations (5(5) $\mathrm{EL}={\displaystyle \frac{\mathrm{ES}{V}_{\mathrm{cr}}}{\mathrm{LTP}}}$(5) ) and (6(6) $\mathrm{LTP}=\sum _{k=1}^n{\displaystyle \frac{\mathrm{△}{A}_k}{A_k}\times {\displaystyle \frac{1}{T}\times 100\text{\%}}}$(6) ) (5) $\mathrm{EL}={\displaystyle \frac{\mathrm{ES}{V}_{\mathrm{cr}}}{\mathrm{LTP}}}$(5) (6) $\mathrm{LTP}=\sum _{k=1}^n{\displaystyle \frac{\mathrm{△}{A}_k}{A_k}\times {\displaystyle \frac{1}{T}\times 100\text{\%}}}$(6) Where EL is the elasticity of change of total ESV concerning LULC change, ESV_cr was explained above Equation (3(3) $\mathrm{ES}{V}_{\mathrm{cr}}={\displaystyle \frac{\mathrm{ES}{V}_{t2}-\mathrm{ES}{V}_{t1}}{\mathrm{ES}{V}_{t1}}\times 100\text{\%}}$(3) ), LTP is the percentage of land conversion denoting LULC change degree, $\mathrm{\Delta }$Ak is the converted area of LULC class ‘k’, and L_i is the initial area of LULC type ‘k’.

2.7. Spatial autocorrelation analyses

The spatial autocorrelation model used in this study mainly represents the degree of spatial concentration of data and is conducive to fully understanding ESV's changing trends and internal mechanisms (F. Zhang, Yushanjiang, and Jing 2019). Global spatial autocorrelation is usually represented by Global Moran’s I, which can reveal overall whether the spatial distribution of a certain attribute in a region is correlated with neighboring areas and the degree of correlation (L. Li et al. 2022), but it is difficult to measure local spatial agglomeration and its spatial heterogeneity. To this end, it is still necessary to use Local Moran’s I to characterize the spatial differentiation characteristics of coupling coordination from the local scale (Shi et al. 2022). The Equations are shown respectively (7) $\mathbf{I}={\displaystyle \frac{n{\sum }_{i=1}^n{\sum }_{i=1}^n{w}_{\mathrm{ij}}(x_i-x^{-})(x_i-x^{-})}{{\sum }_{i=1}^n{\sum }_{i=1}^n{w}_{\mathrm{ij}}{\sum }_{i=1}^n{(x_i-x^{-})}^2}}$(7) (8) ${\mathbf{I}}_{\mathit{i}}={\displaystyle \frac{n(x_i-x^{-}){\sum }_{i-1}^n{w}_{\mathrm{ij}}(x_i-x^{-})}{{\sum }_{i=1}^n{(x_i-x^{-})}^2}}$(8) Where Ii and I represent the global Moran index and the local Moran index respectively; n is the number of spatial units; xi and xj are the observed values of coupling coordination degree of prefecture i and j respectively; x represents the average coupling coordination degree of the study area; Wij is the spatial weight matrix of prefectures i and j, a Queen first-order weight matrix formed by Geoda; the I value range is [−1,1], when I > 0 it is a positive spatial correlation, and the element attributes are clustered and distributed; when I < 0 It is a negative spatial correlation, and the feature attributes present a discretized distribution; when I = 0, it is a random distribution. According to the local Moran index, the coupling coordination level can be divided into four types: HH (High-High clustering), HL (High-Low clustering), LH (Low-High clustering), and LL (Low-Low clustering).

This study uses the Weight tool in GeoDa software to create a spatial weight matrix in Ibb City. It uses the global Moran’s I index and local LISA agglomeration map to analyze the spatial aggregation differences of ESV in Ibb City.

2.8. Sensitivity analysis

The degree of ESV is based on the value of ecological service functions by using the coefficient of sensitivity (CS) to determine whether an equivalent factor is appropriate (Su, Wei, and Lin 2020). In this study, the sensitivity index for each land type was determined by adjusting the equivalent factor specific to that land type by 50%, to evaluate the resulting change in the overall value of ecological services and to determine the coefficient of sensitivity for the given land type (N. Li et al. 2021). To confirm the representativeness of the ecosystem type for each land cover type and the accuracy of the obtained value coefficient (economic value per unit area), the CS was used to measure the accuracy of the value coefficient. If CS > l, it indicates that the estimated ESV is characterized by elasticity concerning the ecological value coefficient, and the accuracy is relatively low. Conversely, if CS < l, the ecosystem value is considered to be inelastic, indicating that the ecological value coefficient is precise and reliable (Haque et al. 2023). The higher the property is, the more accurately the ESV can be estimated. We then calculated the CS by equation (9(9) $\mathrm{CS}=\left|{\displaystyle \frac{(\mathrm{ES}{V}_j-\mathrm{ES}{V}_i)/\mathrm{ES}{V}_i}{{\displaystyle \frac{V{C}_{\mathrm{jk}}-V{C}_{\mathrm{ik}}}{V{C}_{\mathrm{ik}}}}}}\right|$(9) ) (Belay, Melese, and Senamaw 2022; Hu et al. 2020) (9) $\mathrm{CS}=\left|{\displaystyle \frac{(\mathrm{ES}{V}_j-\mathrm{ES}{V}_i)/\mathrm{ES}{V}_i}{{\displaystyle \frac{V{C}_{\mathrm{jk}}-V{C}_{\mathrm{ik}}}{V{C}_{\mathrm{ik}}}}}}\right|$(9) Where, ESV refers to an estimation of the overall ESV offered by a specific LULC type, whereas VC serves as an indicator for comparing the relative value of various land use types. ESVi and ESVj represent the estimated total value of ecosystem services before and after adjusting the ecological value coefficient for the LULC type in question. The variable K represents each LULC type.

3. Results

3.1. LULC changes

This study investigated quantifying the LULC dynamics from 1990 to 2020, and we used the LULC base map from 2005 and transition probabilities from 1990 to 2005 and chose the number of cycles 15 years to simulate the LULC for 2020 by using the CA-Markov model. Following the above process, the future LULC in 2035 and 2050 was predicted. The results showed an overall accuracy of 85.5% and a Kappa coefficient index of 0.80 between the actual LULC and the simulated LULC maps of 2020, which is a good agreement according to (Dey et al. 2021). Moreover, we conducted a validation to ascertain the significance of the association between actual and simulated LULC values for each category in 2020, employing the Chi-Square (χ²) test. The results yielded a χ² value of 26.87 for a degree of freedom (df) equal to 7. The result showed that the validation unit includes five statistical coefficients for kappa: kappa no (K_no) with a value of 0.9132, kappa no (K_location) with a value of 0.8943, kappa location strata (K_{location Strat}) with a value of 0.9304, and standard kappa (K_standard) with a value of 0.8426. This shows that the CA-Markov model has high accuracy to simulate future LULC changes in the study area based on pre-processed data and factors that have an impact on land use changes, such as natural, social, and economic policies. To reduce randomness in the simulation process, this study defines conversion rules for land use types along with influencing factors such as natural conditions, planning policies, climate change, the eco-environment, and other influencing factors. Table 6 and Figures 3 and 4 illustrates the patterns of spatial distribution of LULC from 1990 to 2050 in the study area. From 1990 to 2005, the built-up areas, grassland, barren land, and waterbody increased by 11.8, 56, 359.1, and 2.2 km², at an annual rate of 2.4, 3.7, 23.9, 0.06 km²/year, respectively. Conversely, vegetated land decreased by 451.7 km² with a rate of 4.2 km²/year. During (2005–2020), the barren land decreased by 121 km², at an annual rate of change of 8.1 km²/year. On the other hand, the vegetated land expanded by 63 km², with a rate of 4.2 km²/year, also grassland increased by 27 km², at an annual change rate of 1.8 km². Built-up areas continued to grow, adding 29 km² (1.9 km²/year), and the waterbody area increased by 1.6 km² (0.1 km²/year). The eco-environment has improved due to some environmental restoration projects implemented by some organizations in the region.

Figure 3. Spatial distribution of LULC types in the study area from 1990 to 2050.

Figure 4. Dynamic degree of LULC types of the study area from 1990 to 2050.

Table 6. Area and proportion of different LULC classes.

LULC classes	1990		2005		2020		2035		2050
	Area (km^2)	%	Area (km^2)	%	Area (km^2)	%	Area (km^2)	%	Area (km^2)	%
Vegetation land	1760.4	33.2	1308.7	24.7	1371.7	27.8	1334.9	25.2	1286.7	24.2
Grassland	336.0	6.3	392.3	7.4	419.3	7.9	408.5	7.7	397.6	7.5
Built-up Area	11.8	0.38	47.2	0.94	76.2	1.42	126.4	2.4	184.7	4.1
Barren Land	3190.5	60.1	3549.6	66.9	3428.6	64.8	3424.8	64.6	3423.7	64.5
Waterbody	1.3	0.02	2.2	0.06	3.8	0.08	5.4	0.1	7.3	0.2

The study area has experienced a notable increase in average temperature by 0.25°C per decade, the highest in the area, leading to a decline in biodiversity and vegetation cover. While annual precipitation remained relatively stable in the early 1990s, it approached the multi-year average in the 2000s, peaking in 2005 and 2020, and dropping significantly in 2010. This decline in rainfall has led to water deficits, a critical factor in the degradation of vegetation land, and the increase of barren land. Population trends in the region have shown an increasing trend due to the return of migrants after the Gulf War in 1990, with a growth rate of 4.5% in the first period and 4.4% in the second. Since 2015, population overcrowding due to war and the continuous influx of displaced people has brought about substantial changes in the economic and social fabric, mode of production, and way of life. This shift towards a non-agricultural population has led to a 7.4% decrease in the proportion of agricultural GDP in the region. In addition, national policies and institutions have played a role in shaping LULC changes, where some areas lack a system for registering land titles and documents.

In 2035, cropland, grassland, and barren lands are predicted to decrease totaling a reduction of 36.8, 10.8, and 249.3 km², respectively. Conversely, the built-up area is projected to increase by 50.2 km², driven by rapid population growth in Ibb City, primarily due to rural-urban migration and IDPs. Additionally, the volume of water bodies is anticipated to increase at 1.6 km² due to the expansion of dams and water barriers as outlined in the Water Resources Office's future plan. The decline in vegetation, grasslands, and barren lands is expected to persist in 2050, reaching 1286.7, 397.6, and 3423.7 km², respectively. However, the built-up area is predicted to increase by 184.7 km² and waterbody by 1.9 km².

3.2. Variation in ecosystem services value

3.2.1. Estimation of changes in ESV

The changes in ESV were evaluated from 1990 to 2020 periods, and its predicted amount in the years 2035 and 2050, the total ESV in Table 7, and the area of different LULC types in Figure 5. Please note that we have converted the land area unit of measurement from km² to hectares since calculating the coefficients of ESV. It can be observed in Figure 5 that the total estimated ESV has reduced from US$ 68.5 million in 1990 to US$ 66.8 million in 2005 and US$ 65.2 Mn in 2020, and is predicted to further decrease to US$ 63.8 Mn in 2035 and US$ 61.2 million in 2050.

Figure 5. Trend of total ESV from 1990 to 2050.

Table 7. Estimated ESV for each LULC type of the different reference years (US$ Million/year).

LULC classes	1990 ESV	%	2005 ESV	%	2020 ESV	%	2035 ESV	%	2050 ESV	%
Cultivated land	4.9	7.2	4.2	6.3	4.6	7.1	3.7	5.8	2.8	4.6
Grassland	15.7	22.9	12.4	18.6	9.2	14.1	7.5	11.8	4.6	7.5
Built-up area	0.00	0	0.00	0	0.00	0	0.00	0	0.00	0
Bare Land	0.00	0	0.00	0	0.00	0	0.00	0	0.00	0
Waterbody	47.9	69.9	50.2	75.1	51.4	78.8	52.6	82.4	53.8	87.9
Sum	68.5	100	66.8	100	65.2	100	63.8	100	61.2	100

Table 7 clearly shows that water bodies were the highest contributor to total ESV, followed by grassland, and vegetation land 1990 and 2020 periods. On the other hand, built-up areas and bare land did not contribute to ESV. According to Figure 6, the ESV trend shows a declining change process. A gradual decrease in ESV was observed for vegetation land, and grassland over the study period 1990–2005 and 2005–2020. A similar trend in the provision of services is predicted to continue from 2020 to 2035, except that the provision of cultivated land services has consistently decreased and ranked third. The sudden decline in cultivated land and its environmental and social value is attributed to heavy commercial investments in the area. It is also projected that waterbodies will continue to contribute the most to total ESV by 2050, followed by grasslands and croplands. Since waterbodies were the highest contributor to total ESV, the characteristics of the total ESV change process are dominated by the characteristic change process of ESV derived from waterbodies.

Figure 6. Spatial distribution of ESV changes in the study area from 1990 to 2050.

According to the result, as presented in Table 8 and Figure 6, the annual rate of changes of ESV for different LULC classes during the period 1990–2050. A negative sign in the value indicates a decrease in the ESV. In 1990, total changes to ESV over the first study period were approximately US$ 1.7 Mn (2.5%) of the present value in 1990. The total value of ecosystem services experienced a decrease of US$1.6 million (2.4%,) compared to the year 1990. As previously stated, the total ESV underwent a change of about US$3.3 million between 1990 and 2020, which signifies a 4.8% decrease from its 1990 value. Due to the encroachment of built-up areas on cropland and grassland. On the contrary, the values of the water body showed an increase of approximately US$2.3 million (4.8%). Meanwhile, changes in the total value of ecosystem services are expected to continue to decline by approximately US$ 1.4 million (2.2%) and US$ 2.6 million (4.1%) in 2035 and 2050. As a result, the total value of ES over the full prediction period decreased by US$4.0 million, which means 6.1% of the value that existed in 2020.

Table 8. Total ecosystem service valuation (ESV) (US$ in millions) estimated for each LULC classes.

LULC classes	ESVcr
	1990–2005		2005–2020		1990–2020		2020–2035		2035–2050		2020–2050
	Million US$	%	Million US$	%	Million US$	%	Million US$	%	Million US$	%	Million US$	%
Cultivated land	−0.7	−14.3	0.4	9.5	−0.2	−6.1	−0.9	−19.6	−0.8	−24.3	−1.8	−39.1
Grassland	−3.3	−21.0	−3.2	−25.8	−6.4	−41.4	−1.7	−18.5	−2.9	−38.7	−4.6	−50.0
Built-up area	0.0	0	0.0	0	0.0	0	0.0	0	0.0	0	0.0	0
Bare land	0.0	0	0.0	0	0.0	0	0.0	0	0.0	0	0.0	0
Waterbody	2.3	4.8	1.2	2.4	3.3	7.3	1.2	2.3	1.2	2.3	2.4	4.7
Sum	−1.7	−2.5	−1.6	−2.4	−3.3	−4.8	−1.4	−2.2	−2.6	−4.1	−4.0	−6.1

Although the ESV change for each LULC class shows a significant decrease in the value of cropland and grassland, the ESV change for built-up and bare land is absent because they do not derive the ESV. It is also noted that the LULC classes that occurred during the historical periods from 1990 to 2020 have lost US$ 3.3 Mn, and the projected changes in LULC across the periods 2020–2050 will decrease by US$ 4.0 million if the current condition is not enhanced in terms of land management practices (Table 8). Interestingly, the percentage changes in the estimated ESV for each type of LULC change were comparable across the study period. However, the primary factor contributing to the overall decrease in the total ESV over the entire study period was the conversion of a significant area of cultivated land and grassland into construction land. This conversion was driven by economic development and the indiscriminate implementation of construction projects in the study area. This underscores the significance of conserving and effectively managing natural resources, along with promoting sustainable LULC practices, to uphold and enhance the provision of ecosystem services.

Overall, the projected trends in ecosystem service values suggest that there is a need for continued efforts to protect and manage natural resources in the study area, to ensure the provision of important environmental services, and to support sustainable development.

3.2.3. Changes in the value of individual functions of ecosystem services

In this study, the general trends in the estimated value of individual ecosystem functions (ESV_ft) were calculated for each specific period, and their aggregate changes were based on Equation (4(4) $\mathrm{ES}{V}_{\mathrm{ft}}=\sum _{k=1}^n(V{C}_{\mathrm{fk}}\times A_{\mathrm{ft}})$(4) ). In addition, we calculated the effect of LULC changes on individual ecosystem function change over the historical (1990–2020) and projected (2020–2050) period at a 15-year interval; their total contribution represents approximately US$ 68.5, 66.8, 65.2, and 63.8 61.2 million in 1990, 2005, 2020, 2035, and 2050, respectively.

These 17 ecosystem services were divided into four main ecosystem services: regulating (seven types), provisioning (five types), supporting (four types), and cultural (two types). Understanding how each ecosystem function contributes to the overall ESV is made easier with the help of the value of individual ecosystem functions (Rahman and Szabó 2021). It was revealed that out of the 17 ecosystem services included in the study, 14 of them showed a significant decrease when measured by Equation (4(4) $\mathrm{ES}{V}_{\mathrm{ft}}=\sum _{k=1}^n(V{C}_{\mathrm{fk}}\times A_{\mathrm{ft}})$(4) ). Table 9 shows that throughout the study period, regulating services made the largest contribution to the total ESV, followed by provisioning, supporting, and cultural services.

Table 9. Estimated annual value of ecosystem functions (ESV_f in US $ million per year).

Ecosystem services function	1990		2005		2020		2035		2050
	ESVf	%	ESVf	%	ESVf	%	ESVf	%	ESVf	%
Provisioning services	21.4	31.2	20.7	31.0	20.5	31.4	20.4	32.0	19.9	32.5
Food production	5.8	8.5	7.4	11.1	8.3	12.7	9.1	14.3	9.8	16.0
Water supply	4.9	7.2	4.2	6.3	3.8	5.8	3.4	5.3	3.1	5.1
Raw material	3.8	5.5	3.1	4.6	2.9	4.4	2.6	4.1	2.2	3.6
Genetic resources	2.6	3.8	2.2	3.3	2.1	3.2	2.1	3.3	1.8	2.9
Medical services	4.3	6.3	3.8	5.7	3.4	5.2	3.2	5.0	3.0	4.9
Regulating services	27.5	40.1	26.8	40.0	25.8	39.6	24.9	39.0	23.2	38.0
Climate regulation	3.8	5.5	3.7	5.5	3.2	4.9	2.9	4.5	2.7	4.4
Disturbance regulation	2.9	4.2	2.7	4.0	2.6	4.0	2.3	3.6	2.1	3.4
Water regulation	3.6	5.3	3.5	5.2	3.2	4.9	3.0	4.7	2.7	4.4
Waste treatment	4.8	7.0	4.3	6.4	4.1	6.3	3.8	6.0	3.4	5.6
Erosion control	5.7	8.3	4.9	7.3	4.6	7.1	4.3	6.7	4.1	6.7
Gas regulation	3.9	5.7	3.8	5.7	3.4	5.2	3.2	5.0	1.9	3.1
Biological control	2.8	4.1	3.9	5.8	4.7	7.2	5.4	8.5	6.3	10.3
Supporting services	14.9	21.8	15.3	22.9	15.4	23.6	15.5	24.3	15.6	25.5
soil formation	4.7	6.9	4.5	6.7	4.3	6.6	4.1	6.4	3.7	6.0
Nutrient cycle	3.2	4.5	2.8	4.2	2.5	3.8	2.1	3.3	1.9	3.1
Pollination	4.3	6.3	5.4	8.1	6.3	9.7	7.2	11.3	8.2	13.3
Habitat/refugia	2.7	3.9	2.6	3.9	2.3	3.5	2.1	3.3	1.8	2.9
Cultural services	4.7	6.9	4.0	6.0	3.5	5.4	3.0	4.7	2.5	4.1
Recreation	2.9	4.2	2.6	3.9	2.2	3.4	1.9	3.0	1.7	2.8
Cultural	1.8	2.6	1.4	2.1	1.3	2.0	1.1	1.7	0.8	1.3
Total	68.5	100	66.8	100	65.2	100	63.8	100	61.2	100

In 1990, regulating services 40.1% accounted for the largest annual share of total ESV, followed by provisioning, support, and cultural services which means 31.2%, 21.8%, and 6.9%, respectively, yet, supporting services increased to 22.9% and 23.6%, in 2005 and 2020. On the contrary, provisioning, regulating, and culture services decreased to 31.0%, 40.0%, and 6.0% in 2005; 31.4%, 39.6%, and 5.4% in 2020. The increase of supporting services is predicted to continue at about 24.3% and 25.5% in the periods 2035 and 2050, while provisioning, regulating, and cultural services will continue to decline which means 32.5%, 38.0%, and 4.1%, respectively, in 2050. Table 9 clearly shows that the most significant downward trend in total ESV was regulating services (US$ 1.7 million), followed by cultural services (US$ 1.2 Mn), and provisioning services (US$ 0.9 million) in 2020. The contributions of cultural and Supporting services are predicted to decrease by (US$ 1.0 and 0.6 Mn), but the regulating services continue to increase by (10.1%) across the period 2020–2050 (Figure 7).

Figure 7. Contribution percentages of an individual ecosystem service to total ESV.

Over the periods 1990–2020, the ecosystem service functions of food production, biological control, and pollination contributed significantly to the total economic values of ESV concerning LULC changes which mean 43.1% (US$ 2.5 million), 67.9% (US$ 1.9 million), and 46.5% (US$ 2.0 Mn), respectively. On the contrary, most of the ecosystem services were reduced as follows, from (high to low): water support (US$ 1.1 million ≈ 22.4%), erosion control (US$ 1.1 Mn ≈ 19.3%), and raw materials (US$ 0.9 million ≈ 23.7%), medical services (US$ 0.9 million ≈ 20.9%), waste treatment (US$ 0.7 million ≈ 14.6%), and climate regulation (US$ 0.6 million ≈ 15.8%), as shown in Table 10.

Table 10. Change rate in value (US$ million) of ecosystem service functions from 1990 to 2050.

Ecosystem services function	Overall change
	1990–2005	2005–2020	1990–2020	2020–2035	2035–2050	2020–2050
Provisioning services	−0.7	−0.2	−0.9	−0.1	−0.5	−0.6
Food production	1.6	0.9	2.5	0.8	0.7	1.5
Water supply	−0.7	−0.4	−1.1	−0.4	−0.3	−0.7
Raw material	−0.7	−0.2	−0.9	−0.3	−0.4	−0.7
Genetic resources	−0.4	−0.1	−0.5	0	−0.3	−0.3
Medical services	−0.5	−0.4	−0.9	−0.2	−0.2	−0.4
Regulating services	−0.7	−1.0	−1.7	−0.9	−1.7	−2.6
Climate regulation	−0.1	−0.5	−0.6	−0.3	−0.2	−0.5
Disturbance regulation	−0.2	−0.1	−0.3	−0.3	−0.2	−0.5
Water regulation	−0.1	−0.3	−0.4	−0.2	−0.3	−0.5
Waste treatment	−0.5	−0.2	−0.7	−0.3	−0.4	−0.7
Erosion control	−0.8	−0.3	−1.1	−0.3	−0.2	−0.5
Gas regulation	−0.1	−0.4	−0.5	−0.2	−1.3	−1.5
Biological control	1.1	0.8	1.9	0.7	0.9	1.6
Supporting services	0.4	0.1	0.5	0.1	0.1	0.2
soil formation	−0.2	−0.2	−0.4	−0.2	−0.4	−0.6
Nutrient cycle	−0.4	−0.3	−0.7	−0.4	−0.2	−0.6
Pollination	1.1	0.9	2.0	0.9	1.0	1.9
Habitat/refugia	−0.1	−0.3	−0.4	−0.2	−0.3	−0.5
Cultural services	−0.7	−0.5	−1.2	−0.5	−0.5	−1
Recreation	−0.3	−0.4	−0.7	−0.3	−0.2	−0.5
Cultural	−0.4	−0.1	−0.5	−0.2	−0.3	−0.5
Total	−1.7	−1.6	−3.3	−1.4	−2.6	−4.0

Although the contribution to total ESV for food production, biological control, and pollination remained increasing over the period 1990–2020, their contributions are predicted to decrease i.e. 18.1% (US$ 1.5 million), 34.1% (US$ 1.6 million), and 30.2% (US$ 1.9 Mn), respectively, in 2050. The highest declines in ESV are expected for gas regulation (US$ 1.5 million), water supply (US$ 0.7 Mn), waste treatment (US$ 0.7 million), and raw materials (US$ 0.7 million) over 2020–2050. In the case of ESV supportive and cultural groups, the nutrient cycle and recreation function had a negative role in the total ESV with apparently equal contributions of (US$ 0.7 Mn) over the period 1990–2020. However, their contribution is expected to decrease by (US$ 0.6 million), and entertainment (US$ 0.5 million) in 2050. As mentioned earlier, grassland destruction as a result of uncontrolled urbanization has had a significant impact on habitat services. This impact is directly related to the food supply, which in turn affects biological control and pollination services (Morshed et al. 2022; F. Zhang, Yushanjiang, and Jing 2019). Additionally, the estimated erosion control function of ecosystem services is compelling, given the challenging topography of the region, particularly in the southern part. The fact that most land uses are agricultural may indicate adaptations of populations to a rough terrain environment, particularly erosion control.

3.2.4. Impact of LULC change on ESV

The ESV resilience was calculated for each LULC type using Equation (5, 6) to show the effects of LULC changes on the total ESV. The results revealed that LULC change had varied effects on different ES throughout the study periods. It can be observed that certain ESV experienced an increase, while others underwent a decrease as a consequence of LULC change as illustrated in Figure 8. The elasticity implies that a 1% conversion in land use would result in an average change of US$ 291.7 and 203.8 million in total ESV over the period 1990–2020 and 2020–2050.

Figure 8. Contributing percentages of LULC types (%) and ESV (%) for each study period.

The matrix revealed that the area affected by LULC changes was 777.0 and 223.6 km², respectively, over 1990–2020 and 2020–2050. Across the period 1990–2020, cultivated land, which was mainly the result of urbanization, decreased by 3.4 km² and led to a decrease of US$ 1.3 million. The conversion of cultivated land to grassland and bare land, with an approximate area of 32.7 and 97.5 km² respectively, had a significant impact on the estimated ESV. This LULC change resulted in a decrease in ESV by approximately US$24.7 million for the conversion to grassland and US$68.2 million for the conversion to bare land. Furthermore, the decrease in grassland area contributed to a reduction in ESV by US$83.4 Mn. Conversely, 6.8% and 61.7% of built-up land underwent conversion to cropland and grassland, increasing from US$13.4 and US$13.2 million respectively, as presented in Table 11.

Table 11. Transition matrix of land use and ESVs in the study area from 1990 to 2020 (km²)

LULC classes	LULC Transfer (km^2)
	Cultivated land	Grassland	Built-up Area	Bare Land	Waterbody	Total 2020
Initial State (1990)	Final State 2020
Cultivated land	1238.1	32.7	3.4	97.5	0.0	1371.7
Grassland	136.7	270.6	2.3	9.5	0.2	419.3
Built-up Area	45.4	18.1	4.4	7.9	0.4	76.2
Bare Land	338.6	13.9	1.7	3074.4	0.0	3428.6
Waterbody	1.2	0.7	0.0	1.2	0.7	3.8
Total 1990	1760.4	336.0	11.8	3190.5	1.3	5300
ESVs transfer (US$ million)
Cultivated land	0.0	−24.7	−1.3	−68.2	0.0	−94.2
Grassland	−78.3	0.0	−0.8	−4.7	0.4	−83.4
Built-up Area	13.4	13.2	0.0	3.8	0.3	30.7
Bare Land	−132.1	−9.7	−0.6	0.0	0.0	−142.4
Waterbody	−1.4	−0.2	0.0	−0.8	0.0	−2.4
Total	−198.4	−21.4	−2.7	−69.9	0.7	−291.7

During the predicted period 2020–2050, the transition matrix showed a series of expected decreases in cultivated land by 85.0 km² and grassland by 21.7 km² resulting in a decrease in total ESV by US$ 117.0 and 42.8 million, which is the most significant change in the type of transition that predicts the decline in ESV. Table 12 clearly reports a decrease in agricultural land is also expected by 85 km², due to the conversion of 113.2 km² into bare land and 61.1 km² to grassland, resulting in a decrease of US$ 78.3 and 32.6 million. Meanwhile, converting 50.1% and 13.7% of the built-up area were into cropland and grassland, an increase of US$ 28.3 and 17.4 million, respectively. At the same time, grassland would decrease due to the conversion of 39.5 km² to bare land and 9.4 km² to built-up land resulting in a reduction of US$ 6.7 and 2.6 million in 2050. This fully reflects the enormous impact of LULC changes caused by human activities on ESV. Therefore, it is recommended to implement strict control measures in the conversion of bare land into construction land in future development plans, maintain the existing area of cultivated land, and consider strategies to increase the area of water bodies.

Table 12. Transition matrix of land use and ESVs in the study area from 2020 to 2050 (km²).

LULC Transfer (km^2)
LULC classes	Vegetation land	Grassland	Built-up Area	Bare Land	Waterbody	Total 2050
Initial State (2020)	Final State 2050
Vegetation land	1093.3	61.8	18.4	113.2	0	1286.7
Grassland	40	307.6	9.4	39.5	1.1	397.6
Built-up Area	92.6	25.3	36.7	28.8	1.3	184.7
Bare Land	144.7	23.9	11.7	3243.4	0	3423.7
Waterbody	1.5	0.7	0	3.7	1.4	7.3
Total 2020	1371.7	419.3	76.2	3428.6	3.8	5300
ESVs transfer (US$ million)
Vegetation land	0.0	−32.6	−6.1	−78.3	0.0	−117
Grassland	−34.7	0.0	−2.6	−6.7	1.2	−42.8
Built-up Area	28.3	17.4	0.0	8.6	0.9	55.2
Bare Land	−76.9	−13.2	−4.2	0.0	0.0	−94.3
Waterbody	−2.6	−0.9	0.0	−1.4	0.0	−4.9
Total	−85.9	−29.3	−12.9	−77.8	2.1	−203.8

3.5. Spatial autocorrelation results

3.5.1. Global spatial autocorrelation of ESV

In order to better show the spatial evolution characteristics of ESV, this study used the ‘Queen’ adjacency method, a Moran scatterplot was created to obtain a spatial weight matrix to analyze the ESV of Ibb City from 1990 to 2020, and predicted from 2020 to 2050. The scattered points were mainly distributed in the HH and LL quadrants as shown in Figure 9. Table 13 indicates that Moran's I index of ESV is 0.62, 0.67, and 0.69, and the predicted is 0.73, and 0.78 respectively. After selecting 999 replacements for the number of randomization iterations, the calculated Z scores are all greater than 2.58, and the p values are all 0.01. At the same time, Moran’s I values are all positive, indicating that ESV has global spatial aggregation characteristics and is constantly increasing.

Figure 9. Moran scatter plots of the ESV in the study area.

Table 13. ESV global spatial autocorrelation indices for 1990–2050.

Moran index	Different Years
	1990	2005	2020	2035	2050
Moran’s I	0.62	0.67	0.69	0.73	0.78
Value (P)	0.001	0.001	0.001	0.001	0.001
Score (Z)	8.65	9.34	9.62	10.18	10.88

3.5.2. Local spatial autocorrelation of ESV

A five-phase LISA agglomeration map was obtained by further exploring ESV's local spatial agglomeration characteristics in the study area. It can be seen from the local LISA map in Figure 10 that the insignificant areas are mainly distributed in plain areas, ESV had spatial differences during the study time point, and the agglomeration characteristics showed an LL Mainly agglomeration. In 1990, the study area was still a traditional agricultural city. The study area was mostly LL agglomeration, mainly distributed in the central northern, eastern, and south parts, where the land use types were mostly cultivated land, and grassland. At this time, HH agglomeration was more scattered, mainly concentrated in the central, northern, and other regions; LH and HL types occupy a larger area Small and scattered overall. In 2005, the study area entered a period of rapid development, and the LL agglomeration areas were significantly reduced, mainly concentrated in the east and scattered in the north and southwest parts. At this time, the HH agglomeration areas were further increased and scattered to central and northern parts. The Moran's I index in 2020 is 0.68, which has increased compared with 2005. The proportion of HH types has increased, mainly concentrated in the central and northern and a part of the southern. At the same time, the LL agglomeration areas have decreased further shrink, mainly concentrated in the east and scattered in the northwest part of the study area.

Figure 10. LISA significance map.

It can be seen from Figure 10 that under the ecological security scenario in 2035 and 2050, ecosystem services are generally distributed mainly in HH and LL clusters. The HH types are mainly distributed in the central and northern parts of the study area. The LL agglomeration areas are mainly distributed in the east and scattered in the north and southwest parts; secondly, the HL agglomeration is a more widely distributed district, mainly distributed in the northwest of the study area. It mainly shows that high-value areas of ecosystem services are surrounded by low-value areas, indicating that the ecosystem service capacity of this agglomeration area is higher than that of surrounding areas, and plays a vital role in the supply of ecosystem services in the region; the smallest distribution range is the LH agglomeration area, mainly distributed in the waterbodies of the study area, showing that areas with low value of ecosystem services are replaced by high Surrounded by value areas, it means that the ecosystem service capacity of the agglomeration area is lower than that of surrounding areas, which has a certain inhibitory effect on the improvement of the regional ecosystem service capacity, is not conducive to regional development, and requires focused management and control.

Overall, the global Moran's I index showed a trend of first decreasing and then increasing from 2000 to 2020. Except for insignificant types, the study area is dominated by HH types, which are mainly distributed in the mountainous part with a better natural ecological environment in the north, such areas are rich in natural resources and have a good ecological environment. Soil and water conservation and natural resource protection can effectively improve the ESV of such gathering areas; LL areas have a certain overlap with construction land, and other types of distribution are more scattered and account for a relative proportion. There are few, and focusing on the construction of green infrastructure within the city and the construction of ecological network patterns can effectively improve ESV levels; LH and HL types are scattered near the urban fringe areas of the central city. There is a negative spatial correlation in some areas. There is a certain irrationality in the layout of land use types in such areas. The construction of ecological corridors should be strengthened to maintain the integrity of regional ecological structure-process-function.

3.6. Sensitivity analysis of ecosystem service value

A sensitivity analysis is a useful technique for assessing the impact of changes in input parameters or assumptions on the output of a model or analysis (Chu-Agor et al. 2011; Rahman and Szabó 2021). A sensitivity analysis was calculated to determine the CS value by utilizing Equation (7(7) $\mathbf{I}={\displaystyle \frac{n{\sum }_{i=1}^n{\sum }_{i=1}^n{w}_{\mathrm{ij}}(x_i-x^{-})(x_i-x^{-})}{{\sum }_{i=1}^n{\sum }_{i=1}^n{w}_{\mathrm{ij}}{\sum }_{i=1}^n{(x_i-x^{-})}^2}}$(7) ), intending to evaluate the potential influences of adjusting the service value coefficient for chosen LULC types, the ±50% adjustment in the service value coefficient was based on the recommendation proposed by Costanza et al. (1997; 2014). The CS of the different LULC types in the study area is less than 1 in all cases, the sensitivity of cultivated land and grassland towards total ESV in all indicators was less than 0.2, and the lowest value is 0.05 (Table 14). Generally, the total ESV is inelastic, which proves that the value coefficient selected in this study is more appropriate for the study area. Adjustment to the modified value coefficient for the waterbody had a larger significant impact on the estimated ESV than the value coefficient for other LULC types. When the value coefficient for the waterbody was increased by 50%, the estimated ESV also increased by about 31.7%. The CS values ranged from 0.05 to 0.1 for cropland, to higher than 0.6–0.95 for waterbody when value coefficients for these LULC classes were adjusted by 50%. In contrast, the results of the analysis indicate that while there was a ± 50% change in the grassland value coefficient, the overall change in total ESV was less than ±5.1%. This indicates that the ESV estimates using the ESV coefficients were relatively robust.

Table 14. Estimated total values of ecosystem services and coefficient of sensitivity (CS) after adjusting the ecosystem service value coefficient (VC) by 50%.

Change of VC	1990		2005		2020		2035		2050
	%	CS	%	CS	%	CS	%	CS	%	CS
Cultivated land VC ± 50%	± 2.6	± 0.05	± 3.7	± 0.07	± 4.8	± 0.09	± 6.1	± 0.08	± 7.3	± 0.1
Grassland VC ± 50%	± 5.1	± 0.11	± 4.4	± 0.09	± 3.6	± 0.08	± 2.8	± 0.06	± 2.1	± 0.05
Waterbody VC ± 50%	± 31.7	± 0.6	± 32.6	± 0.7	± 33.7	± 0.9	± 34.6	± 0.8	± 35.4	± 0.95

4. Discussion

4.1. LULC change dynamics

This study provides to measure of LULC changes and their impacts on ecosystem services over the entire study period. Our findings regarding the dynamics of LULC change are in close conformity with those of other studies (Jodder et al. 2023; Rahaman, Kafy, Saha, et al. 2022). Overall, vegetated land and barren land were the major LULC classes and constituted 92% of the total land cover. Urbanization is increasing in the study area, partly due to financial incentives for farmers to build on their land, leading to rapid urban expansion and environmental concerns (Dammag, Jian, et al. 2023; Dey et al. 2021). To achieve Sustainable Development Goal (SDG) 11, which seeks to balance urbanization with environmental protection, the study advocates for strengthened land supervision, legislation to resolve disputes, and the revitalization of land registration systems. Despite past governmental regulations on land use, there's an urgent need for stricter land supervision and sustainable policies, especially given the growing popularity of the khat tree and its effects on groundwater. Population growth, economic shifts, and GDP changes are significant drivers of LULC changes, expected to rise by 4.2%–5% in the coming years, emphasizing the need for careful regional development planning. The study aligns with findings by A. Al Kafy et al. (2022) regarding LULC changes in the tropical catchment of Rajshahi, Bangladesh, where the built-up area is predicted to increase by 9.39 and 13.34 km² in 2029 and 2039, respectively. The city's GDP experienced fluctuations, notably a decline during the Arab Spring revolution but saw a resurgence in the real estate sector recently. The terrain, climate, and overgrazing are other long-term LULC change contributors, suggesting measures like water conservation efforts. Overgrazing and the rise in livestock population have adverse effects on the ecological environment and biodiversity, advocating for ecological resettlement projects (Msofe, Sheng, and Lyimo 2019). Policy factors play a crucial role in regional LULC changes, with short-term changes influenced by socioeconomic factors, and natural factors having a longer-term impact (A. Al Kafy, Faisal, et al. 2022). To ensure sustainable land use, there's a call for increased governmental investment, infrastructure development, and robust policy implementation (Rahaman, Kafy, Faisal, et al. 2022). This study not only addresses the specific characteristics of the study area and the availability of relevant remote sensing and auxiliary data but also incorporates an analysis of SDGs 11, 15, and 17, which could provide insights that can inform sustainable land use planning and policy-making, ensuring the balance between development and conservation, and addressing the specific needs and challenges of the Ibb City community.

4.2. Changes in the total ecosystem service values

Based on remote sensing images and GIS analysis, we conducted a comprehensive study to quantify the changes in ESV in response to LULC change in Ibb city over the entire study period from 1990 to 2050. Our findings observed an increase in certain ecosystem services and a decline in others over the study period was primarily attributed to the patterns and impacts of LULC changes driven by urbanization, economic development, and population growth. The increase in the built-up area led to the conversion of croplands and grasslands, which have high ecological values to urban areas with comparatively lower ecological value. This urban expansion is directly revealed in the decline of provisioning and regulatory services, which are spatially concentrated around urbanized areas. Conversely, certain supporting services showed an increase, due to the relative abundance of water bodies in the study area and their high-value coefficients for individual ecosystem service functions (Liu et al. 2019; Wu 2021). Therefore, preserving natural areas is vital for maintaining ecosystem functions, thus strategies for conserving and restoring cropland, and grasslands are recommended. In light of Ibb City being a key foodshed, the adoption of sustainable agricultural practices, including agroforestry, organic farming, and water-efficient irrigation, is imperative for the enhancement of provisioning services and food production (Morshed et al. 2022; Yuan et al. 2019). Regulatory services, such as water supply and erosion control, are vital in the context of Ibb City’s topography and climate, especially given the challenges posed by urban expansion and climate change. Supporting services like biological control and pollination are significant for maintaining biodiversity and agricultural productivity, which are key to the local economy and food security.

According to our study, the built-up area increased from 0.38% to 4.1% in 2050, at the expense of cropland and grassland. Since the built-up areas have no value for ecosystem service, the total ESV has decreased during the past three decades as a result of urbanization (Chatterjee et al. 2022; Haque et al. 2023). In addition, the impact of climate drought and irrational land use is also greater. If this trend is not effectively controlled, it will have an extremely adverse impact on maintaining ecosystem stability, and the protection and improvement of the ecological environment in the region, as Malaysian (Rahaman, Kafy, Faisal, et al. 2022), Bangladesh (A. A. Kafy et al. 2021), India (Basu et al. 2023). Therefore, implementing climate change adaptation measures, such as developing urban green spaces and water conservation strategies, is crucial. Policymakers must also prioritize biodiversity and ecosystem protection, combating land degradation and biodiversity loss through protected areas and sustainable grazing practices.

In terms of the total value of ecosystem service functions, as LULC intensity has increased, the value of supporting services has shown a consistent year-by-year increase. On the other hand, provisioning services, regulation services, and cultural services have exhibited a downward trend over time. The supporting services contributed approximately 21.8%, 22.9%, 23.6%, 24.3%, and 25.5% of total ESV in the historical years 1990, 2005, and 2020, and the projected years 2035 and 2050. The primary explanation for the significant contribution of supporting services is the abundance of water bodies and their high-value coefficient for each individual ecosystem service function (Grammatikopoulou and Vačkářová 2021). Although the value of functions such as food production, erosion control, and water supply were three major contributors in 1990, the values of erosion control and water supply have shown a decreasing trend over the study period. Specifically, the individual value of erosion control services decreased by 28.1% (US$ 1.6 Mn), while the value of water supply services decreased by 36.7% (US$ 1.8 million) in 2050. In contrast, the individual value of biological control services increased by US$ 3.5 Mn and Pollination by US$2.8 million. When considering the analysis of LULC change, it becomes apparent that differences in the ESV coefficients result in corresponding variations in the impact of LULC change on the estimated ESV. The changing trend of each LULC type aligns with the changing trend of the corresponding ESV (A. Das and Basu 2020; Morshed et al. 2022; F. Zhang, Yushanjiang, and Jing 2019).

The results showed that human activities are the main factors leading to changes in LULC types and ESV, so relevant departments should take corresponding measures to promote sustainable and coordinated regional economic development and strengthen ecological environment protection (Berihun et al. 2021; Chen et al. 2022). In addition, we found there is a trade-off relationship between supply services and other ecosystem services, which is vital for promoting sustainable development in the regional economy, safeguarding the ecological environment, and ensuring the sustainable provision of ecosystem service (Luo et al. 2022), mainly because there is a certain competition between the land used for food and raw material production and other land uses. As a result, reducing LULC types such as cropland and grassland can lead to a significant reduction in total ESV, ultimately resulting in decreased biodiversity, increased land degradation, local warming, and habitat loss (C. Li et al. 2021; Olén et al. 2021). Ibb City has a large area of hills, and its ecosystem is relatively fragile and most sensitive to human disturbance (Dammag et al. 2023). Therefore, to maintain a high level of ESV and enable it to serve humans long-term in a sustainable manner, an integrated approach to land management policies is needed, one that comprehensively considers ecological, social, and economic factors, and addresses the trade-offs between different land uses. This calls for increased governmental investment and infrastructure development, supporting sustainable land use practices. Lastly, establishing a robust system for monitoring and evaluating the effectiveness of these strategies will be key to ensuring their successful implementation and adjustment over time. Engaging local communities and stakeholders in the policymaking process, coupled with educational initiatives, will foster community support and participation. These recommendations aim to strike a balance between developmental needs and ecological conservation, ensuring sustainable resource management and the well-being of the Ibb City community.

4.2. Limitations and future scope of study

This study assessed LULC changes in Ibb's ESV from 1990 to 2050, encountering uncertainties mainly due to data acquisition challenges.(Msofe et al. 2020; Niu et al. 2022). Firstly, the use of BTM for ESV, although practical in data-scarce regions like ours, introduces uncertainties. These arise from applying value coefficients derived from different ecological and socio-economic contexts, which may not accurately reflect the unique characteristics of Ibb City. Therefore, we recommend developing a value coefficient taking into account the local context. For Yemen, there is no specific value coefficient for different types of LULC. So, future studies can focus on developing the value coefficient in light of the local contest of Yemen. Secondly, the CA-Markov model, employed for predicting future LULC changes, inherently assumes that past trends will continue into the future. This assumption may not hold in the face of rapid socio-economic and environmental changes, potentially leading to less accurate future land cover projections (X. Gao et al. 2021; Yang et al. 2020). The static nature of value coefficients used in BTM may not accurately represent the temporal changes in ES in Ibb City. To partially mitigate these limitations, our study incorporated local expert knowledge wherever possible to adjust and validate the transferred values. Nonetheless, we acknowledge that the results obtained through BTM should be interpreted with caution, considering these limitations and uncertainties. The findings provide a valuable, understanding of ESV in Ibb City, highlighting the need for more Future studies should enhance ESV evaluation methods, clarify parameter connotations and values, and conduct long-term field surveys for dynamic parameters like net primary productivity and soil erosion.

5. Conclusions

In this study, our aim was to comprehensively assess the impacts of LULC changes on ESV in Ibb City, Yemen, over a period of three decades (1990–2020) and project these trends into the future (2035 and 2050). The results of the analysis indicated a substantial shift in land cover, characterized by an expansion of built-up areas and barren land, and a corresponding decline in cultivated land and grassland. These changes have led to a notable decrease in the total ESV, from US$ 68.5 million in 1990 to US$ 65.2 million in 2020, with predictions of further decline to US$ 61.2 million by 2050. The sensitivity analysis conducted reinforces the robustness of our ESV estimations and highlights the significant role of specific land cover types, particularly water bodies, in influencing the overall ecosystem service provision. The continuous urban expansion and agricultural land conversion in Ibb City, driven by population pressure and socio-economic activities, pose a critical challenge to the sustainability of regional ecosystems. Our study emphasizes the need for integrated land management strategies that prioritize the protection of high-cover croplands and grasslands, which are crucial for maintaining high ESV and are currently facing rapid degradation. The results underline the importance of incorporating ESV considerations into urban planning and policy-making to balance development needs with ecological preservation. This study contributes valuable insights into the spatial relationship between LULC changes and ESV in a rapidly urbanizing context. It provides a basis for future studies and policy interventions aimed at promoting sustainable land use practices and safeguarding the regional eco-environment in Yemen. The approach and findings of this study are not only relevant to Ibb City but also offer a framework that can be applied to similar data-scarce, rapidly changing urban environments.

Consent to publish

All authors of a manuscript have read and agreed to the published version before it is submitted for publication.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The data presented in this study are available from the corresponding author (B.Q.D.) upon request.

Additional information

Funding

The authors have not received support from any organization for the work submitted.