Modeling the Impact of Land Use Optimization on Non-Point Source Pollution: Evidence from Chinese Reservoir Watershed

: We assessed the effectiveness of land optimization for controlling non-point source (NPS) pollution by combining a multi-objective dynamic planning approach with the application of the Land Use and its Effects at Small regional extents model and the Soil and Water Assessment Tool. The combined modeling approach showed substantial ability to reduce NPS pollution in Shitoukoumen Reservoir, Changchun City, China, reducing the annual total loads of nitrogen and phosphorus in the study area by 8.7 and 10.12%, respectively. The total nitrogen load decreased signiﬁcantly and stabilized at less than 8 kg/hm 2 from a peak level of over 15 kg/hm 2 . Higher total phosphorus loads before land use optimization were concentrated in the central parts of the study area, with the highest values exceeding 2.3 kg/hm 2 , and tended to spread outward but resolved at 1.5 kg/hm 2 after optimization. The results showed that from a macro-perspective, optimization of the spatial distribution and quantitative composition of land use can effectively control NPS pollution. The study also demonstrates the potential effectiveness of the coupled multi-model methodology for mitigating NPS in the future.


Introduction
Our survival relies on water, but shortages (especially of clean water) are of increasing concern.Shortages of clean water are greatly exacerbated by two types of water pollution: point source (from fixed discharge points) and non-point source [1].Non-point source (hereafter NPS) pollution is more difficult to control, mainly because of its diversity, dispersal, and complex inter-relationships [2,3].In China, NPS water pollution is becoming increasingly acute, with no clear sign of improvement [1,4]), so there are increasingly urgent needs to characterize it and develop effective management policies to control it and protect the environment.
Experimental observations provide the most direct data for identifying the main drivers of NPS pollution, but the feasibility of empirical studies is limited by its extensive spatiotemporal distribution.Thus, the use of hydrological models is generally required for studies of such pollution and emissions at watershed (or larger) scales.
Factors influencing NPS pollution can be grouped into four types: topographical, climatic, human agency, and land use [5].Land use practices are key to addressing it as they may influence the other three types.The spatial distribution of these practices affects soil quality, as various land use types change the soil's hydrological characteristics and both the directions and rates of runoff, thereby changing the material cycling processes and energy distribution within the soil.Hence, as shown by many studies, land uses can have profound and extensive impacts on NPS pollution [6][7][8].Both the types and spatial distributions of land use affect its migration and transformations, so there are strong correlations between NPS pollution loads and land use in small-and medium-sized Land 2024, 13, 18 2 of 17 watersheds [9,10].Irresponsible practices can cause severe soil erosion with consequent losses of nitrogen, phosphorus, and other nutrients, followed by watershed-wide pollution due to the migration of large amounts of these substances into water bodies.Thus, there is an urgent need to determine the effects of different land use practices on NPS pollution [11].
The Soil and Water Assessment Tool (SWAT) is a powerful hydrological model with widely verified simulation accuracy and effectiveness that has been extensively used to simulate sediment and NPS pollution under various land uses and soil conditions [2,3,12,13].However, it is difficult to determine correlations between changes in land use and NPS pollution using any single hydrological or land use change model, and more comprehensive modeling is certainly required to optimize land use for pollution control.Therefore, hydrological models must be combined with land use change models [14].The Land Use and its Effects at Small regional extents (CLUE-S) model has verified efficacy for optimizing land allocation in small watersheds, utility for quantitative spatial analyses of multi-scale land use changes, and ability to simulate future land use changes at small and medium scales rapidly and accurately.Thus, it has been frequently used for simulating future land use changes [1,[15][16][17].The overall goal of land optimization is to ensure that land-based resources are rationally allocated, so any methodology applied must include consideration of the scarcity of land resources and solution of a multi-objective decision-making problem (MODP).The MODP methodology has developed from linear planning to a more mature mathematical approach for solving complex problems related to the evaluation, decisionmaking, and prediction of urban, economic, and social development, which inevitably involves multiple intricately linked objectives [18][19][20].We hypothesized that combining the CLUE-S and SWAT models in a MODP approach could provide further improvement of both the rationality and accuracy of land use optimization for robust NPS evaluation and control of NPS.
As a proof-of-concept, we chose to couple a land use model with a hydrological model within the watershed of the Shitoukoumen Reservoir, which stands as the primary water source in Changchun City [21,22] The key objectives encompassed: Specifically, our land use model incorporated the MODP method, aimed at fostering regional sustainable development.

Study Area
The study focused on the Shitoukoumen Reservoir watershed, which (as already mentioned) is a key source of water in Changchun and has a catchment area of nearly 5 km 2 , covering about two-thirds of the Drinking Horse River watershed (Figure 1).The reservoir is situated in the middle reaches of the Drinking Horse River (43 • 58 N, 125 • 45 E) in the transition zone from tablelands of the Songliao plain to the low hills of the Dahei mountain range (central Jilin Province).The study area's elevation ranges between 188 and 1057 m, gradually decreasing from the southeast to the northwest.Its climate is classified as temperate, continental, and humid, with a dry, windy spring, short but warm summer with more rainfall than later in the year, rapidly cooling temperatures and less rainfall in the autumn, and longer periods of colder temperatures in the winter [23].The average annual temperature is about 5.3 • C, while the annual precipitation amounts to 369.9-667.9mm, with half falling in July and August [24].The average annual evaporation is 1658.1 mm, and more than half occurs in April-June, often leading to drought and severe water shortages [25].According to the Harmonized World Soil Database (HWSD), the soils in the watershed are Haplic Phaeozems, Stagnic Phaeozems, Calcaric Fluvisols, Terric Histosols, Cumulic Anthrosols, Calcic Gleysols, Haplic Luvisols, Albic Luvsiols, and Eutric Cambisols.The areas of arable land, woodland, grassland, water, residential land, and underutilized land are reportedly 247,845 hectares, 186,881.75 hectares, 3884.5 hectares, 16,719.25 hectares, 31,158.25 hectares, and 2766.25 hectares, respectively.
Land 2024, 13, x FOR PEER REVIEW 3 of 18 369.9-667.9mm, with half falling in July and August [24].The average annual evaporation is 1658.1 mm, and more than half occurs in April-June, often leading to drought and severe water shortages [25].

Model Approach
Our land use-based simulation of NPS pollution consisted of three steps.First, we used logistic regression analysis to extract land use distribution patterns, and a MODP approach with the CLUE-S model to optimize the spatial distribution of land use types.Next, we constructed a localized SWAT NPS pollution simulation model and verified the model's parameters and accuracy.Finally, keeping all other factors constant, we used a dynamic land-type input method to simulate NPS pollution before and after optimization, and compared the results.In this manner the NPS pollution in the study area was controlled, and the control method's effectiveness evaluated.The study framework is shown in Figure 2.

Model Approach
Our land use-based simulation of NPS pollution consisted of three steps.First, we used logistic regression analysis to extract land use distribution patterns, and a MODP approach with the CLUE-S model to optimize the spatial distribution of land use types.Next, we constructed a localized SWAT NPS pollution simulation model and verified the model's parameters and accuracy.Finally, keeping all other factors constant, we used a dynamic land-type input method to simulate NPS pollution before and after optimization, and compared the results.In this manner the NPS pollution in the study area was controlled, and the control method's effectiveness evaluated.The study framework is shown in Figure 2.

MODP Method
Land use optimization can be defined as the rational allocation of land cover, e.g., the spatial distribution and structure of land uses in a specified area to maximize realization of a specific objective or multiple objectives under certain constraints in response to prob-

MODP Method
Land use optimization can be defined as the rational allocation of land cover, e.g., the spatial distribution and structure of land uses in a specified area to maximize realization of a specific objective or multiple objectives under certain constraints in response to problems such as conflicts between land supply and demand [26].A MODP approach can be used to express specific objectives quantitatively by establishing objective functions and using constraints to achieve focused and integrated effects.Detailed descriptions of MODP have been published by various authors [27,28].The approach we applied consisted of the following steps [29][30][31]: (1) Selection of a time series, based on existing data, reflecting a specific time interval before and after applying the optimal design of the final model; (2) Determination of multi-objective decision variables, with the range of each class selected as the decision variable; (3) Construction of MODP constraints, based on land use characteristics and optimal allocation requirements, taking into account the needs and operability of a comprehensive, targeted management plan incorporating the benefits (social, economic, and ecological) of land use.
The constraint function was expressed as: where a ij is the coefficient of the constraint function; x j is the area of each land use type; j to n represent the coding of land use types, and b i is the constraint constant.
(4) Construction of a multi-objective decision function for the establishment of a land use optimization model, based on specific objectives and designed to maximize the benefits (economic, social, and/or ecological).The objective function was expressed as: where st. is an abbreviation for subject to, which indicates constraints, f i (x) is the ith objective function value, f i (x) is the pre-determined desired value for each target, while d + i and d − i are positive and negative deviation variables, indicating differences between actual values and values that exceed and fall short of the target values, respectively.
An ecological conservation scenario was set as the optimization objective, with the maximum total ecosystem service value, based on equivalents published by Xie et al. [32], as the optimized objective function.

CLUE-S Model
The basic principle of the CLUE-S model is to extract current land distribution rules and use them to simulate future land use distributions under selected constraints.It includes both non-spatial and spatial allocation modules.The former integrates drivers of land use change within a region to calculate the demand for changes in land use each year over a defined time.As the calculation efficiency is low, the MODP approach was added to this module to improve the measurement of land demand [11,33].
We used data on land use for a baseline period for the spatial allocation module.A logistic regression analysis of the identified drivers was carried out [34], then constraints and land-use conversion rules were constructed.We then derived an optimal land use distribution for the study area and expressed it quantitatively using regression equations [35].
The optimal layout for each land evaluation unit within the study area was determined by calculating the probable spatial distribution of every land use type and its optimal utilization.Thus, the spatial allocation involved a combination of empirical spatial analysis and dynamic simulation [21, 36,37] with the following four parts: (1) Patterns of land distribution: we applied a logistic regression analysis to explore relationships between the drivers of each land type's spatial distribution and its probable distribution, which was quantitatively described using the following regression equation: where P i is the probability of land use type i occurring on a particular raster image element, β n is the explanatory variable's coefficient for each regression equation, and X n,i is the value for each driver of different land types.
(2) Rules of land use transfer: a land transfer matrix and the elasticity of land use type conversion (ELAS) were used to set the rules of conversion between land types [38].
The ELAS parameter has a range of 0-1, and the value is positively correlated to the difficulty of converting a specified use type to other types: the smaller the value, the less stable the use and the easier it is to convert.(3) Calculating land demands: areas of land use can be predicted using the non-spatial module, or other models can be introduced [39].In this study, we used the MODP approach.(4) Spatial allocation of land types: the CLUE-S spatial allocation module is based on the probable distributions of the land use types in a given area, using the formula: where TPROP i,u is the total distribution probability of raster cell i for each land use type u, P i,u is the applicability of raster cell u to land type i derived from logistic regression analysis, ELAS u is the conversion elasticity coefficient for land use type u, and ITER u is the specific iteration coefficient.
More detailed explanations of the principles and operation of the model are provided elsewhere [40].

SWAT Model
The SWAT model, presented by Arnold et al. [41], is a powerful distributed hydrological model that was initially applied to simulate long-term effects of various land management factors, such as land use and soil types, on levels of runoff, sediment, and pollutants within a watershed [42].It has a strong physical component because its predecessor included hydrological, soil erosion, and nutrient components, such as crop growth, pesticide use, and daily precipitation [35,43].SWAT favors the simulation of hydrological processes and can calculate daily time steps [44,45] that simulate the hydrological cycle, soil erosion, sediment transport, and infiltration and transfer of various pollutants in a watershed, thus facilitating comprehensive analysis of water resources and their management.Its core element is a hydrological process model than can simulate various aspects of the hydrological cycle, including water balance, evapotranspiration, soil moisture, groundwater, and runoff-related processes.The minimum computational unit is the "hydrological response unit" (HRU), which is based on the soil type, slope, and land use practices.The model simulates water production, sink flow, and pollutant transport separately for each HRU and finally aggregates the output to the watershed outlet [46].The hydrological process is based on the following water balance equation: Land 2024, 13, 18 6 of 17 where SW t and SW o are the final and initial water contents (mm) of the soil, respectively, t is the calculated timespan (day), R day is precipitation (mm), Q sur f is surface runoff (mm), E a is evaporation (mm), W seep is infiltration and lateral flow (mm), and Q gw is groundwater outflow (mm) [47].

Calibration and Validation
To determine whether a selected driver can improve explanations of a land use distribution pattern, receiver operating characteristic (ROC) curves were used to evaluate the CLUE-S model, and the area under the curve (AUC) statistic was used as a composite indicator to reflect the sensitivity and specificity of continuous variables and test the simulation results' accuracy.An AUC value >0.7 indicates that a selected driver has good explanatory power [26,48].We also used the kappa index to evaluate the model's accuracy, with a threshold value of >0.60 for satisfactory accuracy.The index was calculated as follows: where P o is the proportion of simulated results consistent with realistic elements, P c is the proportion of simulated results consistent with realistic elements under completely random conditions, and P p is the proportion of simulated results consistent with realistic elements under ideal conditions.
To improve the SWAT model's operational efficiency and accuracy, the parameters were adjusted and the model validated using SWAT Calibration and Uncertainty Programs (SWAT-CUP 2012) software [49].Sequential Uncertainty Fitting version 2 (SUFI2) was chosen as the calibration algorithm because at the basin scale it has proven efficacy at a watershed scale [1,50].The correlation coefficient (R 2 ) and Nash coefficient (E NS ) were also used to evaluate the model's effectiveness [51,52].The model met the SWAT model accuracy requirements, with R 2 ≥ 0.6 and E NS ≥ 0.5 [52,53], using the formulae: where Q obs i and Q sim i are the ith measured and simulated values, respectively, while Q obs i and Q sim i are the mean values of all measured and simulated values, respectively [54,55].

Model Input Preparation
The input data for the models (CLUE-S and SWAT) included both spatial information and attribute data.Digital elevation model (DEM) data at 15 m × 15 m resolution, and 30 m × 30 m resolution remote sensing monitoring data on land use/land cover were obtained from the Chinese Academy of Sciences' Resource and Environment Data Center (https://www.resdc.cn/,accessed on 5 November 2023).We derived soil-type data, with 1 km × 1 km precision, from the Harmonized World Soil Database (HWSD) (https:// www.fao.org/,accessed on 5 November 2023).The main sources of meteorological data in the SWAT model were the China Meteorological Assimilation Driving Datasets (CMADS) (http://www.cmads.org/,accessed on 5 November 2023), which are public resources that can significantly improve modeling speed and output accuracy [56].We obtained both runoff and pollutant data from the Changchun Ecological Environment Bureau and road data from the Changchun Natural Resources Bureau.

Model Execution
The CLUE-S model was constructed to maximize the total value of ecosystem services, and land use cover types were optimized and input into the model.The parameters "elevation", "slope", "slope direction", and distances to "railroad and highway", "national and provincial roads", and " county and township roads" were selected, and distances to "national and provincial highways", "county and township roads", " towns", "rural settlements", and "rivers" were used as driving factors to extract land distribution patterns.We evaluated the model's accuracy by comparing the difference between the simulated and real situations in 2020, then used the validated model to optimize the spatial distribution of land uses, and input the results before and after optimization into the SWAT model.
We ran the model by dividing the study area into 46 sub-watersheds (Figure 3) with 1081 HRUs based on slope, land use type, and soil type.We used measured runoff and ammonia pollutant data to determine parameter rates and calibrate the model.To ensure its accuracy, a 1-year warm-up period was set, and simulations were conducted in monthly time steps.As the study focused on effects of land use optimization on the control of NPS pollution, land use type was the variable, climate, soil type, slope, etc., were kept constant, and the optimization effect was evaluated by comparing the difference in NPS pollution before and after optimization.

Model Calibration and Performance Evaluation
We evaluated the CLUE-S model's accuracy by comparing differences between simulated and current land uses in 2020 (Figure 4).The total number of rasters was 1,957,019; the number of rasters where the simulated results matched the empirical data was 1,747,318; and the kappa index was 0.868 (>0.60).All the AUC values for each category, shown in Table 1, were greater than 0.8, indicating that the model had the required simulation accuracy and could be applied in the next stage.

Model Calibration and Performance Evaluation
We evaluated the CLUE-S model's accuracy by comparing differences between simulated and current land uses in 2020 (Figure 4).The total number of rasters was 1,957,019; the number of rasters where the simulated results matched the empirical data was 1,747,318; Land 2024, 13, 18 8 of 17 and the kappa index was 0.868 (>0.60).All the AUC values for each category, shown in Table 1, were greater than 0.8, indicating that the model had the required simulation accuracy and could be applied in the next stage.For the SWAT model, the set calibration period spanned 2012-2016, and the validation period consisted of the years 2017-2018.The R 2 and ENS values for runoff and ammonia nitrogen at each hydrological station are shown in Table 2, and the fitted curves are shown in Figures 5 and 6.The R 2 and ENS values for the calibration and validation periods largely exceeded 0.7, and the fitted curves were good, indicating that the model could be used for robust simulation of NPS pollution in the Shitoukoumen Reservoir's watershed.2, and the fitted curves are shown in Figures 5 and 6.The R 2 and E NS values for the calibration and validation periods largely exceeded 0.7, and the fitted curves were good, indicating that the model could be used for robust simulation of NPS pollution in the Shitoukoumen Reservoir's watershed.

Land Use Optimization
The results of land use optimization using the MODP approach, with ecological protection as the goal, are shown in Table 3. Areas of all land types except woodland and

Land Use Optimization
The results of land use optimization using the MODP approach, with ecological protection as the goal, are shown in Table 3. Areas of all land types except woodland and cropland were reduced after optimization.The areas of construction land and underuti-

Land Use Optimization
The results of land use optimization using the MODP approach, with ecological protection as the goal, are shown in Table 3. Areas of all land types except woodland and cropland were reduced after optimization.The areas of construction land and underutilized land were significantly reduced, while the area of woodland significantly increased, by 8504.85 ha.However, grassland, cropland, and water areas did not change significantly.
We put results of the quantitative optimization into the CLUE-S model to optimize the spatial land use distribution, and the results are shown in Figure 7. Areas of woodland in the northwest and around the Shitoukoumen Reservoir increased significantly, while the more concentrated woodland area in the northern part of the study area did not change significantly.The more scattered areas of grassland decreased slightly, while the more concentrated areas of grassland increased slightly.There was no significant change in the area of water, but the construction land area decreased significantly, mainly in the study area's western and central parts, with a slight decrease in the scattered areas in the north.Most of the construction land was converted to woodland.The area of underutilized land around the reservoir also decreased significantly, by 70%, but its baseline was small.more concentrated areas of grassland increased slightly.There was no significant change in the area of water, but the construction land area decreased significantly, mainly in the study area's western and central parts, with a slight decrease in the scattered areas in the north.Most of the construction land was converted to woodland.The area of underutilized land around the reservoir also decreased significantly, by 70%, but its baseline was small.

Comparison of Non-Point Source Pollution
We input the land use parameters before and after optimization into the SWAT model to obtain quantitative estimates of changes in total nitrogen and total phosphorus loads.

Comparison of Non-Point Source Pollution
We input the land use parameters before and after optimization into the SWAT model to obtain quantitative estimates of changes in total nitrogen and total phosphorus loads.The total nitrogen loads in the Shitoukoumen Reservoir's watershed before and after optimization amounted to 7796.83 and 7120.02tons, with an overall output reduction of 8.7%.The monthly rate and amount of change compared with those of the pre-optimization period are shown in Figure 8.The lowest total nitrogen load was 3.4 tons in January, increasing to a peak of 3166.26 tons in August, when the total nitrogen load in a single month represented 40.6% of the annual load.This was followed by a monthly decrease until the load dropped to 4.23 tons in December.
optimization, and the overall reduction the same.The reduction was mostly within 10%, but higher in May, June, and August.August was also the month with the highest reduction, but the total nitrogen load in a single month was 544.91 tons lower than before optimization.The total nitrogen load also increased slightly in February, March, September, and November.
The total phosphorus loads in the watershed before and after optimization were 798.15 and 717.39 tons, respectively, with an overall output reduction of 10.12%.The changes in monthly rates and amounts when compared with those of the pre-optimization period are shown in Figure 9.The total phosphorus load in January, February, November, and December was close to 0, while in March, May, and August it fluctuated repeatedly, rising to a highest value of 311.91 tons in August, which accounted for 43.45% of the total annual load.March, May, and August accounted for 73.6% of the total annual load, representing a more concentrated phosphorus output period.
The trend in total phosphorus load after land use optimization was similar to that before optimization.Overall it decreased, with the reduction mostly within 10%.The reduction was higher in July and the highest in August, but the highest total phosphorus load in a single month was 55.47 tons lower than before optimization.The annual total loads of nitrogen and phosphorus before and after land use optimization were input into ArcGIS to analyze their distributions.The variation in total nitrogen is shown in Figure 10.Among the 46 sub-basins, the output intensity of the total nitrogen load was the highest in sub-basins 2, 34, 35, 39, and 42, and there were clear differences in the output intensities of sub-basins 34 and 39 compared with that of the other sub-basins.These sub-basins are concentrated in the middle and upper reaches of Shuangyang River and middle reaches of the Drinking Horse and Forked Rivers.In these parts, the intensity of human activities is high, with concentrated areas of old agricultural land, relatively low forest and grassland areas, and relatively high arable land and construction land areas.Land use optimization in these areas could thus have the greatest impact on NPS pollution.The spatial distributions of total nitrogen loads before and after optimization were The trend in total nitrogen load after land use optimization was similar to that before optimization, and the overall reduction the same.The reduction was mostly within 10%, but higher in May, June, and August.August was also the month with the highest reduction, but the total nitrogen load in a single month was 544.91 tons lower than before optimization.The total nitrogen load also increased slightly in February, March, September, and November.
The total phosphorus loads in the watershed before and after optimization were 798.15 and 717.39 tons, respectively, with an overall output reduction of 10.12%.The changes in monthly rates and amounts when compared with those of the pre-optimization period are shown in Figure 9.The total phosphorus load in January, February, November, and December was close to 0, while in March, May, and August it fluctuated repeatedly, rising to a highest value of 311.91 tons in August, which accounted for 43.45% of the total annual load.March, May, and August accounted for 73.6% of the total annual load, representing a more concentrated phosphorus output period.
The trend in total phosphorus load after land use optimization was similar to that before optimization.Overall it decreased, with the reduction mostly within 10%.The reduction was higher in July and the highest in August, but the highest total phosphorus load in a single month was 55.47 tons lower than before optimization.
similar, but the intensity after optimization was slightly lower than before.Particularly in the southeast and northeast, in sub-watersheds 20, 21, and 29, the total nitrogen load was reduced to less than 5 kg/hm 2 , while the land use distribution in sub-watershed 2, located at the outlet of the watershed, did not change much, and its total nitrogen load was affected less.The annual total loads of nitrogen and phosphorus before and after land use optimization were input into ArcGIS to analyze their distributions.The variation in total nitrogen is shown in Figure 10.Among the 46 sub-basins, the output intensity of the total nitrogen load was the highest in sub-basins 2, 34, 35, 39, and 42, and there were clear differences in the output intensities of sub-basins 34 and 39 compared with that of the other sub-basins.These sub-basins are concentrated in the middle and upper reaches of Shuangyang River and middle reaches of the Drinking Horse and Forked Rivers.In these parts, the intensity of human activities is high, with concentrated areas of old agricultural land, relatively low forest and grassland areas, and relatively high arable land and construction land areas.Land use optimization in these areas could thus have the greatest impact on NPS pollution.The spatial distributions of total nitrogen loads before and after optimization were similar, but the intensity after optimization was slightly lower than before.Particularly in the southeast and northeast, in sub-watersheds 20, 21, and 29, the total nitrogen load was reduced to less than 5 kg/hm 2 , while the land use distribution in sub-watershed 2, located at the outlet of the watershed, did not change much, and its total nitrogen load was affected less.
The total phosphorus load characteristics in the Shitoukoumen Reservoir Basin are shown in Figure 11.Of the 46 sub-basins, the total phosphorus load output intensity was the greatest in sub-basins 1, 2, 4, 6, 14, 17, and 33, concentrated mainly in the middle and lower reaches of the study area.Spatial distributions of the total phosphorus loads were quite similar before and after optimization, but the intensity was significantly lower after optimization than before, especially in the non-central areas, where the total phosphorus load was less than 1.5 kg/hm 2 after land distribution adjustment.This part of the sub-basin hosts a large area of cropland and construction land that was converted into woodland.In sub-basin 33, the total phosphorus load was still higher than 2.3 kg/hm 2 after optimization, possibly at least partly due to high percentages of arable and construction land in it.Land 2024, 13, x FOR PEER REVIEW 14 of 18 The total phosphorus load characteristics in the Shitoukoumen Reservoir Basin are shown in Figure 11.Of the 46 sub-basins, the total phosphorus load output intensity was the greatest in sub-basins 1, 2, 4, 6, 14, 17, and 33, concentrated mainly in the middle and lower reaches of the study area.Spatial distributions of the total phosphorus loads were quite similar before and after optimization, but the intensity was significantly lower after optimization than before, especially in the non-central areas, where the total phosphorus load was less than 1.5 kg/hm² after land distribution adjustment.This part of the sub-basin hosts a large area of cropland and construction land that was converted into woodland.In sub-basin 33, the total phosphorus load was still higher than 2.3 kg/hm² after optimization, possibly at least partly due to high percentages of arable and construction land in it.

Discussion
The study area is a key source of drinking water but has a chronic problem of NPS pollution [57], so there are urgent needs to treat its water environment.Land use optimization has proven ability to effectively control NPS pollution [1,58].Precipitation, surface

Discussion
The study area is a key source of drinking water but has a chronic problem of NPS pollution [57], so there are urgent needs to treat its water environment.Land use optimization has proven ability to effectively control NPS pollution [1,58].Precipitation, surface runoff, soil physicochemical properties, and other variables all affect this pollution [5,54], and spatial distributions of land use affect soil quality.Thus, the effects of land use types on the surface change the hydrological characteristic, both the direction and rate of runoff, and hence the material cycling in the soil and energy distribution [35].Changes in land use also cause changes in the nature of the subsurface and water quality elements [12] and thus migration paths of pollutants and transformation rates.We found a notable reduction in the total nitrogen load within the Shitoukoumen Reservoir's watershed after land use optimization, with a 8.7% decline in the annual output and the most significant improvement in the southeastern region (where the land was transformed into denser woodland).Previously observable distinctions with the surrounding low total nitrogen load area were no longer significant, illustrating the efficacy of land use optimization for regulating the overall nitrogen output resulting from NPS pollution in the study area.Moreover, land use optimization resulted in a substantial decrease in the total phosphorus load within the watershed, with a 10.12% annual output reduction.Remarkable enhancements were observed in the southeastern area, characterized by concentrated woodland, as well as sub-basins 4 and 6 downstream, further affirming the effectiveness of land use optimization for controlling the total phosphorus output resulting from NPS pollution in the study area.The findings clearly indicate that even modest adjustments, to approximately 4% of the land area and encompassing spatial location optimization, can reduce NPS pollution intensity by approximately 10%, highlighting land use optimization's capacity to effectively mitigate such pollution.
The results also show that combining the CLUE-S land distribution and SWAT hydrological models to simulate land use optimization on a small (watershed) scale and dynamically inputting the land distribution to assess temporal and spatial effects can be used to effectively control NPS pollution and assess the effects of land optimization.This is a better solution than other methods of limiting point source pollution inputs and reducing fertilizer use in surface source pollution control [33,59].In contrast to efforts to control surface pollution by such methods, we propose here an overall scheme for optimizing regional land use from a macro-perspective, which can provide a scientific basis for the future development of land planning and policies in a watershed.While we prioritize ecological preservation in our land use optimization, this does not diminish the policy implications of our proposed plan.as we have established constraints within the optimization process, with the area allocated for land type adjustment representing less than 4% of the total land area.Moreover, the reduction in constructed land primarily focuses on areas located at a significant distance from urban centers.Another implication of our scheme for policymakers is that the spatial distribution of land use can have a significant impact on NPS pollution.Rather than spending a lot of effort to change the quantitative structure of land, perhaps adjusting the spatial distribution of land would be more effective in controlling NPS pollution.In addition, previous studies on land use optimization often only addressed the optimization, not the assessment, and specific effects of the optimization could not be verified, while our approach can provide scientific foundations for assessing optimization effects by comparing surface pollution pre-and post-land use optimization.

Conclusions
In this paper, we propose an NPS control method that combines the MODP methodology with use of the CLUE-S and SWAT models.These tools are complementary and can be used synergistically to both control NPS pollution and evaluate the specific effects of land use optimization.The proposed methodology seemed to perform well in empirical studies and could realize the goal of land use optimization to reduce NPS pollution, which can be applied in other cases.Therefore, the approach could be considered an effective tool for regional land use planning.It is also scientifically sound because we used it to achieve the objectives of assessing the effects of land optimization and controlling NPS pollution simultaneously.Moreover, the resulting optimization results can be presented in the form of land use maps, which can be used as spatially explicit tools to show locations of changes to stakeholders, thereby illustrating ways to improve efforts to meet demands for land in the future and guide land planning efforts.
Regarding the methodology, in conjunction with the MODP approach, the CLUE-S model has robust ability to optimize land use, and in combination with the SWAT model, it seems to be a robust tool for NPS pollution control and land optimization assessment.At the small watershed scale, the models involved in the methodology had satisfactory accuracy and effectiveness, reducing the intensity of NPS pollution by about 10% after adjustment of the use of about 4% of the land area and spatial location optimization in the study area.This indicates that our simple methodology has substantial advantages over the usual NPS pollution control methods.However, although it provides satisfactory spatially explicit results at a small watershed scale, challenges remain at large regional scales and transboundary regions.Therefore, one of the future goals is to collect high-resolution datasets at transboundary scales for more fine-grained land use optimization and NPS pollution control.The integrated modeling approach holds great promise.By integrating the land use model with the hydrological model, we conducted a study on controlling NPS pollution in our research area.We believe our proof-of-concept study can serve as a blueprint for integrated modeling and stimulate further integrated modeling research.A multidisciplinary approach brings forth innovative solutions, particularly concerning fields like environmental science, sociology, and economics.Environmental science offers an in-depth understanding of pollution sources, pathways, and impacts.Sociology delves into the influence of human behavior on the environment and its feedback, while economics addresses viable policies and cost-benefit analyses.The integration of these disciplines furnishes us with more comprehensive and effective methods, fostering sustainable implementation of pollution management and environmental protection measures.Future research can combine the perspectives of multiple fields to propose new solutions for NPS pollution control.
(a) Evaluating the role of the coupled model in optimizing land use and monitoring NPS pollution within the watershed; (b) Contrasting NPS pollution loads pre-and post-land use optimization, conducting an analysis on the impact of land use changes and optimization strategies on NPS pollution distribution and magnitude.(c) Assessing the effectiveness of the proposed coupled multi-model approach for controlling non-point source pollution.

Figure 1 .
Figure 1.The location of the Shitoukoumen Reservoir watershed.

Figure 1 .
Figure 1.The location of the Shitoukoumen Reservoir watershed.

Figure 4 .
Figure 4. Comparison of simulated (a) and actual (b) land use differences before and after optimization.

Figure 4 .
Figure 4. Comparison of simulated (a) and actual (b) land use differences before and after optimization.

For
the SWAT model, the set calibration period spanned 2012-2016, and the validation period consisted of the years 2017-2018.The R 2 and E NS values for runoff and ammonia nitrogen at each hydrological station are shown in Table

Figure 5 .
Figure 5. Fitted curves for the runoff values.

Figure 6 .
Figure 6.Fitted curves for the ammonia nitrogen values.

Figure 6 .
Figure 6.Fitted curves for the ammonia nitrogen values.

Figure 6 .
Figure 6.Fitted curves for the ammonia nitrogen values.

Figure 8 .
Figure 8. Variation in total nitrogen load.

Figure 8 .
Figure 8. Variation in total nitrogen load.

Figure 10 .
Figure 10.Total nitrogen load distributions (a) before and (b) after optimization (Number represents sub watershed coding).

Figure 10 .
Figure 10.Total nitrogen load distributions (a) before and (b) after optimization (Number represents sub watershed coding).

Figure 11 .
Figure 11.Total phosphorus load distributions (a) before and (b) after optimization (Number represents sub watershed coding).

Figure 11 .
Figure 11.Total phosphorus load distributions (a) before and (b) after optimization (Number represents sub watershed coding).

Table 1 .
AUC values for indicated land use types.

Table 2 .
Calibration and validation statistics for the runoff and pollutant data.

Table 2 .
Calibration and validation statistics for the runoff and pollutant data.

Table 3 .
Land use cover optimization.

Table 3 .
Land use cover optimization.