Ecosystem Health Assessment of the Zerendy District, Kazakhstan

1. Introduction

Despite specific difficulties that have arisen in the schedule of the 2030 Agenda, the Sustainable Development Goals (SDGs) are still within reach, but more efforts are required [1]. One of the areas of effort should be research aimed at identifying the internal patterns of “sustainable development”, which can contribute to the earliest possible achievement of the 2030 SDGs. At the same time, despite the passage of more than three centuries since the concept of “sustainability” emerged and the adoption of 17 SDGs, the discussion on the internal structure of “sustainable development” [2,3,4,5] continues. For example, epistemological studies have shown that concepts such as sustainability, resilience, adaptation, transformation, vulnerability and innovation require further understanding and clarification of their application areas [6,7,8,9,10,11]. Importantly, these concepts provide additional information that serves as a basis for developing new approaches and strategies for implementing the SDGs. They also enable the prediction of future research directions in the field of sustainable development [12].

One concept contributing to achieving SDG 2030 is Ecosystem Health (EH), a product of interdisciplinary research [13]. Its aim is to improve the health and well-being of people, animals, and ecosystems and is partly based on the principle of sustainability [14]. The development of EH dates back to 1941. However, despite the constant enrichment and development of the EH concept, there is still no generally accepted definition [15,16,17,18,19,20]. EH is expected to support the stability and integrity of the ecosystem’s internal structure and function, ensure resilience and recovery under external pressures, and provide normal service functions to meet the needs of human social development [20]. According to [21], EH is a relative concept where the relative health of an ecosystem is measured in a specific spatial location and over time. EH combines natural and social sciences [22,23,24,25]. It is believed [24] that achieving EH should become the cornerstone of sustainable development policy, as maintaining it is the main factor in achieving sustainable socio-economic and environmental development. According to [13], research in the field of EH has gone through the stages of anthropocentrism and biocentrism and is currently based on ecocentrism. EH has gradually become an essential issue of our time, as evidenced by the rapid increase in publications on EH over the past ten years [13,26,27]. Several organizations [28,29,30], including individual divisions of the UN [25], and many researchers from different countries and schools are directly involved in the development of EH.

EH requires continuous determination and assessment of its condition [31]. Therefore, the current focus of research has shifted to the methods of EH Assessment (EHA) [15] and the development of its indicators [31]. Various approaches are used for EHA, which can be divided into two large groups: the biological indicator method and the index system method [15]. The first method is based on the use of individual indicator organisms. Due to the lack of information on such indicators, this paper does not consider the biological indicator method. The second method uses indicators combined with socio-economic data, determines their weights and forms a complex index for EHA. Therefore, the index system method is the focus of this paper.

Index systems have been created to synthesize complex information to reflect the actual situation on EH. The prototype of the index method can be considered Vigor, Organization, and Resilience (VOR) [32]. Based on this idea, researchers developed Vigor, Organization, Resilience, and Stability (VORS) [33] and Vigor, Organization, Resilience, Ecosystem Services, and Stability (VORES) [34] and further enhanced it by constructing optimized and improved subindices [35]. The next most frequently used model is the Pressure-State-Response (PSR) model [17,19,36,37] proposed by the Organisation for Economic Co-operation and Development (OECD) [38]. The PSR model is also subject to development. For example, its modifications such as Driving Forces-Pressure-State-Response (DFSR) [20], Driving Forces-Pressure-State-Impact-Response (DPSIR) [39,40,41], Driving Forces-Pressure-State-Impact-Response-Management (DPSIRM) [42], Driving Forces-Pressure-State-Change-Response (DPSCR) [43], Driving Forces-Pressure-State-Exposure-Effect-Action (DPSEEA) [44] and Basic-Pressure-State-Response (BPSR) [19] have been developed. At present, other EHA models have been created. These are Technique for Order Preference by Similarity to Ideal Solution (TOPSIS) [45] and Backpropagation Neural Network (BPNN) [46], which were improved to Genetic Algorithm-Backpropagation Artificial Neural Network (GA-BPANN) [47]. However, it is difficult to give preference to one of these models since, so far, no researcher has used two or more different methods for EHA in the same area on the same time scale [15]. At the same time, all other things being equal, the VOR and PSR methods are improving most rapidly, but as we see from the above, the development of PSR is more intensive. This served as the basis for using the latest modification of PSR in this study – BPSR, which significantly facilitates the comprehensive assessment of the health of the regional ecosystem [19].

One of the disadvantages of the index system is that it is considered to be ignoring detailed information [15]. It is believed that local or regional EHA should focus on studying its temporal and spatial differences [48]. According to [17], each model of the index system is only a basic organizing paradigm for assessing the state of EH. Because each ecosystem has indicators that reflect various aspects of its health, the EHA example, with EHA, cities can use environmental and production functions as internal and economic indicators as external indicators [50].

After establishing an appropriate index system, the Ecosystem Health Index (EHI) is calculated based on the weight assigned to each metric. The choice of weighting method leads to varying ranges of EHI. As for the definition of the indicators, it has been established that there is no uniform standard of EHA. It is determined based on the understanding and interpretation of the connotation of the concept of “ecosystem health” in combination with the results of previous studies and the specific types of ecosystems selected by the researchers, which is entirely subjective [13]. In this case, researchers usually adhere to the principles of integrity, simplicity, and dynamic response, which are not indicators but macro dimensions. Within the framework of macro dimensions, it is relatively easy to understand the drivers of EH, which are social, economic, and environmental indicators [51].

The object of our study is the Zerendy district, a conglomerate of three critical territories. First, it is a city of oblast significance, which is an enclave of this Area of Interest (AOI). Next are the rural divisions of the district engaged in agro-industrial activities, which are the most minor units for which statistical information is collected on a national scale. In addition, a significant part of the territory of the Zerendy district is occupied by a natural national park and a nature reserve or protected areas. Therefore, the direction of our research on EHA is unique due to the coverage of the territory for different target purposes. (Figure 1).

In EHA, of the three areas mentioned, researchers have paid more attention to cities due to their essential role in economic development and the accompanying deterioration in the quality of the urban environment in the context of rapid urbanization [35,50,52,53,54,55,56,57,58].

The peculiarity of rural EHA is that they include a vast diversity. When studying rural health, they are often divided into disciplinary boundaries, which limits the understanding of the health of rural ecosystems as a whole. The many connections between soils, plants, animals and people related to health challenges all monodisciplinary health studies in agricultural ecosystems. Therefore, an interdisciplinary approach is necessary for a more integrated and comprehensive understanding of rural ecosystem health [59,60]. Such studies on EHA for rural regions have typically covered large regions [17,33,61,62,63,64]. However, we have not found any studies that have examined the EHA of the most minor subunits of administrative-territorial divisions, which are below the rural district – such as rural subdivisions.

Of considerable interest are the EHA of territories occupied by national natural parks [65,66], which are divided into strictly protected zones and zones for economic purposes, where work is carried out that is necessary for the functioning and development of the national natural park, as well as the provision of socio-economic infrastructure [67].

Based on the above, the goal of our research was the EHA of the Zerendy district through the construction of a conceptual model based on the BPSR structure and the research objectives:

- identifying changes in individual LULC classes;

- compiling a set of indicators corresponding to the BPSR model;

- performing EHA of individual administrative-territorial divisions of the district;

- analyzing the spatial dependence of the main clusters;

- Identify the main driving forces influencing the district’s EH.

2. Materials and Methods

2.1. Study Area

Zerendy district is located in the north of the Akmola Region of Kazakhstan (52°33′–53°40′N, 68°28′–70°4′E). It includes 22 rural districts (Figure 2) [68], which occupy the third tier in the republican administrative-territorial structure after the oblast and district [69]. Socioeconomic data are collected on them, which can be used for EHA.

The city of Kokshetau is an enclave surrounded by the territory of the Zerendy district and serves as the Centre of the Akmola oblast [70]. The city of Kokshetau is not administratively part of the Zerendy district. However, in this study, the Zerendy district and the city of Kokshetau are considered a single space for assessing the EH since the impact of the city on rural areas is inevitable. The territory of the Zerendy district is located on the western edge of the Kazakh folded country [73], where a relief is characterized by a low-mountain lowland plain. The hydrographic network of the district includes the rivers Shagalaly and its right tributary, Kylshykty, which flow from south to north. In the south of the district, the rivers Zhabay, Arshaly, and Koshkarbay originate from the Ishim River basin. There are also many small lakes in the district territory [71]. The district’s climate is continental, with long, low-snow winters and dry, warm summers. The average temperature in January is (−18°C), and in July (+19°C). The average annual precipitation varies from 350 to 400 mm. Snow cover lasts 140-160 days, but its average height is only 20 cm [74]. The Zerendy district and its enclave uniquely combine urbanized, agricultural and nature conservation areas. The city of Kokshetau, as well as the Kokshetau National Park and the Zerendy Zoological Reserve (Figure 3), have a significant impact on the ecosystems of the district.

If the city’s anthropogenic impact is a pressing concern, the national park and the reserve, in stark contrast, are actively working towards the preservation and restoration of the district’s ecosystem. In general, the heterogeneity of the natural and socio-economic conditions of development in different parts of the district makes it a unique object for EHA.

2.2. Methods

For EHA, we adopted BPSR [19] as the base method. The study was divided into four stages, shown in Figure 4.

2.2.1. Data Collection

Data collection for this study was conducted in two directions. The first direction was spatiotemporal data, including meteorological data classified by LULC classes from Landsat 7/8 images and landscape metrics. The second direction was socioeconomic statistical information by administrative-territorial divisions.

2.2.2. Data Processing

Satellite data from the Earth Engine Data Catalog [75] were used to obtain the spatiotemporal data. Some of their characteristics are given in Table 1.

2.2.2.1. Land Use and Land Cover Classification

Landsat 7/8 images were used to analyze the LULC changes from 2010, 2016 and 2023. The LULC was classified into five classes:

1)

>1) Cropland - CLD,

2)

>2) Pasture - PTE,

3)

>3) Forest - FET,

4)

>4) Water Bodies -WBS,

5)

>5) Urban Land – ULD

using the RF algorithm on the Google Earth Engine (GEE) platform. To improve the classification accuracy were used the following indices: NDVI, EVI, NDWI, NDTI, NDMI, MNDWI [76,77,78,79,80,81], as well as NLT [82], which are described in detail in our previous studies [83,84].

2.2.2.2. Landscape Metrics

Data related to landscape change measures such as Shannon’s diversity index, Contagion index (%), Shannon’s evenness index, Landscape Division Index, Interspersion Juxtaposition Index, and FDIAWM were calculated in Fragstats 4.2 software [85].

2.2.2.3. Statistical Information

The statistical data on administrative-territorial divisions that we used in our study were obtained from the Yearbook of the Bureau of National Statistics Agency for Strategic Planning and Reforms of the Republic of Kazakhstan [86].

2.2.3. Ecosystem Health Assessment

To calculate the EH level, the BPSR system [19] was used, which organizes indicators into four different but interrelated categories: B (Basic), P (Pressure), S (State) and R (Response). In this study, B is the degree of impact of climatic factors, P is the level of pressure exerted by anthropogenic and natural factors on the study area, S characterizes the state of the ecosystem, and R describes the response to pressure. All BPSR indicators we collected and their weighting coefficients are presented in Table 2.

Ecosystem resilience (Rs) was assessed based on the area occupied by LULC classes and calculated using the following formula [87]:

$$\mathit{R}\mathit{s}={\sum }_{\mathit{k}=1}^{\mathit{m}}{\mathit{P}}_{\mathit{k}}\times {\mathit{E}\mathit{C}}_{\mathit{k}}$$

(1)

where ${\mathit{P}}_{\mathit{k}}$_k is the area of the class in the district, and ${\mathit{E}\mathit{C}}_{\mathit{k}}$ is the ecosystem elasticity coefficient of the given class.

Based on studies [88], the ecosystem Elasticity coefficient (EC) of the classes was measured using the following equation:

EC=0.7×Resist+0.3×Resil

(2)

where EC denotes the ecosystem elasticity coefficient; ’Resil’ refers to the Resilience coefficient of the LULC class, and ’Resist’ refers to the Resistance coefficient of the LULC class.

In this study, more emphasis was placed on Resistance, as the area occupied by Cropland expanded rapidly, increasing interventions and disturbances. Therefore, greater weight was assigned to Resistance (0.7) and less to Resilience (0.3). The coefficient applied for each land-use type is listed in Table 3.

Given the differences in measurement and scale of the selected indicators, the data were standardized [89] before analysis using the following equations:

For positive indicators:

$${\mathit{X}}_{\mathit{i}\mathit{j}}={\displaystyle \frac{{\mathbf{x}}_{\mathbf{i}\mathbf{j}}-\mathbf{m}\mathbf{i}\mathbf{n}\left({\mathbf{x}}_{\mathbf{j}}\right)}{\mathbf{m}\mathbf{a}\mathbf{x}\left({\mathbf{x}}_{\mathbf{j}}\right)-\mathbf{m}\mathbf{i}\mathbf{n}\left({\mathbf{x}}_{\mathbf{j}}\right)}}$$

(3)

For negative indicators:

$$X_{ij}={\displaystyle \frac{\max \left(x_j\right)-x_{ij}}{\mathrm{m}\mathrm{a}\mathrm{x}\left(x_j\right)-\mathrm{m}\mathrm{i}\mathrm{n}\left(x_j\right)}}$$

(4)

where ${\mathit{X}}_{\mathit{i}\mathit{j}}$ is the normalized value of indicator j in year i, ${\mathit{x}}_{\mathit{i}\mathit{j}}$represents the value of indicator j in year i; $\mathbf{max}\left({\mathit{x}}_{\mathit{j}}\right)$and $\mathbf{m}\mathbf{i}\mathbf{n}\left({\mathit{x}}_{\mathit{j}}\right)$, denote the maximum and minimum values of indicator j overall years, respectively. All index values range from 0 to 1.

2.2.3.1. Determination of Indicator Weights

This instrument uses the entropy method to determine weighting factors, which measure the uncertainty or variability of data. The greater the variability indicator, the greater its weight. The entropy method allows for an objective consideration of the originality of data and the depreciation of outlier estimates. In addition, this method can correct for limitations associated with the relationships between indicators, making it more efficient and practical [90]. Calculating weighting factors using the entropy method consists of the following steps:

1) Calculation of the share of indicator j in year i:

$${\mathit{P}}_{\mathit{i}\mathit{j}}={\displaystyle \frac{{\mathit{X}}_{\mathit{i}\mathit{j}}}{{\sum }_{\mathit{i}=1}^{\mathit{n}}{\mathit{X}}_{\mathit{i}\mathit{j}}}},(\mathit{i}=1,2,\dots ,\mathit{n};\mathit{j}=1,2,\dots ,\mathit{m})$$

(5)

2) Calculating the entropy level of the index j:

$${\mathit{I}\mathit{E}}_{\mathit{j}}=-\mathit{k}\sum _{\mathit{i}=1}^{\mathit{n}}{\mathit{P}}_{\mathit{i}\mathit{j}}\mathbf{ln}\left({\mathit{P}}_{\mathit{i}\mathit{j}}\right),\mathit{k}={\displaystyle \frac{1}{\mathbf{l}\mathbf{n}\left(\mathit{n}\right)}},{\mathit{I}\mathit{E}}_{\mathit{j}}\ge 0$$

(6)

3) Calculation of the difference coefficient for the indicator j:

$${\mathit{C}}_{\mathit{j}}=1-{\mathit{I}\mathit{E}}_{\mathit{j}}$$

(7)

4) Calculation of the weight of the indicator j:

$${\mathit{W}}_{\mathit{j}}={\displaystyle \frac{{\mathit{C}}_{\mathit{j}}}{{\sum }_{\mathit{j}=1}^{\mathit{m}}{\mathit{C}}_{\mathit{j}}}},\left(\mathit{j}=1,2,\dots ,\mathit{m}\right),\sum _{\mathit{j}=1}^{\mathit{m}}{\mathit{W}}_{\mathit{j}}=1$$

(8)

where ${\mathit{P}}_{\mathit{i}\mathit{j}}$the share of indicator j in the year i, ${\mathit{X}}_{\mathit{i}\mathit{j}}$is the value of indicator j in the year i, ${\mathit{I}\mathit{E}}_{\mathit{j}}$ is the entropy of indicator j, ${\mathit{C}}_{\mathit{j}}$is the difference coefficient of indicator j, ${\mathit{W}}_{\mathit{j}}$is the weight coefficient of the indicator.

2.2.3.2. Calculation of EH Using of Indicator Weights

The calculation of the ecosystem health level was carried out according to the following formula [91]:

$${ \mathit{E}\mathit{H}}_{\mathit{i}}=\sum _{\bullet }^{\mathit{n}}{\mathit{W}}_{\mathit{j}}\times {\mathit{X}}_{\mathit{i}\mathit{j}}(\mathit{j}=\mathit{B}1,\mathit{B}2,\mathit{B}3,\widehat{\mathit{A}}\bullet \widehat{\mathit{A}}\bullet \widehat{\mathit{A}}\bullet ,\mathit{R}9)$$

(9)

where ${\mathit{E}\mathit{H}}_{\mathit{i}}$ is the ecosystem health level for year i, ${\mathit{W}}_{\mathit{j}}$ is the weighting coefficient of category j, ${\mathit{X}}_{\mathit{i}\mathit{j}}$is the normalized value of indicator j for year i.

In turn, the values of the BPSR categories are calculated using the following formulas [19]:

$${\mathit{B}}_{\mathit{i}}=\sum _{\mathit{j}=\mathit{B}1}^{\mathit{n}}{\mathit{W}}_{\mathit{j}}\times {\mathit{X}}_{\mathit{i}\mathit{j}}(\mathit{j}=\mathit{B}1,\mathit{B}2,\mathit{B}3,\mathit{B}4)$$

(10)

$${ \mathit{P}}_{\mathit{i}}=\sum _{\mathit{j}=\mathit{P}1}^{\mathit{n}}{\mathit{W}}_{\mathit{j}}\times {\mathit{X}}_{\mathit{i}\mathit{j}}(\mathit{j}=\mathit{P}1,\mathit{P}2,\mathit{P}3,\widehat{\mathit{A}}\bullet \widehat{\mathit{A}}\bullet \widehat{\mathit{A}}\bullet ,\mathit{P}7)$$

(11)

$${ \mathit{S}}_{\mathit{i}}=\sum _{\mathit{j}=\mathit{S}1}^{\mathit{n}}{\mathit{W}}_{\mathit{j}}\times {\mathit{X}}_{\mathit{i}\mathit{j}}(\mathit{j}=\mathit{S}1,\mathit{S}2,\mathit{S}3,\widehat{\mathit{A}}\bullet \widehat{\mathit{A}}\bullet \widehat{\mathit{A}}\bullet ,\mathit{S}10)$$

(12)

$${ \mathit{R}}_{\mathit{i}}=\sum _{\mathit{j}=\mathit{R}1}^{\mathit{n}}{\mathit{W}}_{\mathit{j}}\times {\mathit{X}}_{\mathit{i}\mathit{j}}(\mathit{j}=\mathit{R}1,\mathit{R}2,\mathit{R}3,\widehat{\mathit{A}}\bullet \widehat{\mathit{A}}\bullet \widehat{\mathit{A}}\bullet ,\mathit{R}9)$$

(13)

where ${\mathit{B}}_{\mathit{i}}$, ${\mathit{P}}_{\mathit{i}}$, ${\mathit{S}}_{\mathit{i}}$, ${\mathit{R}}_{\mathit{i}}$are the category values taking into account the weighting coefficient for year i, ${\mathit{W}}_{\mathit{j}}$is the weighting coefficient of indicator j, ${\mathit{X}}_{\mathit{i}\mathit{j}}$is the normalized value of indicator j for year i.

2.2.4. Analysis Methods and Spatial Correlation Indicators for EH

Moran’s I Index

The Moran’s I Index is used in this study to identify the characteristics of spatial autocorrelation. It allows us to determine whether there are spatial patterns in the data, such as clustering of high or low values of the indicators, and whether one value influences neighbouring values [92].

In calculating Moran’s I index, spatial relationships were defined based on (CONTIGUITY_EDGES_ONLY) only neighbouring polygon features that share an adjacent boundary or overlap. The global Moran’s I index and LISA calculations were performed in the ArcMap software environment version 10.8.

The formula for calculating the global Moran’s I index [93]

$$\mathit{I}={\displaystyle \frac{\mathit{n}}{{\mathit{S}}_0}}{\displaystyle \frac{\sum _{\mathit{i}=1}^{\mathit{n}}\sum _{\mathit{j}=1}^{\mathit{n}}{\mathit{\omega }}_{\mathit{i},\mathit{j}}{\mathit{z}}_{\mathit{i}}{\mathit{z}}_{\mathit{j}}}{\sum _{\mathit{i}=1}^{\mathit{n}}{\mathit{z}}_{\mathit{i}}^2}}$$

(14)

where ${\mathit{z}}_{\mathit{i}}$ is the deviation of the attribute for object i from its mean value (xi –$\overline{\mathit{X}}$), ${\mathit{\omega }}_{\mathit{i},\mathit{j}}\mathbf{i}$s the spatial weight between objects i and j, n is equal to the total number of objects, S0 is the set of all spatial weights, which in turn is equal to:

$${\mathit{S}}_0=\sum _{\mathit{i}=1}^{\mathit{n}}\sum _{\mathit{j}=1}^{\mathit{n}}{\mathit{\omega }}_{\mathit{i},\mathit{j}}$$

(15)

To analyze clusters, outliers and their distribution, Local indicators of spatial association (LISA) were used using ArcMap version 10.8. The formula for LISA [94]:

$${\mathit{I}}_{\mathit{i}}={\displaystyle \frac{{\mathit{x}}_{\mathit{i}}-\overline{\mathit{X}}}{{\mathit{S}}_{\mathit{i}}^2}}\sum _{\mathit{j}=1,\mathit{j}\ne \mathit{i}}^{\mathit{n}}{\mathit{\omega }}_{\mathit{i},\mathit{j}}({\mathit{x}}_{\mathit{j}}-\overline{\mathit{X}})$$

(16)

where ${\mathit{x}}_{\mathit{i}}$is a numeric attribute of object $\mathit{i}$, $\overline{\mathit{X}}$ is the mean value for this attribute, ${\mathit{\omega }}_{\mathit{i},\mathit{j}}$is the spatial weight for the pair i and j, and:

$${\mathit{S}}_{\mathit{i}}^2={\displaystyle \frac{\sum _{\mathit{j}=1,\mathit{j}\ne \mathit{i}}^{\mathit{n}}{({\mathit{x}}_{\mathit{j}}-\overline{\mathit{X}})}^2}{\mathit{n}-1}}$$

(17)

where n corresponds to the total number of objects.

Application of Principal Component Analysis

Principal component analysis (PCA) was used to identify the main factors influencing the health of ecosystems. This method allowed us to reduce the number of original variables to two key components that best reflect the variability of the data [95]. The analysis included 11 variables: population density, agricultural land area per capita, built-up area, pollutant emissions, cultivated green land area, NDVI, SAVI, NDMI, EVI, NDWI indices and waste disposal costs. PCA was performed using IBM SPSS software.

Based on the identified components, a regression model was constructed to predict the EH level. This approach takes into account the influence of independent variables, which ensures accurate modelling and understanding of how various factors affect EH.

The regression equation is as follows:

$$\mathit{E}\mathit{H}=\mathit{\delta }+{\mathit{\alpha }}_1{\mathit{F}}_1+{\mathit{\alpha }}_2{\mathit{F}}_2+\dots +{\mathit{\alpha }}_{\mathit{n}}{\mathit{F}}_{\mathit{n}}$$

(18)

where ${\mathit{\alpha }}_1,{\mathit{\alpha }}_2,\dots ,{\mathit{\alpha }}_{\mathit{n}}$_nare the coefficients that determine the contribution of each component, ${\mathit{F}}_1,{\mathit{F}}_2,\dots ,{\mathit{F}}_{\mathit{n}}$_nare the principal components obtained as a result of PCA, and δ is a constant term of the equation. The regression analysis was performed using IBM SPSS software [96].

3. Results

3.1. Analysis of Land Use and Land Cover Dynamics in the Zerendy District

Significant changes in LULC have occurred in the study area, which can be easily seen from Figure 5 and Table 4.

By 2023, the dominant class in the Zerendy district becomes CLD, replacing PTE and occupying 50.62% of the district area. The overall growth of the CLD class was 27.72% of the total district area from 2010 (22.90%) to 2023 (50.62%).

The leading class in terms of total area change is PTE. The PTE class has been losing areas throughout the period. It has the most significant percentage of lost area, 32.41% of the total area of the district, which cannot but affect the state of the ecosystem.

The FET class has growth dynamics and has increased its area from 2010 (10.57%) to 14.19% per cent in 2023. Thus, the total growth of the class from the initial area was 34.24% for the entire period.

WBS and ULD are the most minor classes in the area. WBS has shown a slight increase, which may be due to seasonal changes. Despite its small area, the ULD class is also characterized by a constant increase. By the end of the observation, its increase from the initial area reached 59.57%. More clearly, the trends in the change of LULC classes are shown in Figure 6a and Figure 6b.

3.2. Spatiotemporal Changes in Ecosystem Health in the Zerendy District

AOI is characterized by noticeable changes in the level of ecosystem health over the period under consideration (Figure 7).

The results of the calculated EH are presented in Table 5.

We have identified 6 EH levels for the Zerendy district: very low level (1-VL), low level (2-L), medium-low (3-ML), medium (4-M), medium-high (5-MH), high (6-H). In 2010, most counties were at a relatively good level of EH: approximately 70% of counties reached 3-ML and above. In particular, Kokshetau had a level of 1-VL (0.37). The counties of Aidabol (0.59) and Sarozek (0.60) were at 5-MH, while Baiterek (0.69) reached 6-H. By 2016, some areas had seen improvements in the ecosystem. For example, in Prirechnoye, the EH level increased from 0.52 to 0.63, and in Sadovoye, from 0.47 to 0.59. However, despite these improvements, Kokshetau’s EH remained at 1-VL (0.38), and in Konysbay, it began to decrease from 0.47 to 0.42. In 2023, the EH condition worsened. The number of 2-L divisions increased. For example, these are Kyzylsaian (0.55) and Kanai Bi (0.55), which moved from 3-ML to 2-L. The average EH values over the entire observation period showed slight fluctuations. In 2010, it was 0.533; in 2016, it increased slightly to 0.542; and in 2023, it decreased again to 0.530. This indicates a general trend towards deterioration of the ecosystem condition in the AOI despite temporary improvements in individual counties.

Figure 8 illustrates these changes, showing the redistribution of EH across divisions for 2010, 2016, and 2023.

The change in the EH level during the period under consideration had different characteristics in the divisions of the Zerendy district. The center of the district, including Kokshetau, is characterized by a low EH level (from 0.35 to 0.40), while the divisions surrounding Baiterek, which are located within the boundaries of the Zerendy reserve, demonstrate a high EH level (more than 0.55).

3.3. Spatial Dependence of Ecosystem Health of Zerendy District

Quantitative measurements of spatial autocorrelation of the ecosystem health level of the Zerendy district using the global Moran’s I index showed that its values in the studied years (2010, 2016, 2023) ranged from 0.42 to 0.59 (Table 6). This indicates a significant spatial autocorrelation of the ecosystem health level in the AOI. That is, the ecosystem health level in each division is under the positive or negative influence of neighbouring territories.

Regarding time dynamics, the global Moran’s I index decreased towards the middle of the period and increased towards the end. Its value decreased from 2010 to 2016, indicating a reduction in the spatial dependence of ecosystem health and a tendency towards disaggregation. From 2016 to 2023, the index value increased significantly, reflecting a relatively high spatial aggregation of ecosystem health. This may be due to the effectiveness of the district development strategy adopted after the COVID-19 recession. The Z-score is significantly higher than zero in all years, indicating the presence of spatial autocorrelation in each year. At the same time, the P-values are extremely low in all cases (less than 0.01), confirming the statistical significance of the results and the low probability of random distribution. The variance remains relatively stable throughout the years, indicating the reliability of the measurements.

Figure 9 shows the results of using LISA to study local spatial autocorrelation and analyze spatial aggregation existing in the divisions of the Zerendy district. Local autocorrelation is primarily represented by clusters of high values, designated as High-High (HH), and clusters of low values, defined as Low-Low (LL). These types of clusters indicate clusters where high and low EH values tend to concentrate in certain territories. At the same time, such types as High-Low (HL) and Low-High (LH), which are characterized by the presence of abnormal values (low values among high or high among low), are not found in the AOI. This indicates a low degree of mixed influences between high and low values of the level of ecosystem health near each other.

The HH type is the primary agglomeration type in the southwestern part of the district, indicating a concentration of divisions with high levels of ecosystem health. The prominent clustering is observed around the Zerendy Reserve and the Baiterek division. This clustering has been increasing over time. Она, especially noticeable in 2023. The LL type is mainly represented in the central and eastern parts of the region, indicating low values of ecosystem health in these regions. Districts with low EH values surround Kokshetau, the oblast’s center. In 2010, clustering in the central part was less pronounced, but by 2023 it had significantly increased. Thus, the LISA analysis shows that in the Zerendy district, there are clear groups with similar values of ecosystem health. This means that areas with high and low values are located next to each other. This is especially noticeable in the central part of the district, where the values are low, and in the southwestern part, where the values are high.

3.4. Identifying the Role of Individual Factors in Ecosystem Health of Zerendy District Based on Principal Component Analysis

PCA analysis concentrated information on changes in eleven key indicators into two principal components (Table 7).

The first principal component (F1) is a characteristic of economic development, mainly reflecting changes in population density, share of built-up area, emissions of pollutants, and pasture area. The second principal component (F2) is a characteristic of climate change, mainly reflecting vegetation indices such as NDVI, SAVI, NDMI, EVI, and NDWI.

A model of statistical relationships between changes in EH levels and principal components was obtained using regression analysis. The results of the regression model are presented in Table 8.

The R2 value shows that the model explains 72.4% of the variability in the EH level based on the principal components. This indicates a high degree of explanatory power of the model. The significance of the model (p < 0.001) also confirms the statistical significance of the obtained results. The EH level showed a negative relationship with the economic development component (F1) and a positive relationship with the climate change component (F2). This indicates that economic development hurts ecosystem health, while climate change, reflected in vegetation indicators, has a positive effect.

4. Discussion

4.1. Selection of Ecosystem Health Assessment Model

The EHA methodology began with developing the Ecosystem Distress Syndrome [97] and the Environmental Degradation Syndrome [98] or EDS models. These models describe in some detail the influence of internal natural factors on an ecosystem without assessing the influence of external factors, especially those emanating from human activities. Since human activity is the main factor threatening ecosystems [99], EHA models were further developed that consider anthropogenic activity – VOR, PSR and others. The VOR model and its modifications are more advanced than the EDS model since EHA is implemented on the vigour, organisation, and resilience framework. At the same time, the VOR, in its initial version, ignored the driving force of the current state and the ecosystem’s response to this driving force. The PSR model captures the interaction between the environment and anthropogenic activities. PSR views human activities as playing a central role in determining the state of the environment and better reflects the essence of EH [100]. At the same time, the shortcomings of both VOR and PSR models are under the close attention of researchers who modify them to improve EHA [13,19,27,100,101]. For example, the PSR structure included base-level indices, which are climatic conditions and soil properties, while maintaining the index of three main categories (Pressure, State, and Response). The effectiveness of this modification was tested on the territory of alpine pasture systems [19]. To more thoroughly assess the internal connections, mutual adaptability, and coordination between interrelated indicators in the PSR model, a modified version has also been put forward, called the Coupling Coordination model (CC–PSR) [100]. However, the CC–PSR model has only been used for water basin EHAs. The PSR model and its modifications have been successfully applied to urban ecosystem EHAs [103,104], nature conservation areas [65,105], rural ecosystems [106,107], and also to entire province ecosystems [108], which include urban, rural, and nature conservation areas simultaneously, as in the case of the Zerendy district. According to [19], BPSR claims to maximally reflect spatial and structural changes in regional ecosystems and the external influences they experience. Therefore, considering the above advantages, BPSR was adopted as a research model for the EHA Zerendy district, which includes urban, rural and nature conservation areas.

4.2. Selection of Index System and Its Comparison with Previous Studies

Creating a comprehensive EHA index system is the key to its most accurate calculation. At present, there is no consensus on the EHA index system, but there is a set of principles that researchers should follow. The proposed index system should first be representative and accessible, then objectively reflect the specifics of the studied ecosystem and be quantifiable [13,19,27,109]. The number of indices in the PSR model is a variable value that depends on the goals and objectives set for the researchers. For example, in one of the early works [100], only 11 indices were used for EHA using the PSR model. In the most recent publication [101], the number of indices in the PSR model was 16, and for researchers using the BPSR model [19], the total number of indices reached 25. In our case, the number of available indices reflecting the specifics of the Zerendy district ecosystem and subject to quantitative determination was 35. This set of indices claims to be a sufficiently representative index system for an objective EHA of the Zerendy district. At the same time, we do not oppose to the view that each ecosystem has indicators reflecting various health aspects [17]. Therefore, another set of indices may be selected for another ecosystem, which is sufficient for EHA.

The BPSR structure we chose included spatiotemporal data and socioeconomic information related to the Zerendy district ecosystem. The base (B) indicators (4) characterise the natural material basis of the ecosystem. In this study, it is the climatic condition of the study area: average annual temperature, average temperature during the growing season, average annual precipitation and average precipitation during the growing season. The pressure (P) indicators (10) consider human activity’s direct or indirect impact on ecosystems. The indicators of this category of indices are Population density, Share of cropland per capita, Share of pasture land per capita, Share of water resources per capita, Number of cattle, Built-up area, NDBI, Amount of atmospheric emissions, Investments in fixed assets and Volume of industrial services produced. The state (S) indicators (10) reflect the internal structural components of the ecosystem necessary for the creation of an ecological environment: Area of cropland, pasture land, NDBI, NDVI, NDTI, SAVI, MNDWI, NDMI, EVI, NDWI and Ecosystem Resilience. The response (R) indicators (11) are related to the ecosystem’s response to stress, changes in functional structure and measures humans take to restore ecosystem functions: Shannon’s diversity index, Contagion index, Shannon’s evenness index, Landscape Division Index, Interspersion Juxtaposition Index, FDIAWM, Forest Area, Funding for environmental protection and Funding for waste disposal. In other studies, the quantitative composition of the first category (B, P, S, R) was usually equal to (B) or lower (P), (S), (R) than in our work [19,101,102,103,104,105,106,107,108,109,110]. Therefore, it can be argued that the set of indices used in this article as a whole for each category (B), (P), (S), (R) are more complete since they maximally cover the spatiotemporal natural, social and economic aspects of the ecosystem of the Zerendy district.

4.3. Spatial and Temporal Changes in Ecosystem Health in the Zerendy District

Assessing spatiotemporal changes in EH is essential for thoroughly understanding ecosystem changes over time. This study used spatial autocorrelation to analyze such changes, including the global Moran’s I index and LISA. These methods have become widely used in environmental health research [88,91] due to their ability to identify spatial patterns and clusters, allowing for considering relationships between adjacent areas.

The global Moran’s I index is widely used to assess spatial autocorrelation in ecosystems, which is confirmed by the results of studies [111,112]. The Moran’s I value determines the degree of EH connection between individual territorial units. This allows us to identify clearly defined zones of concentration of high or low EH values, which is especially important for regions with a diverse land use structure, such as the Zerendy district. The results obtained using Moran’s I index in this study show that the Zerendy district has a pronounced spatial autocorrelation of the EH level.

The LISA method allows the refinement of the global analysis results by identifying local anomalies in the EH distribution, such as clusters of high or low values, which may be due to unique local conditions [111,112]. The use of LISA for the EH analysis in the Zerendy district revealed "High-High" and "Low-Low" clusters. This confirms that the choice of these methods for the analysis of spatial dependence in this study is justified and consistent with best practices in the field of EH studies.

In addition to spatial analysis methods, driver factors that determine the distribution and changes in ecosystem health play an important role. Among such factors, anthropogenic impacts stand out, including population density, degree of development, land use, and industrial development [113,114,115]. These factors directly and indirectly impact the state of ecosystems, changing their ability to recover and adapt to external influences. This study analyzed the impact of anthropogenic factors on the EH level using regression models and the PCA method. The identified factors allowed for a deeper understanding of which aspects of human activity are most critical for the ecosystem’s conservation and restoration, emphasizing the importance of the selected methods for understanding the processes occurring in the Zerendy district ecosystem.

Thus, this study’s results confirm that using the global Moran’s I index and LISA in combination with PCA is a reasonable choice for analyzing spatiotemporal changes in ecosystem health. The combined application of these methods provides a comprehensive understanding of the spatial structure of EH. The consequence of this approach is more accurate decision-making for improving ecosystem health in the context of achieving sustainable development goals.

4.4. Policy Implications

Applying the BPSR methodology and spatial analysis allows for the identification of the most vulnerable areas requiring priority protection and restoration. In particular, attention should be paid to areas with low EH values, where increased measures to protect the environment and reduce anthropogenic impact are required. The results of the presented study can serve as a basis for developing ecosystem management policies in the Zerendy district. The information obtained can be used by local authorities to improve land use planning and establish sustainable development strategies to preserve the district’s ecosystem in the long term. Interested managers can use the results of our assessment to take adequate measures to improve the health of the Zerendy district ecosystem

4.5. Limitations and Future Work

Despite successfully applying the chosen methodology, the study has several limitations. First, the EH assessment was conducted based on available data only, which could affect the accuracy of some indicators. Second, the spatial analysis was limited to time intervals (2010, 2016, 2023), which may only partially reflect the character of changes in the interperiods. In the future, it is recommended that EH be monitored more frequently to improve the accuracy of the analysis. It is also worth considering the possibility of using data with a higher spatial resolution, which will allow for detailed local-level changes and account for specific factors affecting the ecosystem.

5. Conclusions

Changes in LULC, namely the reduction of pasture area and the increase of arable land and urbanised areas, significantly impact the level of EH Zerendy district. As a result, a precise territorial distribution of EH is observed: the southwestern areas located near the national nature park and nature reserve demonstrate a high level of EH, while the central areas with a high degree of urbanisation have a low level of EH. Spatial dependence analysis revealed two main clusters: high EH values are concentrated in the southwestern regions, while low values are found in the central areas. This is confirmed by the global Moran’s I index and LISA, which revealed significant spatial autocorrelation, indicating a stable clustering of EH values. The results of using PCA and regression analysis showed that both economic and climatic factors have a significant impact on the level of EH Zerendy district. The first principal component is represented by economic factors (population density, share of built-up area, emissions, etc.), which have a negative impact on the EH indicator. The second principal component includes climatic factors and the state of vegetation (NDVI, SAVI, etc.), which, on the contrary, have a positive impact on the level of EH.