Suitable Land-Use and Land-Cover Allocation Scenarios to Minimize Sediment and Nutrient Loads into Kwan Phayao, Upper Ing Watershed, Thailand

: Human activity and land-use changes have affected the water quality of Kwan Phayao, Upper Ing watershed, due to the associated high sediment load and eutrophication. This study aims to identify suitable LULC allocation scenarios for minimizing sediment and nutrient export into the lake. For this purpose, the LULC status and change were ﬁrst assessed, based on classiﬁed LULC data in 2009 and 2019 from Landsat images, using the SVM algorithm. Later, the land requirements of three scenarios between 2020 and 2029 were estimated, based on their characteristics, and applied to predict LULC change using the CLUE-S model. Then, actual LULC data in 2019 and predicted LULC data under three scenarios between 2020 and 2029 were used to estimate sediment and nutrient export using the SDR and NDR models. Finally, the ecosystem service change index identiﬁed a suitable LULC allocation for minimizing sediment or/and nutrient export. According to the results, LULC status and change indicated perennial trees and orchards, para rubber, and rangeland increased, while forest land and paddy ﬁelds decreased. The land requirements of the three scenarios provided reasonable results, as expected, particularly Scenario II, which adopts linear programming to calculate the land requirements for maximizing ecosystem service values. For sediment and nutrient export estimation under the predicted LULC for the three scenarios, Scenario II led to the lowest yield of sediment and nutrient exports, and provided the lowest average ESCI value among the three scenarios. Thus, the LULC allocation under Scenario II was chosen as suitable for minimizing sediment or/and nutrient export into Kwan Phayao. These results can serve as crucial information to minimize sediment and nutrient loads for land-use planners, land managers, and decision makers.


Introduction
Freshwater ecosystems are considered as one of the essential natural resources for living organisms. The rate of deterioration of the water quality of freshwater resources, such as lakes, ponds, and rivers, has become a global problem [1]. Lakes are vital components of our planet's hydrological cycle, providing significant social and ecological functions while storing water and supporting aquatic biodiversity [2]; however, they are often the final recipients of nutrients discharged from adjacent uplands and wetlands. The management of a lake means managing its watershed-which are often mismanaged and challenged natural resources. Some problems originate in the lake itself, but most problems originate from activities occurring on the surrounding land [3]. Several studies have shown that land uses play essential roles determining the water quality of lakes. In particular, urban, built-up, and cultivated areas significantly influence the water quality when within the  The problems related to Kwan Phayao have recently been reported by the Department of Fisheries and Royal Irrigation Department, as summarized in Table 1. The Ministry of Agriculture and Cooperative has spent about 4 million USD to increase water storage by sediment dredging, aquatic plant weeding, and flap-gate-weir building [20].

Materials and Methods
The research methodology consisted of data collection and preparation, in five significant components ( Figure 2).

Data Collection and Preparation
The required input data for data analysis included GIS, remote sensing, and relevant data, as summarized in Table 2.

Data Collection and Preparation
The required input data for data analysis included GIS, remote sensing, and relevant data, as summarized in Table 2. Annual and monthly rainfall between 2011 and 2019 TMD C-factor and P-factor [22] Predicted rainfall between 2020 and 2029 [24] Nitrogen/Phosphorus loading [25] Nitrogen retention coefficients [25] Phosphorus retention coefficients [25]  What's Best in Excel CLUE-S model [26] SDR and NDR models of InVEST software NatCap Ecosystem Services Change Index (ESCI) [27] Note

LULC Classification and Change Detection
We downloaded Landsat 5 TM data in 2009 and Landsat OLI data in 2019 from the USGS website (www.earthexplorer.usgs.gov: accessed date 1 May 2019), which were applied to classify LULC data by conducting supervised classification using the support vector machine (SVM) algorithm of the EnMap-Box software. Two training areas for ten LULC types were separately prepared in order to define an optimal hyperplane for LULC classification using SVM. The scaled reflectance data of Landsat data (visible, NIR, and SWIR bands) and additional bands, Normalized Difference Vegetation Index (NDVI) [28], Normalized Difference Moisture Index (NDMI) [29], Soil-Adjusted Vegetation Index (SAVI) [30], Modified Normalized Difference Wetness Index (MNDWI) [31], and DEM were applied to classify the LULC map in 2009 and 2019. The optimized model parameters were provided by a grid search: namely, Gaussian radial basis function kernel, which required the variables gamma (γ), which defines the width of the Gaussian, and the regularization parameter (C), which controls the trade-off between the maximization of the margin between the training data vectors and the decision boundary plus the penalization of training errors [32].
After classification, the LULC maps in 2009 and 2019 were assessed, in terms of thematic accuracy, based on reference information from very high spatial resolution images from Google Earth in 2009/2010 and a field survey in 2020, respectively. The number of sample points used for the accuracy assessment was 788, based on multinomial distribution, and sample points were allocated using the stratified random sampling technique suggested by [33]. Final LULC maps in 2009 and 2019 were further used to detect LULC change, using a post-classification comparison algorithm to describe from-to change information between 2009 and 2019, as suggested by [34,35].

Land Demand Estimation of Three Different Scenarios
The land requirements under three different scenarios for LULC prediction between 2020 and 2029 were estimated based on their characteristics, as follows: (1) Scenario I (Trend of LULC evolution): The land requirements were calculated based on the annual change rate of LULC between 2009 and 2019 from the transition area matrix using the Markov Chain model.
(2) Scenario II (Maximization of ecosystem service values): The land requirements were calculated based on LULC allocation for maximizing ESV using Linear Programming (LP) with the Simplex method. The objective function and constraints were solved to maximize ecosystem service values using Equation (1), based on the coefficient values for the ESs of each LULC type (Table 3).

(3) Scenario III (Economic crop zonation):
The suitable zonation of four economic crops (paddy field, field crop, para rubber, and perennial trees and orchards) of LDD in 2018 was updated using the existing LULC data for 2019 in order to estimate the land requirements.

LULC Prediction of Three Different Scenarios
LULC prediction between 2020 and 2029 for the three scenarios was conducted using the CLUE-S model. In practice, the selected driving factors of LULC change (soil drainage, distance to stream, distance to water body, distance to village, slope, distance to road, distance to fault line, annual rainfall, elevation, income per capita, and population density at subdistrict level), as suggested by [37], were first examined in terms of their multicollinearity, using the variance inflation factor (VIF) to prevent the correlation of driving factors. As a general rule of thumb, the VIF value should not exceed 10 [38,39]. Then, binary logistic regression analysis was performed to identify the significant driving factors for specific LULC type allocations (Equation (2)): Log P i 1 − P i = β 0 +β 1 X 1,i +β 2 X 2,i . . . . . . +β n X n,i , where Pi is the probability of a grid cell for the considered land-use type at location i, and the X values are the location factors. The coefficients (β) are estimated through logistic regression, using the actual land-use pattern as the dependent variable [40]. After that, two sets of local parameters for LULC prediction using the CLUE-S modelnamely, the conversion matrix and elasticity of LULC change-were prepared. These parameters were considered and set up based on the transition probability matrix of LULC data between 2009 and 2019. The conversion matrix, which indicates the LULC change opportunity among LULC types, is assigned a value of 1 when it is allowed, or 0 when it is not allowed. Meanwhile, the elasticity values imply the probability of land-use change, which ranges from 0 (easily converted) to 1 (irreversible change), and are set up according to the transition probability matrix in the past period [26,41]. In this study, elasticity values were assigned according to the transition probability matrix of LULC change between 2019 and 2029 by the Markov Chain model, as suggested by Ongsomwang and Iamchuen [42].
Finally, the conversion matrix, the elasticity of LULC change, and land requirements under different scenarios were combined to predict LULC change data between 2020 and 2029, according to the driving factors of LULC change for specific LULC type location preference.

Ecosystem Services Assessment: Sediment and Nutrient Export
The base year LULC in 2019 and the predicted LULC between 2020 and 2029 of three scenarios, as primary input data, were used to estimate sediment and nutrient exports through the SDR and NDR models in the InVEST software suite.

Sediment Export Estimation
The sediment export is the amount of sediment eroded in the watershed from overland sources and delivered to the stream. In principle, the SDR model is first applied to calculate the amount of annual soil loss, using the Revised Universal Soil Loss Equation (RUSLE) of [43]: where A i is the annual soil erosion (ton. ha −1 yr −1 ), R i is the rainfall erosivity (MJ mm ha −1 h −1 y −1 ), K i is the soil erodibility (ton·ha·hr (MJ·ha·mm) −1 ), LS i is the slope lengthgradient factor, C i is the crop-management factor, and P i is a support practice factor for erosion control. The required input data for annual soil erosion are summarized below.
(1) The Rainfall erosivity factor (R) was calculated based on monthly rainfall data, as suggested by [44], as: where R is the rainfall erosivity factor (MJ mm ha −1 h −1 y −1 ), P i is the monthly rainfall (mm), and P is the annual rainfall (mm). We collected monthly rainfall data between 2011 and 2019 from TMD for model calibration and validation and actual sediment export estimation in 2019. Simulated rainfall data between 2020 and 2029 were collected from the National Center for Atmospheric Research (NCAR), for sediment export estimation under three scenarios.
(2) The Soil erodibility (K) erodibility was extracted from the soil series data of LDD, whereas its value under the slope complex in the soil map was extracted from the geology unit (Table A1 in Appendix A).
(3) The Slope length gradient factor (LS) was calculated from the DEM with a method developed by [45], as follows: where S i is the slope factor for a grid cell, calculated as a function of the slope radians θ (S = 10.8·sinθ + 0.03 where θ < 9%, or S = 16.8·sinθ − 0.50, where θ ≥ 9%); A i−in is the contributing area (m 2 ) at the inlet of a grid cell, which is computed using the d-infinity flow direction method; D is the grid cell linear dimension; x i = |sin α i | + |cos α i |, where α i is the aspect direction for grid cell i, and m is the RUSLE length exponent factor. (4) The crop management (C) and support practice (P) factors for erosion control, according to LULC type, are summarized in Table A2 in Appendix A. Then, the model calculates the sediment delivery ratio using a connectivity index (IC), threshold flow accumulation, and maximum SDR to indicate sediment retention. The SDR value was calculated, as suggested by [46], as: where SDR max is the maximum theoretical SDR, set to an average value of 0.8 [47], and IC 0 and k are calibration parameters that define the shape of the SDR-IC relationship as a Sigmoid function. Finally, the sediment reaches the stream at the outlet of the Upper Ing watershed. The total sediment export was calculated from the sum of the amount of annual soil loss multiplied by the sediment delivery ratio [48]: where E i is the sediment that erodes from any LULC that reaches the stream.

Nutrient Export Estimation
Under the NDR model, the nutrient loads are first defined based on the LULC map and associated loading rates. Then, the model calculates LULC-based loads and the runoff potential index to approximate modified loads, which are divided into sedimentbound (surface flow) and dissolved parts (subsurface flow). After that, nutrient delivery is computed for surface NDR and subsurface NDR based on the properties of pixels belonging to the same flow path (particularly the slope and retention efficiency of the land use) [48]. The NDR model requires specific factors and parameters to run the model. Some factors in this model are the same as those in the SDR model, including DEM, LULC, TFA, and K b . Other parameters, including nitrogen and phosphorus loads, nutrient runoff proxy, maximum retention efficiency, and critical length, were added. The specific required factors of the NDR models are summarized below.
(1) The nutrient runoff proxy (RP) is used to calculate the runoff potential index. The runoff proxy was interpolated, based on annual rainfall, using the inverse distance weighted (IDW) method.
(2) Land-use and land-cover (LULC) data represent the influence of the nutrient delivery to the stream. LULC data were used as input data to assign the nutrient loading for each LULC class (loads), as summarized in Table 4.
(3) The maximum retention efficiency (eff) indicates the proportion of nutrient retention by vegetation. The value for each LULC class varies between zero and one (see Table 4).
(4) The critical length (crit_len) is the distance that a patch of LULC retains nutrients in its maximum capacity. The critical length ranges from 10 to 300 m [49,50]. In this study, the value was first set to the pixel resolution, with a value of 30 m (see Table 4).
(5) The subsurface proportion (proportion_subsurface) value is the proportion of dissolved nutrients that travels by surface and subsurface flow. We set its value to zero, which indicates that all nutrients are delivered by surface flow (see Table 4).
After that, the NDR model was computed to transport nutrients by surface flow. The surface nutrient delivery ratio is the product of a delivery factor, representing the ability of downstream pixels to transport nutrients without retention, and was calculated using Equation (8) [48]: where NDR 0,i is the proportion of a nutrient that is not retained by downstream pixels (which is based on the maximum retention efficiency of the land between a pixel and the stream), IC i is a topographic index, and IC 0 and k are calibration parameters that define the shape of the NDR-IC relationship. Finally, the total nutrient export at the outlet of the watershed is estimated, from the sum of the product of the load and the NDR [48], as: where x expi is the nutrient export from any LULC, as the product of the load and the NDR. Furthermore, the performance of the SDR and NDR models was examined based on data observed by the PCD in the calibration period (2011)(2012)(2013)(2014)(2015) and validation period (2016-2018) with standard scale, using the coefficient of determination (R 2 ) and percent bias (PBIAS), as suggested by [51,52] (see Equations (10) and (11) and Table 5).
where X o is the observed export value at station i, X o is the average of observed export value over the validation period, X s is the simulated export value at station i, and X s is the average simulated export value over the validation period. Furthermore, i is the number of stations, and n is the total count of data pairs. The value of R 2 varies from 0 to 1.
where X o i is an observed export value at time step i, and X s i is a simulated export value at time step i.

Suitable LULC Allocation Scenario to Minimize Sediment and Nutrient Export
To assess the state of change in ecosystem services (i.e., sediment and nutrient export due to LULC change), the ecosystem services change index (ESCI) was applied to assess the ecosystem service states (ES), as proposed by [27]: where ESCI x is the ecosystem service change index of service x, and ES CURx j and ES HISx i are the current and historic ecosystem service state values of service x at times j and i, respectively. In this study, the ecosystem service change indices of sediment and nutrient export ecosystem service values for LULC in 2019 (as the base year), and those of the predicted LULC between 2020 and 2029 were separately calculated in a pairwise manner using the ESCI, then averaged to identify the suitable LULC allocation under the proposed scenarios to minimize sediment or/and nutrient export.

LULC Classification and LULC Change Detection
The results of LULC classification in 2009 and 2019 using the SVM algorithm are presented in Table 6 and Figure 3. As a result, the top three dominant LULC types in 2009 and 2019 were forest land, paddy field, and perennial trees and orchards. Meanwhile, the top three least dominant LULC types were miscellaneous land, para rubber, and rangeland in 2009 and miscellaneous land, wetland, and para rubber in 2019.
According to a simple comparison of LULC change in the area, the annual change rate and percentage of change between 2009 and 2019 are reported in Table 7. The dominant increasing areas of LULC types were perennial trees and orchards, para rubber, and rangeland, with annual change rates of 1.95, 1.64, and 1.05 km 2 , respectively. Contrarily, the major decreasing areas of LULC classes in the same period were forest land and paddy fields, with annual change rates of 3.98 and 2.04 km 2 , respectively. The primary cause of change areas was the conversion of paddy fields into perennial trees and orchards and the expansion of agricultural areas into forest land.

Driving Factors of LULC Change
The multicollinearity test among the physical and socio-economic factors (independently available) was conducted using the VIF, as reported in Table A3 in Appendix B. As a result, all selected driving factors were found to be insignificantly correlated, as

Driving Factors of LULC Change
The multicollinearity test among the physical and socio-economic factors (independently available) was conducted using the VIF, as reported in Table A3 in Appendix B. As a result, all selected driving factors were found to be insignificantly correlated, as the VIF values did not exceed 10. They were further used to analyze specific LULC type location preferences by binary logistic regression analysis, as shown in Table 8. According to the result, the most significant driving factor for all LULC type allocation in the study area was soil drainage. The second most important vital driving factors were slope and annual rainfall. The third most crucial driving factors for the LULC type allocation area were distance to the road and elevation. In the meantime, other factors played crucial roles in specific LULC types, as shown in the table. Nevertheless, the distance to fault line and the income per capita at the subdistrict level were insignificant driving factors for all LULC type allocations in the study area.

Local Parameter of CLUE-S Model for LULC Prediction
The conversion matrix of LULC change for LULC prediction of three scenarios, which were assigned by considering the transition probability matrix of LULC data between 2009 and 2019 (Table 9), is presented in Tables 10-12. Meanwhile, elasticity values, as probability values for the urban and built-up area, paddy field, field crop, para rubber, perennial trees and orchards, forest land, waterbody, rangeland, wetland, and miscellaneous land were 1.00, 0.84, 0.22, 0.37, 0.67, 0.92, 0.93, 0.38, 0.45, and 0.29, respectively.

Land Requirement Estimation and LULC Prediction under Scenario I
The land requirement estimation for Scenario I (Trend of LULC evolution) was calculated based on the annual rate of LULC change from the transition area matrix between 2019 and 2029, using the Markov Chain model presented in Table 13. The significant increase in land requirements in the predictive period was observed in perennial trees and orchards, para rubber, water body, rangeland, urban and built-up area, field crop, and miscellaneous land. In contrast, the land requirements decreased for forest land, paddy field, and wetland.    The results of LULC prediction data for Scenario I, which were simultaneously allocated based on elasticity values and conversion matrix of LULC change (Tables 9 and 10), the land requirements (Table 13), and driving factors on LULC change for specific LULC type location preference (Table 8), are presented in Table 14 and Figure 4. Meanwhile, the deviation between the estimated land requirements and predicted LULC data in 2029 under Scenario I is reported in Table A4 in Appendix B.

LULC in 2019
Urban and built-up area (UR) Note: 0 is not allowed, and 1 is allowed. As a result, the LULC prediction between 2020 and 2029 was determined in terms of the driving factors of LULC change and local parameters (conversion matrix and elasticity values), particularly land requirements, which were estimated based on the transition area matrix between 2019 and 2029 using the Markov Chain model. The derived predicted LULC data correspond to the definition of Scenario I, which allows for LULC change (decreased or increased area) according to the trend of LULC evolution from 2009 to 2019, to represent socio-economic development in the study area.  As a result, the LULC prediction between 2020 and 2029 was determined in terms of the driving factors of LULC change and local parameters (conversion matrix and elasticity values), particularly land requirements, which were estimated based on the transition area matrix between 2019 and 2029 using the Markov Chain model. The derived predicted LULC data correspond to the definition of Scenario I, which allows for LULC change (decreased or increased area) according to the trend of LULC evolution from 2009 to 2019, to represent socio-economic development in the study area.
There was a slight difference between the required land area and the predicted area of each LULC type in 2029 under Scenario I (see Table A4 in Appendix B). The deviation values varied between −0.0002% (0.02 km 2 ; underestimation) to 0.0002% (0.02 km 2 ; overestimation). The deviation values depend on the iterative driving factor that determines the highest probability that each spatial unit will be converted to a specific land-use type in the following year [53,54]. Nevertheless, the summation of deviation values (i.e., considering the trade-off between overestimation and underestimation among LULC types) was 0.00%. Therefore, the LULC prediction under Scenario I can be considered acceptable for estimating sediment and nutrient export in the study area. There was a slight difference between the required land area and the predicted area of each LULC type in 2029 under Scenario I (see Table A4 in Appendix B). The deviation values varied between −0.0002% (0.02 km 2 ; underestimation) to 0.0002% (0.02 km 2 ; overestimation). The deviation values depend on the iterative driving factor that determines the highest probability that each spatial unit will be converted to a specific land-use type in the following year [53,54]. Nevertheless, the summation of deviation values (i.e., considering the trade-off between overestimation and underestimation among LULC types) was 0.00%. Therefore, the LULC prediction under Scenario I can be considered acceptable for estimating sediment and nutrient export in the study area.

Land Requirement Estimation and LULC Prediction under Scenario II
The land requirement estimation under Scenario II (Maximization of ecosystem service values), which was estimated based on the annual rate change between classified LULC data in 2019 and allocated LULC data in 2029 after maximization of ecosystem services according to the objective function (see Equation (1)) and constraining decision variables (see Table A5 in Appendix B) using LP through the Simplex method, is reported in Table 15.
As a result, the increased LULC classes under this scenario were urban and built-up area, para rubber, perennial trees and orchards, and wetland, with increasing annual rates of 0.38, 0.87, 1.95, and 1.34 km 2 , respectively. In contrast, the decreased LULC classes in 2029 were paddy fields, forest land, rangeland, and miscellaneous land, with annual rates of decrease of 0.99, 2.06, 1.36, and 0.14 km 2 , respectively. Meanwhile, field crops and water bodies were unchanged.
The results of LULC prediction data between 2020 and 2029 under Scenario II, which were simultaneously allocated based on elasticity values and the conversion matrix of LULC change (Tables 9 and 11), the land requirements (Table 15), and driving factors of LULC change for specific LULC type location preference (Table 8), are presented in Table 16 and Figure 5. Meanwhile, the deviation between the estimated land requirements and predicted LULC data in 2029 under Scenario II is reported in Table A6 in Appendix B.  As a result, the significant LULC types with increasing area were urban and built-up area, para rubber, perennial trees/orchards, water body, and wetlands, but the several LULC types with decreased area were paddy fields, field crops, forest land, rangeland, and miscellaneous land. In particular, wetland-which provides the highest coefficient value for ecosystem services-is expected to increase in the future. Wetland areas increased from 16 Table A6 in Appendix B). Hence, the LULC prediction under Scenario II can be considered As a result, the significant LULC types with increasing area were urban and built-up area, para rubber, perennial trees/orchards, water body, and wetlands, but the several LULC types with decreased area were paddy fields, field crops, forest land, rangeland, and miscellaneous land. In particular, wetland-which provides the highest coefficient value for ecosystem services-is expected to increase in the future. Wetland areas increased from 16.41 km 2 in 2019 to 29.72 km 2 in 2029, coming from paddy fields (9.50 km 2 ), rangeland (2.83 km 2 ), and miscellaneous land (0.98 km 2 ) in 2019. This result implies the efficacy of linear programming for determining reclaimed and other areas changing into the wetland.
Like Scenario I, there was a slight difference between the required land area and the predicted area of each LULC type in 2029 under Scenario II. The deviation values varied from −0.0013 to 0.0033%; nevertheless, the summation of deviation values was 0.00%. (See Table A6 in Appendix B). Hence, the LULC prediction under Scenario II can be considered acceptable for estimating sediment and nutrient export in the study area.

Land Requirement Estimation and LULC Prediction under Scenario III
Land requirement estimation under Scenario III (Economic crop zonation) was estimated based on areas of suitability classes for economic crops according to the LDD and Markov Chain model, as shown in Table 17.
As a result, the increased LULC classes were paddy field, field crop, water body, rangeland, urban and built-up area, and miscellaneous land, with annual increase rates of 5.96, 1.21, 0.80, 0.61, 0.38, and 0.10 km 2 , respectively. In contrast, the decreased LULC classes were perennial trees and orchards, forest land, wetland, and para rubber, with annual decrease rates of 5.23, 3.65, 0.16, and 0.01 km 2 , respectively.
The results of LULC prediction data between 2020 and 2029 under Scenario III, which were simultaneously allocated based on elasticity values and the conversion matrix of LULC change (Tables 9 and 12), the land requirements (Table 17), and driving factors of LULC change for specific LULC type location preference (Table 8), are presented in Table 18 and Figure 6. In the meantime, the deviation between the estimated land requirements and predicted LULC data in 2029 under Scenario III is reported in Table A7 in Appendix B. As a result, under Scenario III, economic crop areas, precisely paddy field, field crop, para rubber, and perennial trees and orchards were located based on suitability classes of economic crop zonation. Paddy fields expanded from 220. 74

Sediment Export Estimation Using SDR Model
In general, the SDR model reports soil erosion, sediment retention, sediment deposition, and sediment export. For this study, sediment export was selected to describe ecosystem services in the study area by each LULC allocation scenario (Scenarios I-III).
The statistical performance of the model under the calibration and validation periods provided a high correlation between the observed and estimated sediment export, with R 2 values of 0.697 and 0.824, respectively, indicating the good and very good fitting performance rate of the model for nitrogen export estimation, as suggested by [51,52]. These findings are consistent with the previous study of [55], who applied the SDR model to analyze sediment at the Rmel river basin, and found that the correlation between the estimated and observed sediment export in the calibration process was relatively high, with R 2 values of 0.84 and 0.706.
Meanwhile, the PBIAS values under the calibration period provided a good fit between the estimated and observed sediment export, with a value of −27.03%, while the PBIAS value under the validation period delivered an unsatisfactory fit between the estimated and observed sediment export, with a value of 65.60%, in accordance with [51,52]. The optimum local model parameters and domain values under the calibration period of the SDR model are summarized in Table 19.

Default
Minimum Maximum Adjusted Optimum As with Scenario I, there was a slight difference between the required land area and the predicted area of each LULC type in 2029 under Scenario III. The deviation values varied from −0.0023 to 0.0015%; however, the summation of deviation values was 0.00%. (See Table A7 in Appendix B). Hence, the LULC prediction under Scenario III can be considered acceptable for estimating sediment and nutrient export in the study area.

Sediment Export Estimation Using SDR Model
In general, the SDR model reports soil erosion, sediment retention, sediment deposition, and sediment export. For this study, sediment export was selected to describe ecosystem services in the study area by each LULC allocation scenario (Scenarios I-III).
The statistical performance of the model under the calibration and validation periods provided a high correlation between the observed and estimated sediment export, with R 2 values of 0.697 and 0.824, respectively, indicating the good and very good fitting performance rate of the model for nitrogen export estimation, as suggested by [51,52]. These findings are consistent with the previous study of [55], who applied the SDR model to analyze sediment at the Rmel river basin, and found that the correlation between the estimated and observed sediment export in the calibration process was relatively high, with R 2 values of 0.84 and 0.706.
Meanwhile, the PBIAS values under the calibration period provided a good fit between the estimated and observed sediment export, with a value of −27.03%, while the PBIAS value under the validation period delivered an unsatisfactory fit between the estimated and observed sediment export, with a value of 65.60%, in accordance with [51,52]. The optimum local model parameters and domain values under the calibration period of the SDR model are summarized in Table 19.  As a result, the highest sediment export occurred on miscellaneous land, with an average value of 765.57 tons/km 2 (54.97%), while the lowest sediment export came from forest land, with a value of 0.56 tons/km 2 (0.04%). Additionally, urban and built-up areas, water bodies, and wetlands do not create sediment export, according to the C and P coefficients in the biophysical table.

Sediment Export Estimation of Predicted LULC under Scenario I
Estimates of total and average sediment export under the predicted LULC between 2020 and 2029 under Scenario I (Trend of LULC evolution) and mean rainfall erosivity are presented in Table 21. As a result, the highest total and average sediment exports were about 58,798 tons and 65.97 tons/km 2 , respectively, in 2026. The lowest total and average sediment exports were about 32,373 tons and 36.32 tons/km 2 , respectively, in 2022. These results indicate that the primary influence of rainfall erosivity and LULC types in the RUSLE model affect sediment export, as has been suggested by many researchers [56][57][58][59][60]. The influence of rainfall erosivity on sediment export in this study was confirmed by simple linear regression analysis, resulting in Figure 7, with an R 2 value of 0.6424. This finding indicates that rainfall erosivity can explain the linear relationship with the sediment export by about 64%.

Sediment Export Estimation of Predicted LULC under Scenario II
The estimated total and average sediment export under predicted LULC between 2020 and 2029 and mean rainfall erosivity, according to Scenario II (Maximization of ecosystem service values), are presented in Table 22. As a result, the highest total and average sediment exports were about 48,115 tons and 53.98 tons/km 2 , respectively, in 2026. The lowest total and average sediment exports were about 30,190 tons and 33.87 tons/km 2 , respectively, in 2022.
As with Scenario I, these results indicated that the primary influence of the rainfall erosivity factor and LULC types in the RUSLE model affected sediment export. The influence of rainfall erosivity on sediment export under Scenario II was relatively stronger than that in Scenario I, as shown in Figure 8. According to the R 2 value of 0.9589, the rainfall erosivity can explain the linear relationship with sediment export by about 96%.
Moreover, the contribution of the predicted LULC under Scenario II to sediment export between 2020 and 2029 revealed that miscellaneous land decreased every year under this scenario; still, it caused the highest average sediment export until 2027, with values between 1088.47 tons/km 2 in 2020 and 512.54 tons/km 2 in 2027. After that, the field crop yielded the highest average sediment export in 2028 and 2029, with values between 447.46 tons/km 2 and 470.16 tons/km 2 , respectively. In contrast, forest land generated the lowest average sediment export with values between 0.64 tons/km 2 in 2022 and 1.02 tons/km 2 in 2026.

Sediment Export Estimation of Predicted LULC under Scenario II
The estimated total and average sediment export under predicted LULC between 2020 and 2029 and mean rainfall erosivity, according to Scenario II (Maximization of ecosystem service values), are presented in Table 22. As a result, the highest total and average sediment exports were about 48,115 tons and 53.98 tons/km 2 , respectively, in 2026. The lowest total and average sediment exports were about 30,190 tons and 33.87 tons/km 2 , respectively, in 2022.
As with Scenario I, these results indicated that the primary influence of the rainfall erosivity factor and LULC types in the RUSLE model affected sediment export. The influence of rainfall erosivity on sediment export under Scenario II was relatively stronger than that in Scenario I, as shown in Figure 8. According to the R 2 value of 0.9589, the rainfall erosivity can explain the linear relationship with sediment export by about 96%.
Moreover, the contribution of the predicted LULC under Scenario II to sediment export between 2020 and 2029 revealed that miscellaneous land decreased every year under this scenario; still, it caused the highest average sediment export until 2027, with values between 1088.47 tons/km 2 in 2020 and 512.54 tons/km 2 in 2027. After that, the field crop yielded the highest average sediment export in 2028 and 2029, with values between 447.46 tons/km 2 and 470.16 tons/km 2 , respectively. In contrast, forest land generated the lowest average sediment export with values between 0.64 tons/km 2 in 2022 and 1.02 tons/km 2 in 2026.

Sediment Export Estimation of Predicted LULC under Scenario III
The estimated total and average sediment export under the predicted LULC between 2020 and 2029 and mean rainfall erosivity according to Scenario III (Economic crop zonation) are presented in Table 23. As a result, the highest total and average sediment exports were about 76,068 tons and 85.34 tons/km 2 , respectively, in 2026. The lowest total and average sediment exports are about 34,520 tons and 38.73 tons/km 2 , respectively, in 2022.
As with Scenario I, these results indicate that the primary influence of rainfall erosivity factor and LULC types in the RUSLE model affected sediment export; however, the influence of rainfall erosivity on sediment export under Scenario III was lower than that under Scenario I, as shown in Figure 9. The simple linear equation indicates a moderately positive correlation between the mean rainfall erosivity factor and average sediment export under Scenario III, with an R value of 0.5235 and R 2 value of 0.274. According to the R 2 value, the rainfall erosivity can explain the linear relationship with sediment export by only 27%.
The contribution of the predicted LULC under Scenario III to sediment export between 2020 and 2029 reveals that miscellaneous land caused the highest average sediment export, with values between 1658.33 tons/km 2 in 2022 and 7372.09 tons/km 2 in 2029, while forest land generated the lowest average sediment export, with values between 0.65 tons/km 2 in 2022 and 1.08 tons/km 2 in 2026.

Sediment Export Estimation of Predicted LULC under Scenario III
The estimated total and average sediment export under the predicted LULC between 2020 and 2029 and mean rainfall erosivity according to Scenario III (Economic crop zonation) are presented in Table 23. As a result, the highest total and average sediment exports were about 76,068 tons and 85.34 tons/km 2 , respectively, in 2026. The lowest total and average sediment exports are about 34,520 tons and 38.73 tons/km 2 , respectively, in 2022.
As with Scenario I, these results indicate that the primary influence of rainfall erosivity factor and LULC types in the RUSLE model affected sediment export; however, the influence of rainfall erosivity on sediment export under Scenario III was lower than that under Scenario I, as shown in Figure 9. The simple linear equation indicates a moderately positive correlation between the mean rainfall erosivity factor and average sediment export under Scenario III, with an R value of 0.5235 and R 2 value of 0.274. According to the R 2 value, the rainfall erosivity can explain the linear relationship with sediment export by only 27%.  The average sediment export between 2020 and 2029 under the three scenarios then compared, as shown in Figure 10, which shows that the predicted LULC under nario II (Maximization of ecosystem service values) delivered the lowest annual sedi export, compared to Scenarios I and III, with an average value of 42.27 tons/km 2 d the increasing areas of wetland and decreasing areas of rangeland and miscellaneous In contrast, the areas of miscellaneous land in scenarios I and III were increased, base the annual change rate of the Markov chain model. Though miscellaneous land sho only a minor increase, it caused high soil loss and sediment export.
Moreover, Scenario III delivered the highest annual sediment export when comp to the other scenarios, as the paddy field and field crop areas increased, according to nomic crop zonation by the LDD in 2018. The contribution of the predicted LULC under Scenario III to sediment export between 2020 and 2029 reveals that miscellaneous land caused the highest average sediment export, with values between 1658.33 tons/km 2 in 2022 and 7372.09 tons/km 2 in 2029, while forest land generated the lowest average sediment export, with values between 0.65 tons/km 2 in 2022 and 1.08 tons/km 2 in 2026.
The average sediment export between 2020 and 2029 under the three scenarios was then compared, as shown in Figure 10, which shows that the predicted LULC under Scenario II (Maximization of ecosystem service values) delivered the lowest annual sediment export, compared to Scenarios I and III, with an average value of 42.27 tons/km 2 due to the increasing areas of wetland and decreasing areas of rangeland and miscellaneous land. In contrast, the areas of miscellaneous land in scenarios I and III were increased, based on the annual change rate of the Markov chain model. Though miscellaneous land showed only a minor increase, it caused high soil loss and sediment export.
Moreover, Scenario III delivered the highest annual sediment export when compared to the other scenarios, as the paddy field and field crop areas increased, according to economic crop zonation by the LDD in 2018.
the increasing areas of wetland and decreasing areas of rangeland and miscellaneous land. In contrast, the areas of miscellaneous land in scenarios I and III were increased, based on the annual change rate of the Markov chain model. Though miscellaneous land showed only a minor increase, it caused high soil loss and sediment export.
Moreover, Scenario III delivered the highest annual sediment export when compared to the other scenarios, as the paddy field and field crop areas increased, according to economic crop zonation by the LDD in 2018.

Nutrient Export Estimation Using NDR Model
Similar to sediment export, nutrient (N and P) exports were selected to describe ecosystem services in the study area by each LULC allocation scenario (Scenarios I-III).
The statistical performance of the NDR model for nitrogen export estimation under the calibration and validation periods showed satisfactory and very good fitting between the observed and estimated nitrogen export, with R 2 values of 0.575 and 0.895, respectively. At the same time, the PBIAS values under the calibration period provided a very good fit performance rate of the model, with a value of −20.42%, and PBIAS values under the validation period provided a good fit performance rate, with a value of 33.39%, in accordance with [51,52].

Nutrient Export Estimation Using NDR Model
Similar to sediment export, nutrient (N and P) exports were selected to describe ecosystem services in the study area by each LULC allocation scenario (Scenarios I-III).
The statistical performance of the NDR model for nitrogen export estimation under the calibration and validation periods showed satisfactory and very good fitting between the observed and estimated nitrogen export, with R 2 values of 0.575 and 0.895, respectively. At the same time, the PBIAS values under the calibration period provided a very good fit performance rate of the model, with a value of −20.42%, and PBIAS values under the validation period provided a good fit performance rate, with a value of 33.39%, in accordance with [51,52].
Meanwhile, the statistical NDR model performance for calibration and validation of phosphorus export provided a very good and good fit between the observed and estimated phosphorus export, with R 2 values of 0.828 and 0.643, respectively. At the same time, the PBIAS values for calibration and validation provided very good and good fitting performance rates of the model, with values of 12.57 and 30.21%, respectively, in accordance with [51,52]. The optimum local model parameters and domain values under the calibration period of the NDR model are summarized in Tables 24 and 25.  The estimated total and average nutrient (N and P) load and export of actual LULC in 2019 are presented in Table 26. Meanwhile, the amount of nutrient export from each LULC type in 2019 is presented in Table 27.
As a result, the total nitrogen and phosphorus loads in 2019 were about 1,422,800 kg (1596.23 kg/km 2 ) and about 308,268 kg (345.84 kg/km 2 ), and the total nitrogen and phosphorus exports were about 193,308 kg (216.87 kg/km 2 ) and about 41,979 kg (47.10 kg/km 2 ), respectively.  In the meantime, the highest total nitrogen and phosphorus exports occurred at the paddy field, with values of about 95,829 and 26,135 kg, respectively. In contrast, the lowest total nitrogen export came from waterbody areas, with a value of 2.22 kg, and the lowest total phosphorus export came from miscellaneous land, with a value of 0.14 kg. However, the highest average nitrogen and phosphorus export appeared on field crop and perennial trees and orchards, with values of 454.14 and 127.83 kg/km 2 , respectively.

Nutrient Export Estimation of Predicted LULC under Scenario I
The estimated total and average nitrogen and phosphorus export of predicted LULC between 2020 and 2029 and mean annual rainfall under Scenario I (Trend of LULC evolution) is presented in Table 28. As a result, the highest total and average nitrogen export occurred in 2029, with values of about 210,908 kg and 236.62 kg/km 2 , respectively, while the lowest total and average nitrogen export occurred in 2020 with a value of about 197,973 kg and 222.10 kg/km 2 , respectively. Likewise, the highest total and average phosphorus export occurred in 2029, with a value of about 47,221 kg and 52.98 kg/km 2 , respectively, while the lowest total and average phosphorus export occurred in 2020, of about 43,358 kg and 48.64 kg/km 2 , respectively.
These results indicate that the changes in LULC types and areas affect parameters in the biophysical table, leading to different nitrogen and phosphorus export. In contrast, the annual rainfall, as a runoff proxy, was not observed as being sensitive to the estimated data, due to its calculation to modify the load, in order to account for runoff potential by relating the precipitation per cell to the average over the raster, as suggested by [61]. Therefore, non-linear regression analysis was considered in order to reconfirm the suggestion of [61] and the previous study of [62]. Figure 11 shows the results of non-linear regression analysis between mean annual rainfall (mm) and average nutrient export (kg per km 2 ). The best fit of the non-linear regression by Trend Analysis in the MS Excel software was a sixth-order polynomial equation with R 2 values of 0.9638 for nitrogen export and 0.9615 for phosphorus export. As a result, it was confirmed that annual rainfall, as a runoff proxy, is not sensitive to the estimated nutrient export under the NDR model. Furthermore, the contribution of the predicted LULC under Scenario I on nutrient export indicated that the highest total nitrogen and phosphorus export occurred on paddy fields. In contrast, the lowest total nitrogen export occurred on water bodies, and the lowest total phosphorus export occurred on miscellaneous land; however, the highest average nitrogen and phosphorus export appeared on field crops, as well as perennial trees and orchards. gestion of [61] and the previous study of [62]. Figure 11 shows the results of non-linear regression analysis between mean annual rainfall (mm) and average nutrient export (kg per km 2 ). The best fit of the non-linear regression by Trend Analysis in the MS Excel software was a sixth-order polynomial equation with R 2 values of 0.9638 for nitrogen export and 0.9615 for phosphorus export. As a result, it was confirmed that annual rainfall, as a runoff proxy, is not sensitive to the estimated nutrient export under the NDR model. Furthermore, the contribution of the predicted LULC under Scenario I on nutrient export indicated that the highest total nitrogen and phosphorus export occurred on paddy fields. In contrast, the lowest total nitrogen export occurred on water bodies, and the lowest total phosphorus export occurred on miscellaneous land; however, the highest average nitrogen and phosphorus export appeared on field crops, as well as perennial trees and orchards.
These findings suggest that the change in LULC types associated with the biophysical table parameters affects nitrogen and phosphorus export. The LULC data of Scenario I were simulated based on the annual rate of LULC change from the transition area matrix between 2009 and 2019 using the Markov Chain model, which did not represent dramatic change under this scenario; the minor change in area also changed the load amounts and export of nutrients.

Nutrient Export Estimation of Predicted LULC under Scenario II
The estimated total and average nitrogen and phosphorus export of predicted LULC between 2020 and 2029 and mean annual rainfall under Scenario II (Maximization of ecosystem service values) is presented in Table 29.
As a result, the highest total and average nitrogen export occurred in 2020, with values of about 196,965 kg and 220.97 kg/km 2 , respectively, while the lowest total and average These findings suggest that the change in LULC types associated with the biophysical table parameters affects nitrogen and phosphorus export. The LULC data of Scenario I were simulated based on the annual rate of LULC change from the transition area matrix between 2009 and 2019 using the Markov Chain model, which did not represent dramatic change under this scenario; the minor change in area also changed the load amounts and export of nutrients.

Nutrient Export Estimation of Predicted LULC under Scenario II
The estimated total and average nitrogen and phosphorus export of predicted LULC between 2020 and 2029 and mean annual rainfall under Scenario II (Maximization of ecosystem service values) is presented in Table 29. As a result, the highest total and average nitrogen export occurred in 2020, with values of about 196,965 kg and 220.97 kg/km 2 , respectively, while the lowest total and average nitrogen export occurred in 2022, with values of about 195,426 kg and 219.25 kg/km 2 , respectively. Likewise, the highest total and average phosphorus export were in 2029, with values of about 43,798 kg and 49.14 kg/km 2 , respectively, while the lowest total and average phosphorus export occurred in 2022, with values of about 42,961 kg and 48.20 kg/km 2 , respectively. These results indicate that the changes in LULC types and areas affect parameters in the biophysical table, leading to variations in nitrogen and phosphorus export.
As with Scenario I, the annual rainfall, as a runoff proxy, was not a sensitive factor for estimating nutrients. The best fit of the sixth-order polynomial equation between annual rainfall and nutrient export had R 2 values of 0.4920 for nitrogen export and 0.8834 for phosphorus export (Figure 12). As with Scenario I, the annual rainfall, as a runoff proxy, was not a sensitive factor for estimating nutrients. The best fit of the sixth-order polynomial equation between annual rainfall and nutrient export had R 2 values of 0.4920 for nitrogen export and 0.8834 for phosphorus export (Figure 12).
Moreover, the contribution of the predicted LULC under Scenario II on nutrient export demonstrated that the highest total nitrogen and phosphorus export occurred on the paddy fields. In contrast, the lowest nitrogen export occurred on water bodies, while the lowest phosphorus export occurred on miscellaneous land; however, the highest average nitrogen export appeared on paddy fields and field crops, while the highest average phosphorus export occurred on perennial trees and orchards. Additionally, these findings suggest that the change in LULC types associated with the biophysical table parameters affects nitrogen and phosphorus export. Notably, the LULC data under this scenario influenced nutrient export due to the LULC data of Scenario II, simulated based on the annual rate of LULC change from transition area matrix between 2009 and 2019 for some LULC types and the LP to maximize ecosystem service Moreover, the contribution of the predicted LULC under Scenario II on nutrient export demonstrated that the highest total nitrogen and phosphorus export occurred on the paddy fields. In contrast, the lowest nitrogen export occurred on water bodies, while the lowest phosphorus export occurred on miscellaneous land; however, the highest average nitrogen export appeared on paddy fields and field crops, while the highest average phosphorus export occurred on perennial trees and orchards.
Additionally, these findings suggest that the change in LULC types associated with the biophysical table parameters affects nitrogen and phosphorus export. Notably, the LULC data under this scenario influenced nutrient export due to the LULC data of Scenario II, simulated based on the annual rate of LULC change from transition area matrix between 2009 and 2019 for some LULC types and the LP to maximize ecosystem service values, by reducing the area of rangeland and miscellaneous land and increasing the area of wetland. Therefore, significant LULC change was observed under this scenario, due to the area change, load amounts, and export of nutrients.

Nutrient Export Estimation of Predicted LULC under Scenario III
The estimated total and average nitrogen and phosphorus export according to the predicted LULC between 2020 and 2029 and mean annual rainfall under Scenario III (Economic crop zonation) are presented in Table 30.
As a result, the highest total and average nitrogen export under Scenario III were about 229,756 kg and 257.76 kg/km 2 , respectively, occurring in 2029, while the lowest total and average nitrogen export were about 200,387 kg and 224.81 kg/km 2 , respectively, occurring in 2020. Likewise, the highest total and average phosphorus export were about 51,149 kg and 57.38 tons/km 2 , respectively, occurring in 2029, while the lowest total and average phosphorus export were about 43,956 kg and 49.31 kg/km 2 , respectively, occurring in 2020. These results indicate that the changes in LULC types and areas affect parameters in the biophysical table, which leads to different nitrogen and phosphorus export, as observed in Scenarios I and II. The best fit of the sixth-order polynomial equation between annual rainfall and nutrient export had R 2 values of 0.9530 for nitrogen export and 0.9469 for phosphorus export (Figure 13). As a result, the highest total and average nitrogen export under Scenario III were about 229,756 kg and 257.76 kg/km 2 , respectively, occurring in 2029, while the lowest total and average nitrogen export were about 200,387 kg and 224.81 kg/km 2 , respectively, occurring in 2020. Likewise, the highest total and average phosphorus export were about 51,149 kg and 57.38 tons/km 2 , respectively, occurring in 2029, while the lowest total and average phosphorus export were about 43,956 kg and 49.31 kg/km 2 , respectively, occurring in 2020. These results indicate that the changes in LULC types and areas affect parameters in the biophysical table, which leads to different nitrogen and phosphorus export, as observed in Scenarios I and II. The best fit of the sixth-order polynomial equation between annual rainfall and nutrient export had R 2 values of 0.9530 for nitrogen export and 0.9469 for phosphorus export (Figure 13).  Furthermore, the contribution of the predicted LULC under Scenario III on nutrient export indicated that the highest total nutrient and phosphorus export occurred on paddy fields. In contrast, the lowest total nitrogen export occurred on water bodies, while the lowest total phosphorus export occurred on miscellaneous land. Meanwhile, the highest average nitrogen and phosphorus export occurred on field crops, as well as perennial trees and orchards.
These findings suggest that the change in LULC types associated with the biophysical table parameters affects nitrogen and phosphorus export. In particular, the LULC data under this scenario influenced nutrient export, according to the LULC data of Scenario III simulated based on the annual rate of LULC change from transition area matrix between 2009 and 2019 for some LULC types and the economic crop zonation, particularly the increase in paddy field and field crop areas and the decrease in perennial tree and orchard areas, which represented dramatic change under this scenario, leading to changes in the load amounts and export of nutrients. Hence, significant LULC change was observed under this scenario.
The average nutrient (N and P) export between 2020 and 2029 for the three scenarios was compared, as shown in Figure 14, which shows that the predicted LULC under Scenario II (Maximization of ecosystem service values) delivered the lowest annual nitrogen and phosphorus export, compared to Scenarios I and III, due to the allocation of area using LP to maximize the ecosystem service values by increasing wetland areas. Increases in wetland areas, such as wetland restoration and constructed riparian wetlands, can reduce nitrogen [63] and phosphorus [64] export into water bodies. Furthermore, based on the biophysical table, wetland provides the highest maximum retention efficiency and provides low nitrogen and phosphorus load. The decreased paddy field area under scenario II affected the nutrient export, as this LULC type supplies the highest nitrogen and phosphorus load with low maximum retention efficiency.
Nevertheless, scenario III (Economic crop zonation) generated higher nutrient export than other scenarios, as the paddy field and field crop areas were increased according to their LDD suitability classes. These areas provide the highest nitrogen and phosphorus load, but low maximum retention efficiency. This result is in agreement with [65], who applied the NDR model to calculate nutrient export under two different scenarios. They found that the cropland balance policy negatively impacted water purification by increasing nitrogen export, which was 8.36% higher than that in the no strict cropland protection scenario.
The average nutrient (N and P) export between 2020 and 2029 for the three scenarios was compared, as shown in Figure 14, which shows that the predicted LULC under Scenario II (Maximization of ecosystem service values) delivered the lowest annual nitrogen and phosphorus export, compared to Scenarios I and III, due to the allocation of area using LP to maximize the ecosystem service values by increasing wetland areas. Increases in wetland areas, such as wetland restoration and constructed riparian wetlands, can reduce nitrogen [63] and phosphorus [64] export into water bodies. Furthermore, based on the biophysical table, wetland provides the highest maximum retention efficiency and provides low nitrogen and phosphorus load. The decreased paddy field area under scenario II affected the nutrient export, as this LULC type supplies the highest nitrogen and phosphorus load with low maximum retention efficiency.
Nevertheless, scenario III (Economic crop zonation) generated higher nutrient export than other scenarios, as the paddy field and field crop areas were increased according to their LDD suitability classes. These areas provide the highest nitrogen and phosphorus load, but low maximum retention efficiency. This result is in agreement with [65], who applied the NDR model to calculate nutrient export under two different scenarios. They found that the cropland balance policy negatively impacted water purification by increasing nitrogen export, which was 8.36% higher than that in the no strict cropland protection scenario.

Suitable LULC Allocation Scenario to Minimize Sediment Export
The ESCI values of sediment export in ten periods-as well as its average under the three scenarios-were compared, in order to identify a suitable LULC allocation scenario for minimizing sediment export, in terms of ecosystem service change, as summarized in Table 31 and Figure 15.

Suitable LULC Allocation Scenario to Minimize Sediment Export
The ESCI values of sediment export in ten periods-as well as its average under the three scenarios-were compared, in order to identify a suitable LULC allocation scenario for minimizing sediment export, in terms of ecosystem service change, as summarized in Table 31 and Figure 15. sediment export of this scenario were also the lowest, with an average ESCI value of 0.4260. Therefore, the LULC allocation of Scenario II was chosen, in order to minimize the sediment export into Kwan Phayao from the Upper Ing watershed. Moreover, the average ESCI for sediment export under the three allocation LULC scenarios was tested, in terms of the difference of the mean, using the t-test statistic. The results demonstrated significant differences among average ESCI values on sediment export of three scenarios at the 95% confidence level. See detail in Table A8 in Appendix B. Figure 15. Comparison of ESCI on sediment export in ten periods under three scenarios.

Suitable LULC Allocation Scenario to Minimize Nutrient Export
The ESCI values for nutrient (N and P) export in ten periods and its average of three scenarios were compared, in order to identify a suitable LULC allocation scenario for minimizing nutrient export, in terms of ecosystem service change, are reported in Tables 32  and 33 and shown in Figures 16 and 17. According to the result, the LULC allocation of Scenario II (Maximization of ecosystem service values) generated the lowest sediment export every year between 2020 and 2029, with an average sediment export of 37,677.19 tons. The cumulative ESCI values on sediment export of this scenario were also the lowest, with an average ESCI value of 0.4260. Therefore, the LULC allocation of Scenario II was chosen, in order to minimize the sediment export into Kwan Phayao from the Upper Ing watershed.
Moreover, the average ESCI for sediment export under the three allocation LULC scenarios was tested, in terms of the difference of the mean, using the t-test statistic. The results demonstrated significant differences among average ESCI values on sediment export of three scenarios at the 95% confidence level. See detail in Table A8 in Appendix B.

Suitable LULC Allocation Scenario to Minimize Nutrient Export
The ESCI values for nutrient (N and P) export in ten periods and its average of three scenarios were compared, in order to identify a suitable LULC allocation scenario for minimizing nutrient export, in terms of ecosystem service change, are reported in Tables 32 and 33   Note: * The average value from data between 2020 and 2029.   Moreover, the average ESCI on nitrogen and phosphorus export of three LULC allo- Moreover, the average ESCI on nitrogen and phosphorus export of three LULC allocation scenarios was used to test the mean differences using the t-test statistic. The results revealed significant differences among average ESCI values for nutrient export under the three scenarios at the 95% confidence level. See detail in Tables A9 and A10 in Appendix B.

Suitable LULC Allocation Scenario to Minimize Sediment and Nutrient Export
The average ESCI values for sediment and nutrient (N and P) exports in ten periods of three LULC allocation scenarios were compared, in order to identify suitable LULC allocation scenarios to minimize sediment and nutrient export, in terms of ecosystem service change, as shown in Figure 18 and summarized in Table 34.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 33 Figure 18. Comparison of average ESCI value of sediment and nutrient export on ecosystem ser among the three scenarios.    Therefore, the LULC allocation of Scenario II (Maximization of ecosystem service values) was selected as the suitable LULC allocation scenario to minimize sediment and nutrient exports into Kwan Phayao from the Upper Ing watershed. These findings can serve as crucial information to allocate LULC in the Upper Ing watershed by land-use planners, land managers, and decision-makers in order to minimize sediment and nutrient loads into Kwan Phayao in the future.

LULC Classification Using SVM Algorithm
The overall accuracy and Kappa hat coefficient for the thematic accuracy of the LULC map in 2009 and 2019 were 90.86 and 87.00%, and 89.59 and 85.85%, respectively. Kappa hat coefficient values of more than 80% represent a substantial agreement between the classification map and the reference data [66]; likewise, the overall accuracy of the LULC maps-higher than 85%-indicates that they can provide an acceptable result [67].
In addition, the overall accuracy and Kappa hat coefficient of the current study were consistent with other researchers who have classified LULC data based on Landsat imagery using the SVM algorithm [68][69][70][71]; however, to apply SVM for LULC classification, users must select sample points between the boundaries of LULC classes precisely, where mixed pixels are common, in order to ensure accurate classification [72]. Therefore, selecting the SVM algorithm's training points for LULC classification under the EnMap BOX software requires time and skill.

Land Requirement Estimation under Three Different Scenarios
Land requirement is essential information for LULC prediction using the CLUE-S model. In this study, three different scenarios including (1) Scenario I (Trend of LULC evolution), (2) Scenario II (Maximization of ecosystem service values), and Scenario III (Economic crop zonation) are estimated based on their characteristics. The land requirement of Scenario I is estimated based on the annual change rate of LULC between 2009 and 2019 from the transition area matrix using the Markov Chain model. The land requirement of Scenario I is dictated by the accuracy of LULC data in 2009 and 2019. Meanwhile, the land requirement of Scenario II is estimated based on LULC allocation for maximizing ESV using Linear Programming with the Simplex method. The land requirement of Scenario II depends on the efficiency of linear programming for determining reclaimed and other areas changing into the wetland. At the same time, the land requirement of Scenario III is estimated based on areas of suitability classes for economic crops (paddy field, field crop, para rubber, and perennial trees and orchards) and the Markov Chain model. The land requirement of Scenario III is delimited by the economic crop zonation and the accuracy of LULC data in 2019. Consequently, the definition of each Scenario should first be well defined and local government agencies and stakeholders should be consulted before land requirement estimation.

LULC Prediction by CLUE-S Model
According to binary logistics regression analysis, the derived AUC values for each LULC type allocation under the CLUE-S model exhibit good fit (0.80) and excellent fit (0.99) between the predicted and real LULC transition, as mentioned by Chen et al. [73]. These results imply that the significant driving factors on LULC change, including soil drainage, distance to stream, distance to water body, distance to village, slope, distance to road, annual rainfall, elevation, and population density at subdistrict level, are suitable to apply for each LULC type allocation under CLUE-S model.
As a result of LULC prediction under three scenarios by the CLUE-S model, the predicted LULC data of three scenarios can deliver realistic results as an expectation. The deviation values between the required land area and the predicted area of each LULC type under three scenarios are very small and vary from −0.0023 to 0.0033% or −0.23 km 2 (underestimation) to 0.33 km 2 (overestimation). In principle, the deviation value depends on iteration driving factors of each LULC type, which indicate the maximum different allowance between the land requirement and land allocation of LULC type under the CLUE-S model [53,54].
Consequently, the CLUE-S model can be an effective tool to predict LULC data according to specific policies as the scenario. Essentially, the suitable multiple linear equations from the logistics regression analysis for each LULC type allocation, a land requirement of different scenarios assigned by policy transformation, and model parameters (elasticity and LULC conversion matrix) are very important for predicting LULC under the CLUE-S model.

Sediment Export Estimation
According to the results of sediment export estimation, it was found that the primary influence of rainfall erosivity and LULC types in the RUSLE model affect sediment export under three different scenarios. The rainfall erosivity can explain the linear relationship with the sediment export from about 27% in Scenario III (Economic crop zonation) to 96% in Scenario II (Maximization of ecosystem service values). These findings are consistent with the previous studies [56][57][58][59][60].
The contribution of the predicted LULC under three scenarios on sediment export between 2020 and 2029 indicated that miscellaneous land causes the highest average sediment export. Meanwhile, forest land generates the lowest average sediment export. This finding was consistent with the previous study of Srichaichana et al. [14]. They found that miscellaneous land (bare land and abandoned mine) created the highest average sediment export, with a value of 659.72 tons/km 2 . At the same time, evergreen forest generated the lowest average sediment export, with a value of 0.01 tons/km 2 , in the Klong U-Tapao watershed, Songkhla Province, Thailand. Likewise, Degife et al. [74] found that the highest sediment export per unit of area was observed from miscellaneous land (bare land), while the highest contribution of the total sediment that reached the surrounding water bodies was from cultivated land (40.7%). Similarly, Zhou et al. [75] found that decreases in miscellaneous land (bare land) significantly reduced sediment export in the Qiantang River Basin, China. In contrast, increases in agricultural land, such as cropland and garden plots, increased sediment export in the studied watershed.

Nutrient Export Estimation
The nutrient (N and P) export estimation under three different scenarios indicates that the change in LULC types and areas affects parameters in the biophysical table, leading to different nitrogen and phosphorus export. The contribution of LULC type under three scenarios indicates that the highest total nitrogen and phosphorus exports occur on paddy fields. In contrast, the lowest total nitrogen export occurs on water bodies, and the lowest total phosphorus export occurs on miscellaneous land. However, the highest average nitrogen and phosphorus export mostly appeared on field crops and perennial trees/orchards. These findings agree with [76], who indicated that the highest nitrogen and phosphorus exports occurred on cultivated land. There was very little nitrogen and phosphorus export on forest land, water areas, and unused land. This is similar to the result of [77], who analyzed nutrient load and delivery from different scenarios and found the most significant load rate per unit area and low retention efficiency of cultivated crops. Furthermore, agriculture was the leading cause of nutrient release in the watershed under the different LULC scenarios. These findings are similar to those of [78], who found that cropland and agroforestry influenced roughly 90% of the nutrients exported, while water bodies were identified as sinks.
Moreover, as a runoff proxy, the annual rainfall is not a sensitive factor for estimating nutrients, as suggested by [61] and the previous study of [62]. The best fit of the sixth-order polynomial equation between annual rainfall and nutrient (N and P) export under three different scenarios in this study reconfirm the relationship as mentioned earlier. In this study, the R 2 value of the relations between annual rainfall and nitrogen export varies from 0.4920 to 0.9638, while the R 2 value of the relations between annual rainfall and phosphorus export diverges from 0.8834 to 0.9615.

Conclusions
Land-use and land-cover (LULC) classification and change detection between 2009 and 2019 was successfully conducted using a supervised classification with a support vector machine and post-classification comparison change detection algorithms. The overall accuracy and Kappa hat coefficient of the LULC maps in 2009 and 2019 were higher than 85%. Then, the land requirements under three scenarios-Scenario I (Trend of LULC evolution), Scenario II (Maximization of ecosystem service values), and Scenario III (Economic crop zonation)-were estimated based on their characteristics. Time-series LULC data between 2020 and 2029 were then effectively predicted using the CLUE-S model. After that, the actual LULC data in 2019, as base data, and the predicted LULC data under the three scenarios between 2020 and 2029 were used as significant inputs for ecosystem service assessment, in terms of sediment and nutrient export, using the SDR and NDR models. Finally, a suitable LULC allocation scenario was successfully identified in order to minimize sediment and nutrient export using the ecosystem services change index: the most suitable LULC allocation scenario to minimize sediment or/and nutrient export into Kwan Phayao was Scenario II (Maximization of ecosystem service value).
In conclusion, the integration of remote sensing data with an advanced classification method (support vector machine classifier), GIS data with linear programming, and advanced geospatial models (CLUE-S model, SDR and NDR model) can be used as an efficient tool to assess ecosystem services at the watershed level-particularly sediment and nutrient (N and P) export-and can be further applied to identify a suitable LULC allocation scenario to minimize sediment and nutrient export into certain lakes. The results of the current study can serve as crucial information for land-use planners, land managers, and decision-makers in order to reduce sediment and nutrient export into Kwan Phayao in the future.

Acknowledgments:
The authors would like to thank the Suranaree University of Technology for providing a scholarship to Jiraporn Kulsoontornrat and the University of Phayao for supporting the facilities to undertake this research. Special thanks from the authors go to the anonymous reviewers for their valuable comments and suggestions that improved our manuscript from various perspectives.

Conflicts of Interest:
The authors declare no conflict of interest.