Modeling the Influence of Changes in the Edaphic Environment on the Ecosystem Valuation of the Zone of Influence of the Ozogoche and Atillo Lake Systems in Ecuador

: Ecosystem valuation (EV) of soil resources is essential for understanding changes in environmental services in monetary terms. A lack of this information, which includes economic indices, hinders the optimal management of natural resources. This study evaluated the influence of changes in the edaphic ecosystem on the EV of the zone of influence of the Ozogoche and Atillo lake systems in Ecuador. The classification was carried out through spectral indices and support vector machines (SVMs), and the EV was determined through opportunity costs including environmental service provisioning and indirect use. The land use and EV classification methods were performed efficiently; the degradation trend was constant. The Modified Water Difference Index was the most efficient in the extraction of water bodies, with an accuracy of 91%. The SVMs algorithm, in recognizing coverage in general, had an overall accuracy of 85%. The adjustment made to the SVMs algorithm to improve the selection of hyperparameters was effective; a robust architecture of the algorithm in terms of automation was achieved. Between 2000 and 2020, moorland, water and wetland degraded by 19%, 2% and 3.4%, respectively. In 2000, the EV as a function of avoided CO 2 content was USD 8.00 × 10 6 ; in 2010 and 2020, it was USD 6.00 × 10 6 and USD 5.00 × 10 6 , respectively.


Introduction
The high Andean systems are of volcanic origin, and it is believed that more than 50 volcanoes in the Andes of Ecuador gave rise to their development; they are located at an altitude between 3200 and 3600 masl and are home to deep, cold lagoons with a characteristic dark blue color [1,2].The montane cloud forests, shrubby moorlands and herbaceous moorlands are characteristic ecosystems of their environment [3].They are home to 99 plant species with a high level of endemism; among the genera with the greatest diversity are Lysipomia, Centrophogon, Brachyotum, Brachyotum campii, Chusquea falcata, Chuquiraga jussieui, Oritrophium peruviamun, Podocarpus oleifolius, Oreopanax rosei, Brachyotum azuayense and Miconia bullata, and their diversity is attributed to the fact that these taxa probably developed at this site [4].Species with higher endemism are found in areas with constant rainfall, low temperature levels and high relative humidity [5].
The extreme climatic conditions of the Andean areas are attributed to changes in temperature, relative humidity, atmospheric pressure, air density and solar radiation.Annual rainfall values in the study area have great variability throughout the year.The dry weather reaches its highest level in the months of December and January.During this time, sunshine is high and temperature changes are very varied.Humidity tends to decline with temperature variability; consequently, areas with higher humidity tend to have less thermal fluctuation.In the first 30 cm of soil depth, water occupies 61% of its volume.It is estimated that soil weight can double due to the effect of water retention [1].
The importance of environmental services (ESs) for Andean ecosystems lies mainly in water regulation, carbon storage and species habitat [6].Andean areas are considered essential systems for maintaining the balance of microclimates and reducing greenhouse gas emissions [7].The Andean systems as geographical units have been affected by environmental factors and human activities; the increase in temperature and agricultural and livestock activities have caused land-cover degradation [8].On the other hand, the increase in population has increased the demand for food, causing the expansion of the agricultural sector and, consequently, the levels of water consumption for the population and production [9].
ES degradation due to land-use changes has manifested in several of South America's territories.For example, the high Andean plateau of Bombón in Peru is facing an increasing loss of soils due to accelerated anthropogenic land use and climate change.The scarcity of information on the historical development of degradation complicates the development of protection plans for vulnerable areas [10].In the Calchaquíes summits in Argentina, ESs have been affected by climatic conditions, vegetation variations and water balance [11].In the mountainous areas of Colombia, in recent years, ESs have been intensively degraded due to soil degradation, affecting water regulation.A previous study stated that this is a consequence of climate change processes, mining, agriculture and livestock [12].
There are direct and indirect effects on Andean soil uses related to cattle ranching activities.The direct effects are associated with the degradation of the flora of the area, and the indirect effects are generated by the compaction of the soil layer due to the weight of the cattle; the presence of heavy animals in the natural ecosystem causes problems in the porosity of the cover, which hinders the hydraulic circulation through the system.Regarding agriculture, the impact on Andean soils is usually linked to the release of nutrients immobilized in the soil cover.After the first harvest, the soil decreases its nutrient concentration and its structure becomes compact and dry, which causes runoff to flow, dragging with them soil particles [2].
The identification of the factors involved in land changes is important to develop real land-use change models that allow for the analysis of changes in phenology, the detection of vulnerable areas and the analysis of the state of natural resources [3].Methodologies to detect different land uses are closely linked to remote sensing, geographic information systems and computer processes that allow the processing of data sets.There are different modeling approaches to determine land-use changes, for example, dynamic and static models.Dynamic models are related to temporal and dynamic changes in land use, relating changes in past transitions and their trajectories.Static models are used to examine land-use change factors based on local knowledge.Other examples are inductive and deductive models.Inductive models are elaborated from correlations between land-use transitions and a selected set of variables that provide information on possible land change scenarios.Deductive models use theory as a basis for establishing patterns of change, developing the model based on the relationships between humans and the environment [4].
Changes in land use are accompanied by variations in their ecosystem value, which is why it is important to address methodologies that evaluate changes in the ecosystem valuation of natural resources; this type of method is based on assigning a monetary value to natural resources related to the production of environmental services [5].The environmental services (ESs) provided by the Andean zones to their environment according to the type of service are provisioning services, which are of direct use, and regulating services, cultural services and support services, which are of indirect use.Provisioning services are related to the food products that people obtain from the resource.Regulating services represent the contribution to the conservation of air and soil quality.Cultural services provide stability in the sense of belonging and spirituality, landscape aesthetics, customs, heritage and identity.Support services are made up of the processes necessary for the production of all other ecosystem services [13].
There are different methodologies for natural resource ESs, based on two approaches: (1) the anthropocentric approach and (2) the intrinsic approach [14].The first approach studies the value of natural resources from the utility that ESs provide to human beings, and the second approach analyzes the value of natural resources by affirming that every natural good is worth its own existence.The anthropocentric valuation approach divides the EV by the type of value into the categories of use value and option value.The use values allow determining the economic income obtained from provisioning services and are divided into the direct use value and the indirect use value.Indirect use values help to determine the value of support, regulatory and cultural ESs.Option values provide costs linked to future values.The intrinsic valuation approach is studied with the non-use value and the existence value.This vision is predominant in some Native peoples and communities with a total conservation vision [6].
Monitoring land-cover changes in Andean systems and analyzing the variability of their ecosystem valuation (EV) is a challenge that must be addressed in a comprehensive manner in order to develop sustainable strategies in accordance with the reality of the conservation status of the resource [15,16].Several studies have assessed land-cover change using the maximum similarity ranking and minimum distance ranking methods.The maximum similarity classification methodology is based on the assumption that the reflectivity data for each of the land-cover categories follow a multivariate normal distribution.The method allows setting probability thresholds for each category and rejecting pixels with low probability.However, the basic assumption of normality is not always met, so it is always necessary to verify it [17,18].The classification by minimum distance determines the types of coverage from the distance between the different pixels, where each pixel will be assigned to a class with respect to its minimum distance.This method is not very reliable because it overclassifies the image, i.e., no pixel is left unclassified [19,20].
Other researchers have used the parallelepiped method, which establishes limits for each class based on the maximum and minimum reflectivity values of each of the bands of a series of rectangles.The drawback of this method is that it does not allow modeling the dispersion of the control zones, due to the high correlation between the bands [21].Artificial neural networks have also been used in coverage monitoring.This type of methodology has been developed from the weighted sum of the data of the input variables to the network; by means of a threshold function, an output signal is produced, and the results are adjusted by means of back propagation algorithms.Some of the most commonly used neural networks are random forest (RF) and the classification and regression tree (CART).The disadvantages of neural networks are that in the regressions predictions cannot be made beyond the range of values of the training set, the predictions are not continuous in nature and the predictions do not allow constant control of the model [22].
The support vector machines (SVMs) methodology bases its land-use recognition technique on supervised learning classification and identifies the best hyperplane to achieve a clear separation between a set of data.Some researches [23][24][25] have focused on analyzing the performance of SVMs against other classifiers such as RF, CART and the parallelepiped method.The authors claim that the SVMs methodology outperforms classifiers.However, its performance is limited by the quality of its hyperparameter settings.
In relation to the analysis of ecosystem valuation based on land-use transitions, several studies have used methodologies that combine remote-sensing tools and monetary indices.For example, the contingent valuation method assigns an economic value to natural resources through the simulation of a hypothetical market by interviewing consumers, in which they state their maximum willingness to pay for a given good [26].The drawback of this method is that the responses of the interviewees do not always represent reality.Other studies have determined the monetary value of environmental services using the revealed preferences methodology, which is based on observable data from actual choices in markets.This approach analyzes the value of ecosystem goods and services based on people's choices according to value-theoretic principles.This means that the ES included will belong to a subset of the total value of the service [27].The scope of this method does not include non-use values.
Researchers have also used the valuation method considering the cadastral appraisal and the opportunity cost methodology.The valuation method considering cadastral valuation calculates the value of the natural asset using the physical and geographical characteristics of the property that refer to the existence value of the land as a resource of agricultural potential.It is common to consider that land with poor soil quality is equally distributed, and moving to subsequent categories increases inequity [28].The opportunity cost (OC) methodology assumes that the cost of using a resource can be approximated by the income foregone from other uses of the resource, i.e., it estimates how much income must be sacrificed and decides whether it is worthwhile [28].
For the land-cover study, this study used the SVMs method, due to its efficiency in land-use recognition, and made efforts to reduce the margin of error in the performance of the method based on evolutionary algorithms and effective validation error rates.Ecosystem valuation based on land-cover changes was determined using the opportunity cost methodology, including the relationships between the direct and indirect ESs of land use.
Based on the above, this research focused on the following objectives: (a) evaluate the conservation status of the areas of influence of the Ozogoche and Atillo del Ecuador lake systems, using remote-sensing techniques, geographic information systems and market indices; (b) optimize the performance of the support vector machines (SVMs) method, adjusting the configuration of its hyperparameters to improve land-cover recognition; (c) determine the ecosystem value of the studied resource in monetary terms, associating environmental services of direct and indirect regulatory use.
The information generated by this study will be of vital importance to understand the causes of the changes in the Andean systems in monetary and environmental terms, which will allow for the development of management plans and conservation policies aimed at protection and sustainable management from an economic approach.

Study Area
The research focused on the land cover of the Achupallas and Cebadas parishes, which are located in the province of Chimborazo, Ecuador, and constitute a zone of direct influence on the Ozogoche and Atillo lake systems (Figure 1).Its altitudinal range is between 2400 and 4730 masl, the average temperature is around 8 • C and the average annual rainfall is 1000 mm of water [1].The soils in the area of influence are characterized by a dark color and a pH between 3.9 and 5.4.It is estimated that the soil's weight can double due to the effect of water retention [2].

Workflow
The study was carried out through a combination of two stages (Figure 2): Stage 1: determination of land cover using spectral indices and support vector machines (SVMs) algorithm; Stage 2: evaluation of ecosystem valuation through the analysis of land-use changes and natural value in economic terms.

Data Processing
The work used images from Landsat 7 and Landsat 8 satellites through NASA's official website, https://earthexplorer.usgs.gov,accessed on 5 March 2023.Several scenes of the study area were downloaded from 2000 to 2020.Aerosols in the atmosphere were minimized by considering less than 10% cloudiness.The geometric accuracy of the images was verified using mapping information from the Ecuadoran Military Geographic Institute [29].The atmospheric correction and geometric verification procedures were performed using QGIS 3.4 software, which was downloaded from https://qgis.org/es/site/forusers/download.html,accessed on 25 February 2023 [30].

Workflow
The study was carried out through a combination of two stages (Figure 2): Stage 1: determination of land cover using spectral indices and support vector machines (SVMs) algorithm; Stage 2: evaluation of ecosystem valuation through the analysis of land-use changes and natural value in economic terms.

Data Processing
The work used images from Landsat 7 and Landsat 8 satellites through NASA's official website, https://earthexplorer.usgs.gov,accessed on 5 March 2023.Several scenes of the study area were downloaded from 2000 to 2020.Aerosols in the atmosphere were minimized by considering less than 10% cloudiness.The geometric accuracy of the images was verified using mapping information from the Ecuadoran Military Geographic Institute [29].The atmospheric correction and geometric verification procedures were performed using QGIS 3.4 software, which was downloaded from https://qgis.org/es/site/forusers/download.html,accessed on 25 February 2023 [30].

Workflow
The study was carried out through a combination of two stages (Figure 2): Stage 1: determination of land cover using spectral indices and support vector machines (SVMs) algorithm; Stage 2: evaluation of ecosystem valuation through the analysis of land-use changes and natural value in economic terms.

Data Processing
The work used images from Landsat 7 and Landsat 8 satellites through NASAʹs official website, https://earthexplorer.usgs.gov,accessed on 5 March 2023.Several scenes of the study area were downloaded from 2000 to 2020.Aerosols in the atmosphere were minimized by considering less than 10% cloudiness.The geometric accuracy of the images was verified using mapping information from the Ecuadoran Military Geographic Institute [29].The atmospheric correction and geometric verification procedures were performed using QGIS 3.4 software, which was downloaded from https://qgis.org/es/site/forusers/download.html,accessed on 25 February 2023 [30].

Reference Sites
Eight hundred control points were established for land-cover classification.The years included in the analysis were 2000, 2010 and 2020.These years were chosen because they represent the time with the greatest advance in the agricultural frontier [31].Seventy-five percent of the control points were field-checked.Land covers were established from the land-use catalog defined by the Ministry of Environment, Water and Ecological Transition of Ecuador [32].The following land covers were included: Moorland (Mo), Crop (C), Pasture (Pz), Wetland (We), Forest Plantation (PF), Water (W) and Soil (S).

Variables
The following set of variables were selected in order to study the land coverages included in the study: Normalized Difference Water Index (NDWI), Modified Difference Water Index (MNDWI), Automated Water Extraction Index (AWEI), altitude (DEM) and soil organic carbon (SOC).The altitudinal information (DEM) was taken from the official cartographic geoportal of the Military Geographic Institute of Ecuador through the site, https://www.geoportaligm.gob.ec/portal/,accessed on 10 April 2023 [33].GSOC data were selected from the Global Soil Organic Carbon Map version 1.5 data sources through the site https://www.fao.org/documents/card/en/c/ca7597,accessed on 10 April 2023 [34], and spectral data were calculated by means of Python in QGIS 3.4 software, through band clustering determined by the programming module [30].

Classification Methods
The classification was carried out in two stages.In the first stage, terrestrial cover was delimited from water bodies.A threshold of zero was established to separate terrestrial areas from aquatic areas [35].This stage was carried out due to the complexity of creating a precise boundary between terrestrial and aquatic cover.The second stage proceeded to classify the different soil types by means of the support vector machines (SVMs) methodology.The two classifications were carried out by means of the TerrSet v.19.0.8 software, which was downloaded on 15 May 2023 via https://clarklabs.org/download/[36].

Spectral Water Indices
The Normalized Difference Water Index (NDWI) was used to detect lake systems in wetland areas and to establish pond dimensions.Aquatic systems usually reflect positive values, whereas vegetation and soil usually reflect zero or negative values [37].The NDWI for Landsat 7 and 8 images was calculated by the following Equation ( 1): The Modified Difference Water Index (MNDWI) is efficient for extracting water bodies.The resulting values representing water features have positive values, and non-water features have negative MNDWI values [38].The MNDWI for Landsat 7 and 8 imagery was calculated by the following Equation ( 2): The Automated Water Extraction Index (AWEI) was used to differentiate water and non-water areas in shaded locations.Water sources reflect values above zero, while areas containing no water reflect negative values [39].The AWEI for Landsat 7 and 8 images was calculated by the following Equation (3): where SWIR 1 (Landsat 7: band 5; Landsat 8: band 6); SWIR 2 (band 7); NIR is near infrared (Landsat 7: band 4; Landsat 8: band 5); GREEN (Landsat 7: band 2; Landsat 8: band 3); BLUE (Landsat 7: band 1; Landsat 8: band 2).

Support Vector Machines (SVMs) Classification
Support vector machines (SVMs) are a non-parametric supervised classification method that can develop multiband images of optimal resolution from information from remotesensing tools.The essential basis of SVMs is the development of a hyperplane.By forming a hyperplane, the linked data are classified into different coverages [40].An SVM produces a concentrated boundary that divides the different classes along the hyperplane (Figure 3).If the points representing the categories are located above the hyperplane, they will be categorized as plus one, and if they are located below the hyperplane they will be categorized as minus one.Control zones that are close to the optimal hyperplane are called support vectors.The best way to train the SVMs algorithm is by means of linearly separable classes, as detailed in Equations ( 4) and ( 5) [41].
produces a concentrated boundary that divides the different classes along the hyperplane (Figure 3).If the points representing the categories are located above the hyperplane, they will be categorized as plus one, and if they are located below the hyperplane they will be categorized as minus one.Control zones that are close to the optimal hyperplane are called support vectors.The best way to train the SVMs algorithm is by means of linearly separable classes, as detailed in Equations ( 4) and ( 5) [41].Training points (Y , X ), (Y , X ) …(Y , X ) …Y ∈ [−1, +1] are recognized as linearly separable if there is a vector W and a scalar B for the inequalities [42], as described in Equations ( 4) and ( 5).
These inequalities are valid for the whole training set (Y , X ); k is the proportion of samples denoted as (Y , X ); x ∈ Rn, an n-dimensional space; y ∈ [−1, +1] represents the class label; W is the linear hyperplane.
To avoid multiple evaluations of the objects to be analyzed in the support vector machines (SVMs) algorithm, the error in the performance of the method was minimized in relation to the evolutionary algorithms of the hyperparameters and the bounds of the validation error bounds were improved.The SVMs algorithm was processed by the following expressions [42]: Training points ] are recognized as linearly separable if there is a vector W k and a scalar B k for the inequalities [42], as described in Equations ( 4) and ( 5).
These inequalities are valid for the whole training set (Y k , X k ); k is the proportion of samples denoted as (Y k , X k ); x ∈ Rn, an n-dimensional space; y ∈ [−1, +1] represents the class label; W k is the linear hyperplane.
To avoid multiple evaluations of the objects to be analyzed in the support vector machines (SVMs) algorithm, the error in the performance of the method was minimized in relation to the evolutionary algorithms of the hyperparameters and the bounds of the validation error bounds were improved.The SVMs algorithm was processed by the following expressions [42]: where x ∈ R d = Appl.Sci.2024, 14, x FOR PEER REVIEW 8 of 24 where x ∈ R = ℓ is a vector composed of elements to be analyzed; w is the vector of adjustable parameters computed from a set of training points; Φ is the mapping function over space;  represents the training data set.The signs of the vectors evaluated by means of Equation ( 6) were essential to guarantee a correct classification of the different classes.The settings of the parameter vector, w, were performed by training the classifier.The training sets (x ) and the associated classes were represented separately by (y ) ∈ (−1, +1) for i =1, ..., l; the vector of weights, w, was expressed by means of the linear combination of the different classes [43].The term α was found through the quadratic function, by means of Equation (8).
Equation ( 8) is subject to the constraint C ≥ α ≥ C, for i = 1, ..., l, where C is the value to be selected by the user, and was used to make the classifier tolerant to the lack of separability of the training data in space.The values α ≥ 0 were defined as the support vectors and represented the closest points to the decision boundary hyperplane in the space H.The inner product Φ ⋅ (x ) ⋅ Φ ⋅ (x ) in H was performed by using the function K(x , x ), known as kernel.Through this function it was possible to choose the high-level parameters, usually called hyperparameters [44].The selection of the hyperparameters was verified by means of the evolutionary algorithm of covariance matrix adaptation (CMA); each individual represented a dimensional vector with real values, which were updated by recombinations.This process involved calculating the center of mass of the is a vector composed of elements to be analyzed; w is the vector of adjustable parameters computed from a set of training points; Φ is the mapping function over space; x i represents the training data set.
The signs of the vectors evaluated by means of Equation ( 6) were essential to guarantee a correct classification of the different classes.The settings of the parameter vector, w, were performed by training the classifier.The training sets (x i ) and the associated classes were represented separately by ( y i ) ∈ (−1, +1) for i = 1, . .., l; the vector of weights, w, was expressed by means of the linear combination of the different classes [43].The term α was found through the quadratic function, by means of Equation (8).
Equation ( 8) is subject to the constraint C ≥ α i ≥ C, for i = 1, . .., l, where C is the value to be selected by the user, and was used to make the classifier tolerant to the lack of separability of the training data in space.The values α i ≥ 0 were defined as the support vectors and represented the closest points to the decision boundary hyperplane in the space H.The inner product Φ•(x i )•Φ• x j in H was performed by using the function K(x i , x j ), known as kernel.Through this function it was possible to choose the high-level parameters, usually called hyperparameters [44].The selection of the hyperparameters was verified by means of the evolutionary algorithm of covariance matrix adaptation (CMA); each individual represented a dimensional vector with real values, which were updated by recombinations.This process involved calculating the center of mass of the individuals in the population and recombining the distributed random vectors with zero mean.

Performance of Water Spectral Indices Using Pearson's r
The statistical measure of Pearson's r allows the strength and direction of a linear association between different images to be determined; if the association is not linear, then the coefficient is not adequately represented [45].Pearson's r value between two images was defined by the following equation: where µ A and µ B are the mean values of the two images (A, B), respectively.The correlation coefficient can take a range of values from +1 to −1.A value of 0 indicates that there is no association between the two variables.A value greater than 0 indicates a positive association; that is, as the value of one variable increases, so does the value of the other.A value less than 0 indicates a negative association; that is, as the value of one variable increases, the value of the other decreases [45].

Validation of Coverage Classification
The validation of the land-cover classification was carried out by means of the statistical measures of producer accuracy (PA), user accuracy (UA), overall accuracy (OA) and the Kappa index.The PA allowed establishing references by means of the predictions of the classes included in the study.Based on the initial predictions, the proportion of correct and incorrect classifications was determined [46].The UA made it possible to recognize the probability that a classified category effectively represents the class determined in the field.The OA contributed to the development of a base of successfully classified categories.The Kappa index made it possible to evaluate the relationship between the data observed in the image and the data identified by means of the classifier [47].
After the land-cover classification, a multitemporal analysis was developed by means of IDRISI's Crosstab module through the TerrSet software [48], generating change matrices between the different categories analyzed.For the case study, the changes were established in two periods (2000, 2010; 2010, 2020).The transition matrices generated by the changes in coverage made it possible to estimate losses, gains and permanence for each category.In addition, it was possible to establish the types of coverage transitions, which can be systematic or random.The equations for validating the classification and evaluating the temporal changes of the hedges are detailed in Table 1 [49].

Ecosystem Valuation
The ecosystem valuation (EV) was developed using the opportunity cost methodology based on the main productive activities and soil carbon concentration.The determination of opportunity cost values was based on indicators such as the net benefits (NBs) of each land-cover use and the net present value (NPV) of productive trajectories [50].
In order to obtain the aforementioned indicators, the costs, income, rate of return and net profit of the main land uses in the study area were defined.The land uses considered were crops and pasture.The calculations related to crops were calculated based on the relationship between the yield and the costs of the crop, the net profit was defined by means of the income and the realized investment rate and the rate of return was determined from the profit obtained by the worker.The monetary value of the use of pasture hedges was calculated by means of the dairy sector; most inhabitants have opted for this sector due to the stability of its price.The average values of milk generation per unit of pasture were determined by means of the nutritional parameters of the cattle and the properties of the Kikuyo [51] pasture, which is characteristic of the study site.The yield and profit values of the pastures were calculated using the same relationships as in the cultivated land uses.
The cost, income and profit base for crops and pasture was developed based on official information from the Ministry of Agriculture and Livestock (MAG) [51], the National Autonomous Institute of Agricultural Research (INIAP) [52], the Ministry of Environment, Water and Ecological Transition of Ecuador (MAATE) [32] and the Consumer Price Index (CPI) [53]; price variability and projections during the study period were considered.The equations used to develop the aforementioned calculations are detailed below [51].
Here, NP i,x is the net profit (USD); i is the land-use; Q is the quantity produced; P is the price (USD); CT is the total cost (USD).
Considering all the characteristic productive trajectories in each of the analysis zones, as well as the net benefits of each of the uses that comprise them, the net present value (NPV) was estimated.
Here, NPV is the net present value USD ha ; i is the use; x is the study area (ha); j is the production trajectory; (1 + r) t is the discount factor for year t; r is the discount rate; t is time in years.
The formulas described above made it possible to estimate the set of trajectories characteristic of each study area and to establish a comparative analysis of the different opportunity costs over time.The discount rate used to evaluate the costs was 12% [54]; this value is generally used in Ecuador.
To relate opportunity costs to carbon content, information from the organic carbon map (GSOC) [34] and the elevation map (DEM) [33] was used.Using the methodology of average carbon volumes [55], average carbon values were estimated for the different sequences of land use; it was assumed that with changes in land use there will be an emission equivalent to the difference between the total carbon contained in the study area and that which is maintained in the alternative use.Based on this difference, the amount of carbon released into the atmosphere by each type of land use was estimated; the carbon content released was transformed into CO 2 terms, considering that 3.67 tons of CO 2 emitted is equivalent to one ton of carbon [56].The equation used to relate opportunity costs to organic carbon degradation in soils is detailed below [55].
where OC is the opportunity cost ( USD tCO 2eq ); NPV is the net present value ( USD ha ); A is the carbon concentration ( C(t) ha ); B is a conversion factor of 3.67 ( co 2eq (t) C(t) ).

Performance of Water Spectral Indices by Means of Pearson's r
The efficiency of the spectral indices to detect water masses was measured by means of the Pearson's r statistical measure.Table 2 details the correlations between the different indices.The highest Pearson's r was observed between the NDWI and MNDWI in 2000, 2010 and 2020, with values of 0.95, 0.97 and 0.98, respectively.The lowest Pearson's r values were observed between the NDWI and AWEI in 2000, 2010 and 2020, with values of 0.80, 0.85 and 0.88, respectively.All Pearson's r correlations were equal to or greater than 0.80; therefore, their interpretation was identified as very good.

Accuracy of the Detection of the Analyzed Land Cover
Tables 3 and 4 detail the accuracy of the recognition of lakes through spectral indices and land cover.Accuracy was measured by means of producer accuracy, user accuracy, overall accuracy and Kappa index.In 2000, 2010 and 2020, through the analysis of the producer accuracy (PA) and user accuracy (UA), it was found that, both in the learning phase and the validation phase, the NDWI and the SVMs classifiers were the most accurate in detecting water bodies.Although the NDWI and SVMs classifiers obtained the best results, the NDWI presented fewer errors than the SVMs algorithm in the recognition of water and soil surfaces based on the control points added by the operator.The same performance trend was present in the recognition of classes based on the total number by the program.The efficiency of the SVMs algorithm in the recognition of the terrestrial categories, W, C, Mo, Pz, S, We and PF, was efficient in both PA and UA.
The overall accuracy (OA) indicated that the proportion of objects correctly classified by NDWI, NDWI and AWEI for the validation phase was 93%, 86% and 80%, respectively, with a Kappa of 92%, 85% and 78% for 2000, 2010 and 2020, respectively.For its part, the overall ranking of the SVMs algorithm in the classification coverages involved in the validation phase was 85%, 86% and 83% with a Kappa of 80%, 85% and 83% for 2000, 2010 and 2020, respectively.

Determination of Land-Cover Changes
Changes in land cover are closely linked to the economic development of the inhabitants.The loss of natural resources, such as moorland soils, wetlands and water bodies, increased during the study period, and the degradation trend continues.Tables 5 and 6 below detail the percentages of changes in the different categories analyzed; Tables 7 and 8 describe the main systematic losses and gains; Figure 4 shows the variability of the mapping; and Figure 5 shows the combination of transitions of all the coverages.From 2000 to 2010 (Table 5), the increase in land cover was presented as follows.Pasture soils (Pz) had the greatest increase, increasing from 11% to 16.4%, with a persistence of 9%, a gain of 7.4% and a loss of 2%.Crop cover (C) increased from 10.7% to 14.5%, with a persistence of 8%, a gain of 6.5% and a loss of 2.7%.Soil cover (S) increased from 1.2% to 2.1%, with a persistence of 0.9%, a gain of 1.2% and a loss of 0.3%.Forest plantations (PF) presented a slight increase from 1% to 1.3%, with persistence mostly represented by 8%, a gain of 0.5% and a loss of 0.2%.The degradation of the analyzed coverages was presented in the following order.Moorland soils (Mo) decreased from 54.4% to 46.1%, with a persistence of 45%, a gain of 1.1% and a loss of 9.4%.Water represented by lakes (W) degraded from 5.3% to 4.6%, with a persistence of 4%, a gain of 0.6% and a loss of 1.3%.Wetland (We) slightly lost coverage, from 16.4% to 15%.From 2000 to 2010 (Table 5), the increase in land cover was presented as follows.Pasture soils (Pz) had the greatest increase, increasing from 11% to 16.4%, with a persistence of 9%, a gain of 7.4% and a loss of 2%.Crop cover (C) increased from 10.7% to 14.5%, with a persistence of 8%, a gain of 6.5% and a loss of 2.7%.Soil cover (S) increased from 1.2% to 2.1%, with a persistence of 0.9%, a gain of 1.2% and a loss of 0.3%.Forest plantations (PF) presented a slight increase from 1% to 1.3%, with persistence mostly represented by 8%, a gain of 0.5% and a loss of 0.2%.The degradation of the analyzed coverages was presented in the following order.Moorland soils (Mo) decreased from 54.4% to 46.1%, with a persistence of 45%, a gain of 1.1% and a loss of 9.4%.Water represented by lakes (W) degraded from 5.3% to 4.6%, with a persistence of 4%, a gain of 0.6% and a loss of 1.3%.Wetland (We) slightly lost coverage, from 16.4% to 15%.
The changes in land cover from 2010 to 2020 (Table 6) were very similar to the transitions of the previous analyzed period.Pasture (Pz) maintained its predominance of gain with respect to the other coverages, increasing from 16.4% to 23%, with a persistence of 12%, a gain of 11% and a loss of 4.4%.In this period, crops (C) increased more than the previous period, increasing from 14.5% to 20%, with a persistence of 11%, a gain of 9% and a loss of 3.5%.Soils maintained a growth trend, increasing from 2.1% to 3.5%, with a persistence of 1.6%, a gain of 1.9% and a loss of 0.5%.Forest plantations (PF) had little growth, as the area expanded from 1.3% to 1.5%.The trend of land-cover degradation was accentuated in paramo ecosystems, with a gain lower than in the 2000-2010 period and a higher loss, decreasing from 46.1% to 35.7%.Water bodies maintained the trend of loss, gain and persistence, decreasing from 4.6% to 3.3%.Wetland (We) slightly reduced in this period, degrading by 2.5%, with a persistence of 12.5% and a gain of 0.5%.
The western zone is the area that shows the greatest degradation, possibly because accessibility to this area is not complicated, which has facilitated population growth, bringing with it an increase in the productive sector.Pasture soils represent the cover that has caused the greatest environmental impact on the natural systems of the study area, such as the moorlands and wetlands.This may have occurred because of the set of laws that provide stability to the prices of the dairy sector.After pasture, crops also represent, in good proportion, one of the main land covers that have caused degradation.The main  According to several authors [57,58], if the first-half category systematically gains against the existing second-half category and the second-half category systematically loses against the first-half category, the results indicate a systematic transition process from the first-half category to the second-half category.Therefore, it is critical to identify The changes in land cover from 2010 to 2020 (Table 6) were very similar to the transitions of the previous analyzed period.Pasture (Pz) maintained its predominance of gain with respect to the other coverages, increasing from 16.4% to 23%, with a persistence of 12%, a gain of 11% and a loss of 4.4%.In this period, crops (C) increased more than the previous period, increasing from 14.5% to 20%, with a persistence of 11%, a gain of 9% and a loss of 3.5%.Soils maintained a growth trend, increasing from 2.1% to 3.5%, with a persistence of 1.6%, a gain of 1.9% and a loss of 0.5%.Forest plantations (PF) had little growth, as the area expanded from 1.3% to 1.5%.The trend of land-cover degradation was accentuated in paramo ecosystems, with a gain lower than in the 2000-2010 period and a higher loss, decreasing from 46.1% to 35.7%.Water bodies maintained the trend of loss, gain and persistence, decreasing from 4.6% to 3.3%.Wetland (We) slightly reduced in this period, degrading by 2.5%, with a persistence of 12.5% and a gain of 0.5%.
The western zone is the area that shows the greatest degradation, possibly because accessibility to this area is not complicated, which has facilitated population growth, bringing with it an increase in the productive sector.Pasture soils represent the cover that has caused the greatest environmental impact on the natural systems of the study area, such as the moorlands and wetlands.This may have occurred because of the set of laws that provide stability to the prices of the dairy sector.After pasture, crops also represent, in good proportion, one of the main land covers that have caused degradation.The main crop that has been developed in the area is potatoes, perhaps because this type of crop adapts efficiently to the climatological and topographical characteristics of the Andean zones.
Soils have increased to a greater extent in the areas near the lakes; water bodies have decreased as a result of the increase in soil cover.This dynamic of land-cover changes could be occurring due to notable changes in temperature, possibly as a result of climate change.Temperature variability causes a decrease in the hydrological capacity of soils, causing dryness, the release of carbon into the atmosphere and the instability of related ecosystems.The alterations to the natural resources of the study area can be attributed to the lack of adequate agricultural and livestock planning.
According to several authors [57,58], if the first-half category systematically gains against the existing second-half category and the second-half category systematically loses against the first-half category, the results indicate a systematic transition process from the first-half category to the second-half category.Therefore, it is critical to identify systematic transitions based on gains and losses to verify the simultaneous impact and to recognize the key signals of coverage changes in the area.The main systematic transitions in terms of gains and losses that occurred between 2000 and 2020 are described as follows: Soil cover was gained and replaced water cover.Pasture cover was gained and replaced crops, forest plantations, soil, wetland and moorland.Pasture cover change is predominant over all the changes that occurred in the rest of the analyzed covers; although water was replaced by soil cover, soil was replaced by pasture.
The main transition forces are represented by changes from moorland to pasture and changes from water to soil.The least representative transition force values are related to changes from pasture to other types of cover, which can occur because, although pasture cover loses surface area, the growth of its area is much greater.The transitions in which the pasture category did not intervene were also systematic, but the lack of simultaneous occurrences determines that the values obtained neither show the key information to study the change in the analyzed territory nor do they provide a clear idea related to the coverages that were exposed to a greater pressure of change.
The transitions that have occurred indicate that the agricultural frontier is advancing and replacing natural ecosystems.There is no systematic transition in terms of gains that indicates that any natural cover is recovering in a noticeable way, so natural ecosystems are losing the environmental services they provide to their surroundings.The natural landscape has not had an efficient recovery period in recent years, and, if the transitions maintain their trend, it is possible that there will be a crisis in the social and environmental balance, as the soils will decrease their productive capacity, directly affecting economic development and, consequently, the social factor.

Evaluation of the Variability of the Ecosystem Valuation
The calculation of the ecosystem valuation was segmented according to three altitude ranges, which were 2500 to 3000, 3001 to 3500 and 3501 to 4000 masl.The economic indices that were related to carbon concentrations were determined from the main productive trajectories, which are associated with the most characteristic crops of the Andean zones and pastures in terms of dairy production in the area.Table 9 shows the changes in ecosystem valuation during the study period in terms of tons of carbon stored in the soil, tons of carbon dioxide emitted and economic values.Table 10 and Figures 6 and 7 show the general changes in ecosystem valuation.The EV decreases as a function of the loss of organic carbon in soil, loss of cover productivity and increase in investment costs.
Ecosystem valuation (EV) synergies are closely linked to changes in land cover, because these changes cause direct effects on the environmental services of natural resources.Ecosystem valuation was calculated for the following periods: 2000-2010; 2010-2020.These periods were chosen because they represent the time when the greatest agricultural expansion occurred in the Andean zone.
The results show that the degradation of ecosystem valuation occurred in an ascending manner according to altitudinal ranges.The lowest zones were the most affected, perhaps because these places are the ones that are most in contact with the population; this dynamic was maintained throughout the study period.In general, the decrease in ecosystem valuation occurred as follows: In 2000, the EV as a function of avoided carbon dioxide content was USD 8.00 × 10 6 .In 2010, it was USD 6.00 × 10 6 , and, in 2020, it was USD 5.00 × 10 6 .Ecosystem valuation (EV) synergies are closely linked to changes in land cover, because these changes cause direct effects on the environmental services of natural    Ecosystem valuation (EV) synergies are closely linked to changes in land cover, because these changes cause direct effects on the environmental services of natural Changes in ecosystem valuation in the different altitudinal levels occurred progressively.In 2000, at an altitude between 2500 and 3000 m above sea level, the carbon content was 150 tons per hectare, the carbon dioxide concentration was 551 tons per hectare and its EV in relation to avoided carbon dioxide content was USD 140 per hectare.At an altitude between 3001 and 3500 m above sea level, the carbon storage was 220 tons per hectare, the number of tons of carbon dioxide avoided was 807 and the EV was USD 96 per hectare in relation to carbon dioxide emissions avoided.At an altitude between 3501 and 4000> masl, the tons of carbon per hectare was 324 and the tons of carbon dioxide was 1189; the ecosystem value at this altitude was USD 65 per hectare.
In 2010, at an altitude between 2500 and 3000 m above sea level, the concentration of carbon and carbon dioxide was 130 and 477 tons, respectively; the ecosystem valuation decreased to USD 110 per hectare in terms of carbon dioxide mitigated.At an altitude between 3001 and 3500 m above sea level, there was carbon storage of 200 tons, with 734 tons of carbon dioxide, and the ecosystem valuation was valued at USD 71.At an altitude between 3501 and 4000 m above sea level, the ecosystem valuation was USD 45 with a concentration of 318 tons of carbon and 1167 tons of carbon dioxide avoided.The trend was maintained in 2020, where the ecosystem valuation was USD 98, 58 and 33 based on carbon dioxide avoided, with a carbon concentration of 100, 170 and 300 tons and a concentration of buffered carbon dioxide of 367, 624 and 1101 tons in the altitudinal ranges of 2500-3000, 3001-3500 and 3501-4000> masl, respectively.
Although the trend of loss of ecosystem valuation is constant, there are areas that maintain the natural value of the resources, as the study area is located next to protected areas, which is why the owners of these areas receive economic incentives to use the soils with a sustainable approach.Many of the inhabitants have opted to include ecotourism activities in their economic development, taking advantage of the richness of the Andean soils and lake systems.Through this type of strategy, soils have been able to move into a phase of restoration and recovery [4].The application of indicators based on opportunity costs can multiply the number of benefits achieved from limited budgets and even inform priorities in conservation strategies.Environmental impacts can be determined by opportunity costs through an analysis of the return on economic investment in conservation programs [59].

Discussion
Studies analyzing the influence of land-use transitions on the variability of ecosystem valuation (EV) in mountainous areas with complex topographic characteristics can be efficiently performed using satellite land information and economic indices from existing markets [60,61].Some authors [50,59] state that the historical detection of land-cover changes as a function of EV is still a challenging path, due to the lack of reliable and extrapolatable information on the economic value of environmental services (ESs); the scarcity of data causes complications in understanding the importance of natural resources from an economic perspective.
Despite the difficulties in developing reliable ES studies, this study, as well as other research [25,26], has shown that the use of remote-sensing tools and geographic information systems effectively facilitates the recognition of the trend of land-cover changes over time.The quantitative information of edaphic transitions allows the association of economic values linked to the productive activities of the study site by means of soil trajectory analysis methods, as is the case of the opportunity cost methodology [24,62].
For this case study, the land-cover classification was carried out in two phases: in the first phase, a boundary between water bodies and land systems was determined; in the second phase, the different land covers were detected and classified.In the determination of the extraction accuracy of water bodies, it was detected that accuracy problems can be relevant in areas where the background soils include areas of low albedo or shade.Shadows in satellite images cause errors due to reflectance patterns similar to those of water zones; this resemblance can reduce the extraction of aquatic systems.
In low-reflectance zones where there are shaded areas without water, simple detection methodologies will not be able to satisfactorily classify pixels in images with water and without water [36,63].Therefore, three different water indices were included in the present study: MNDWI, NDWI and AWEI.By means of Pearson's r, it was determined that the correlation between the indices was good; therefore, each index provided important information for the detection of water bodies.The index with the highest classification efficiency was the MNDWI, with an overall average accuracy of 91% and a Kappa of 89%.The application of zero as a threshold value can generate a more optimal accuracy in the extraction of water bodies.
In the second phase of classification, the control zones were sufficient to successfully classify the coverages of the study area; the SVMs classifier used in this phase obtained an overall accuracy of 84.81% and a Kappa of 82.94%.The adaptation of the covariance matrix in the adjustment of the SVMs hyperparameter selection improved the efficiency of the model, as it allowed the development of a robust architecture of the algorithm in terms of automation.The information obtained through this classification allowed for understanding of the state of conservation and degradation of Andean ecosystems.From the analysis, it was possible to understand that factors linked to the economic growth of the inhabitants are transcendental to understand the degradation of the natural resource.Changes in land use in the Andean Mountain Range should be understood as socioecological and economic changes, as viable conservation strategies can only be generated from an environmental, social and economic approach [27,64,65].
The variability of the EV of the study area was manifested as follows: In 2000, the EV, as a function of carbon dioxide content avoided, was USD 8.00 × 10 6 .In 2010, it was USD 6.00 × 10 6 and in 2020, it was USD 5.00 × 10 6 .The results per hectare in relation to carbon dioxide mitigated in 2000, 2010 and 2020 were USD 30, 226 and 188, respectively.The EV per hectare data obtained in this study are close to those developed in other research [5,66]; this may be due to the similarity of topographic, climatological and political characteristics between the two areas studied.In another study conducted in the Andean Highlands, the changes in the EV data compared to those generated in this study are important; this may be due to the difference in agricultural production, productivity level, conservation of environmental services and sowing planning [67].
The conservation or degradation of the ecosystem valuation of the studied cover implies that landowners must desist from their exploitation.The lost income will represent opportunities in terms of tons of carbon dioxide avoided to conserve the ecosystems; these values represent the opportunity costs linked to the payment for environmental services (PESs) [65].Based on the results of the present study, it was possible to spatially determine the changes in the EV, which made it possible to evaluate the most vulnerable areas and to locate estimated spatial values for the PESs.Identifying priority areas optimizes time and funds and increases the efficiency of environmental management tools.The natural ecosystems most affected by the loss of ecosystem valuation were moorlands, followed by water bodies and wetlands.This trend was maintained during the entire study period, i.e., during the last 20 years, so it could be expected that the development of the Andean landscape will maintain these criteria in the future if effective environmental protection measures are not taken.
This study adds knowledge in the following ways: Firstly, it has been shown that changes in ecosystem valuation can be spatially segmented according to environmental services of direct use, such as land use for food services, and environmental services of indirect use, such as carbon sequestration.The mapping of agricultural income based on environmental services is a useful spatial indicator for the responsible planning of natural landscapes.Secondly, an indicator of agricultural income based on monetary rent has been developed, which is temporally and spatially stable; based on opportunity costs, natural resources can be managed through rotational information of the agricultural system.Thirdly, a method of adjustment of the SVMs was developed for the classification of coverages and a more effective threshold was determined to extract water bodies, minimizing problems caused by the concentration of shadows.
Although the methodology used in the study is effective for extrapolation to areas with similar environmental and political characteristics, the technique may be susceptible to diminishing efficiency if drastic changes in production-related values, investment prices and expected returns occur.In addition, it is important to consider that the ecosystem valuation mapping information determined in the study period can be understood from different approaches, depending on the target sector on which the assessment is being made.

Conclusions
The methodology proposed in this study to evaluate the influence of the edaphic environment on the ecosystemic valuation of the zone of influence of the Ozogoche and Atillo lake systems in Ecuador was adequate to meet the objectives of the study.It was conceptually demonstrated via a quantifiable method that the information developed through the evaluation of land use is crucial to analyzing the variability of ecosystem valuation and to determining the conservation status of natural resources in monetary and environmental terms.Through this methodology, it will be possible to obtain, in a practical way, data on the conservation status of the resource over time, allowing us to solve problems related to the scarcity of data and leading to the understanding of changes in the area from a socioecological approach, i.e., covering the environmental impacts of human activities on natural systems.The basis of the developed method allows for replication of the methodology.The performance of the classifiers in recognizing the analyzed land cover proved to be efficient.The extraction of water bodies by means of the MNDWI had an accuracy rate of 91%, and the support vector machines (SVMs) algorithm for the recognition of coverage in general had an overall accuracy of 85%.The adjustment made to the SVMs algorithm to improve the selection of hyperparameters was efficient, as it allowed the development of a robust algorithm architecture in terms of automation.
Between 2000 and 2020, moorlands, water bodies and wetlands have degraded by 19%, 2% and 3.4%, respectively.The ecosystem value of the study area in 2000 in terms of avoided carbon dioxide content was USD 8.00 × 10 6 ; in 2010, it was USD 6.00 × 10 6 ; and in 2020, it was USD 5.00 × 10 6 .The natural ecosystems that have suffered the greatest impact from the loss of cover are the moorlands, followed by water bodies and wetlands; this trend was maintained throughout the study period, that is, during the last 20 years, so it could be expected that the development of the Andean landscape will maintain these criteria in the future.The most affected areas are those that are most accessible to the inhabitants of the area, which is why anthropogenic activity is recognized as an important factor in the degradation of resources.
The proposed opportunity cost method was efficient in evaluating the trajectories of food production associated with the environmental services (ESs) of carbon supply and concentration and associated with indirect and support ESs, so it was determined that the information coming from the dynamics of environmental services is an important indicator to understand the real state of the resource and consequently allows for understanding the necessary strategies to undertake sustainable actions of ecological management.

Figure 1 .
Figure 1.Location of the study area.

Figure 2 .
Figure 2. Flowchart of the main stages of the study.

Figure 1 .
Figure 1.Location of the study area.

Figure 1 .
Figure 1.Location of the study area.

Figure 2 .
Figure 2. Flowchart of the main stages of the study.

Figure 2 .
Figure 2. Flowchart of the main stages of the study.

Table 1 .
Statistical measures used to analyze the performance of classification models and temporal analysis of changes in coverages.

Table 1 .
Cont.P j+ − P jj Degradation : L ij = P +j − P jj Persistence := P +j − Gain D ij : correctly classified pixels; C i : total classified classes corresponding to a coverage; N: total number of failures; Pr (a) : proportion of data matching in both strata; Pr (e) : probability that data agree in both strata; P +j : gains at time 2; P jj : no change; P j+ : degradation at time 2.

Table 3 .
Accuracy of lake classification using spectral indices.

Table 4 .
Accuracy of lake and land-cover classification using the SVMs classifier.

Table 5 .
Changes in land cover between 2000 and 2010.

Table 6 .
Changes in coverage between 2010 and 2020.

Table 7 .
Systematic transitions in terms of gains from hedging changes.Crop wins.Crop does not replace moorland.

Table 8 .
Systematic transitions in terms of losses caused by coverage changes.
Crop loses.Soil does not replace the crop.

Table 8 .
Cont.Crop loses.The moorland does not replace the crop.Wetland is lost.Moorland does not replace wetland.Moorland is lost.Crop does not replace moorland.
Moorland is lost.Water does not replace moorland.

Table 9 .
Changes in ecosystem valuation at altitudinal level.

Table 10 .
General variability of ecosystem valuation.

Table 10 .
General variability of ecosystem valuation.

Table 10 .
General variability of ecosystem valuation.