Assessment of Non-point Source Total Phosphorus Pollution from Different Land Use and Soil Types in a Mid-high Latitude Region of China

The transport characteristics of phosphorus in soil and the assessment of its environmental risk have become hot topics in the environmental and agricultural fields. The Sanjiang Plain is an important grain production base in China, and it is characterised by serious land use change caused by large-scale agricultural exploitation. Agricultural inputs and tillage management have destroyed the soil nutrient balance formed over long-term conditions. There are few studies on non-point source phosphorus pollution in the Sanjiang Plain, which is the largest swampy low plain in a mid-high-latitude region in China. Most studies have focused on the water quality of rivers in marsh areas, or the export mechanism of phosphorus from specific land uses. They were conducted using experimental methods or empirical models, and need further development towards mechanism models and the macro-scale. The question is how to find a way to couple processes in phosphorus cycling and a distributed hydrological model considering local hydrological features. In this study, we report an attempt to use a distributed phosphorus transport model to analyse non-point source total phosphorus pollution from different land uses and soil types on the Sanjiang Plain. The total phosphorus concentration generally shows an annually increasing trend in the study area. The total phosphorus load intensity is heterogeneous in different land use types and different soil types. The average total phosphorus load intensity of different land use types can be ranked in descending order from paddy field, dry land, wetlands, grassland, and forestland. The average total phosphorus load intensity of different soil types can be ranked in descending order: paddy soil, bog soil, planosol, meadow soil, black soil, and dark brown earth. The dry land and paddy fields account for the majority of total phosphorus load in the study area. This is mainly caused by extensive use of phosphate fertilizer on the cultivated land. This has important implications for future agricultural management and non-point source control in this agricultural area of the mid-high latitude region.


Introduction
Phosphorus is one of a large number of elements that are necessary for plant growth and development.It is also the main limiting factor in crop yield [1], and is of significance in global food security.In many parts of the world, in order to meet the growing populations' demand for food, it is necessary to add mined rock phosphate or artificial phosphorus fertilizer to the cultivated land to maintain food production [2].A large amount of phosphorus accumulation in the soil enters into water bodies via soil erosion, causing the phosphorus content to become elevated and the water bodies to undergo eutrophication.When the total phosphorus content of water is more than 0.02 mg/L, it is likely to lead to algal blooms [3], and therefore non-point source (NPS) pollution of phosphorus caused by soil erosion has generated wide public concern.
Phosphorus is lost from soils through several processes including leakage, surface runoff, and soil erosion.The loss of phosphorus is closely related to the physical structure of the soil, field capacity, rainfall, vegetation coverage, fertilizer application, and irrigation methods.Phosphorus load from surface runoff is three or four times that from the subsurface flow.Thus, the main part of the phosphorus output is suspended total phosphorus, while the percentage of dissolved total phosphorus and proportion of orthophosphate are very low [4].Migration of dissolved phosphorus through the soil occurs in two ways: one is a slow flow through normal (matrix) soil pores, and the other is rapid (preferential) flow involving macropores.The two types together contribute to the seepage erosion of soil phosphorus [5].
The Sanjiang Plain is one of the largest grain production bases in China and has experienced rapid land cover change, making it a critical area in assessing the environmental impact of LUCC and proposing strategies for minimizing its impact on non-point source (NPS) pollution.Since the 1950s, the Sanjiang Plain has been subjected to massive agricultural development.A large number of marsh wetlands were drained and transformed into fertile land for cultivation, and deforestation and destruction of local grasslands have taken place.Consequently, significant changes have taken place in the mode of regional land use.The application of pesticides and chemical fertilizers, and changes in farming methods, have destroyed the long-standing soil nutrient equilibrium, and accelerated the migration rate of soil nutrients across the water/soil interface.These nutrients are mainly nitrogen and phosphorus that reach water bodies as pollutants, and have resulted in serious environmental pollution.
In recent years, there are many mechanism models that can be used to simulate and describe the migration and conversion of nutrients, such as the PSYCHIC model, the CREAMS model, the ANSWERS model, the AnnAGNPS model, the EPIC model, and the SWAT model [6][7][8][9][10][11][12].Among the above models, the SWAT model is the most widely used, and has been applied in North America, Europe, Africa, South Asia and other regions, as well as in China, and other countries.Several case studies analysing the impact of land cover changes on pollution via NPS nutrients have been conducted using SWAT model [13][14][15].Although the SWAT model is suitable for long-term continuous non-point phosphorus simulation in agricultural areas [16], the SWAT model, which simplifies the model runs by dividing the hydrological response unit, belongs to a semi-distributed model.In addition, the SCS curve method, which is used in the hydrology module of SWAT, is an empirical model, and it has difficulties simulating the effects of best management practices (BMPs) on runoff because of its reliance on curve number methodology [17].Models such as the ANSWERS model, which can only simulate at the field-scale or small watershed scale, are not suitable for large-scale phosphorus cycle simulations.Although remote sensing technology cannot directly detect water quality, it can obtain the spatial information of NPS pollution model parameters.Spatial data can clearly meet the needs of digital simulation on NPS pollution.Thus, coupling multi-source remote sensing information, data assimilation, and a distributed NPS pollution model is the main trend in NPS pollution model development.
The Ecohydrological Assessment Tool (EcoHAT) [18][19][20][21][22][23][24][25] is a type of ecological hydrology simulation system, which can be used to comprehensively simulate regional ecological and hydrological processes based on physical and chemical mechanisms.With EcoHAT, hydrological cycles, nutrient cycles and plant growth processes are integrated, while nutrient cycling and growth of vegetation are considered in relation to the hydrological cycle.EcoHAT, coupled with a remote sensing model to invert the land surface parameters, provides spatial data for ecological and hydrological process simulation.Driven by remote sensing data, EcoHAT makes full use of remote sensing technology in solving large-scale soil phosphorus, NPS pollution problems with time-series simulations.Meanwhile, EcoHAT integrates GIS spatial analysis tools, which can effectively analyse the simulation results [19].
The NPSs of phosphorus pollution are important, such as leakage (loss of phosphorus), which may cause adverse effects in surface and ground water.The Sanjiang Plain is an important agricultural region with various types of land use and soil.Analysis of the amount of phosphorus load under different land use and soil types is helpful for protecting the water bodies, and for controlling NPS pollution in this agricultural area of the mid-high latitude region [26].However, there are few studies on non-point pollution by phosphorus in the Sanjiang Plain, and most of these focused on water quality of the river in marsh areas, or on the export mechanisms of phosphorus from specific land uses or specific soil types [27].The studies were done using experimental methods or empirical models, and the research requires further advances towards mechanism model development at the macro-scale.Questions remain about how to find a way to couple the processes in phosphorus cycling and a distributed hydrological model, with consideration for local hydrological features.Furthermore, studies of the non-point-source pollution by phosphorus in the Sanjiang Plain require solutions that can only be provided by access to the spatial data needed to drive the pertinent models, and to expand their scale.
To ensure good soil quality during agricultural development, the trend of changes in soil phosphorus due to land use change must be assessed.Understanding the characteristics of changes in soil phosphorus content is important for improving agricultural management practices and making land use practice more sustainable.The aims of this study are as follows: (1) to identify the temporal-spatial pattern of land use change occurring from 2000 to 2010 (11a) in the Sanjiang Plain; (2) to explore the transformation of soil phosphorus from six main soil types and five types of land use; and (3) to provide suggestions for future farmland management and establishment of sustainable land use systems while preserving soil phosphorus contents.
The climate of the Sanjiang Plain, and hence the study area, is "continental monsoon climate", with four distinct seasons, and with heat and rainfall occurring together.The average annual precipitation is 406 mm, and most rainfall occurs from July to August, a period that accounts for 70% of annul precipitation.The yearly average temperature ranges from 2 to 8 • C. The highest temperatures are registered from June to August, and the average annual summer temperature is 21.6 to 27.2 • C.Under these conditions, vegetation grows from April through September.Average monthly wind speeds vary from 3.1 to 4.1 m•s −1 .The research focused only on the 5 km riparian buffer zone around the Sanjiang Plain.The main land use types in this area are dry land, paddy field, forestland, pasture, and wetland (Figure 2).Due to widespread agricultural and domestic NPS pollution, the Sanjiang Plain is a suitable area for studying NPS pollution problems.Water 2016, 8, 505 5 of 17

Data
The data (Table 1) used in the model came from multi-source remote sensing data (e.g., Landsat TM or MODIS), in situ experimental data, meteorological data, and survey data.Meteorological data came from public, free, data platforms (e.g.,China meteorological data net [28]).Multi-source remote sensing data were also from free public platforms (e.g., USGS: United States Geological Survey [29]).Other data were from in situ experiments and statistics.All original input data were processed by programs written in IDL (Interactive Data Language).Some data were resolved using ENVI4.8(Harris Corporation, Melbourne, FL, USA) or ArcGIS10.1 (Environmental Systems Research Institute, Redlands, CA, USA).All the input data were converted to images, which were then projected into a uniform coordinate system to render the disparate sets of spatial data consistent.The spatial resolution of the input data was also a key factor, because only if the input data have the same spatial resolution can the P simulation model run smoothly.However, the different input data have very different spatial resolutions.For instance, the spatial resolution of the Landsat TM (Thematic Mapper) data we used is about 30 m, while the spatial resolution of MODIS (Moderate-Resolution Imaging Spectral Radiometer) data is 1000 m.We tried to use GLDAS (Global Land Data Assimilation System) data to improve our spatial-temporal (three-hourly) resolution of meteorological data, but the GLDAS data had a coarse spatial resolution of 25 km.To obtain uniform 1000-m spatial resolution, we needed to input data with different spatial resolutions into the model, so we resampled, and used some downscaling methods [30].Validation of the downscaling result was satisfactory.After a three-year land use overlay analysis, the area and distribution of permanent dry land (D-D), permanent paddy field (P-P), permanent forestland (F-F), and permanent wetland (W-W) were identified for the period from 2000 to 2010.In addition, the spatial information on four types of land use conversion in the same period were analysed: dry land reclaimed from forestland (F-D), paddy field reclaimed from dry land (D-P), dry land reclaimed from wetland (W-D), and paddy field reclaimed from wetland (W-P) (Table 2).
According to the USDA soil classification system [31], there are various types of soil in the Sanjiang Plain.These include black, dark brown, swamp, meadow, white pulp, paddy, peat, alluvial, grey brown, rocky, volcanic ash and brown coniferous forest soils.The six main soil types (dark brown, swamp, white pulp, meadow, black and paddy soil) account for 98.7% of the soil area over the entire Sanjiang Plain.The spatial distribution of soil types is shown in Figure 3. Dark brown soils are mainly distributed in hilly areas.White pulp soils are mainly distributed in grade I valleys at all levels on the terrace.The black soil is distributed at the top of the western piedmont platform.Meadow soils are mainly distributed in the floodplain of the main rivers.The paddy soil is found in river flood plains and on grade I terraces.Swamp soils are mainly distributed in the central low-lying areas, and are perennially wet.The lower part of the upper peat soil is mineralized to form grey soils, but if the peat layer thickness is greater than 50 cm, it is defined as a peat soil.The Sanjiang Plain vegetation types mainly include forest, swamp meadow plants and marsh plants, which belong to the Changbai flora.The low hilly land around the forest is mainly found on the plain, and its vegetation types are divided into coniferous forest, mixed needleleaf-broadleaf, and deciduous broadleaf forest.Natural wetland vegetation is mainly distributed in Heilongjiang.Along the Songhua and Ussuri Rivers and their tributaries, there are flood lands, lakes and lake depressions with mainly aquatic grass, wet meadows and shrubs.
Water 2016, 8, 505 6 of 17 Meadow soils are mainly distributed in the floodplain of the main rivers.The paddy soil is found in river flood plains and on grade I terraces.Swamp soils are mainly distributed in the central low-lying areas, and are perennially wet.The lower part of the upper peat soil is mineralized to form grey soils, but if the peat layer thickness is greater than 50 cm, it is defined as a peat soil.The Sanjiang Plain vegetation types mainly include forest, swamp meadow plants and marsh plants, which belong to the Changbai flora.The low hilly land around the forest is mainly found on the plain, and its vegetation types are divided into coniferous forest, mixed needleleaf-broadleaf, and deciduous broadleaf forest.Natural wetland vegetation is mainly distributed in Heilongjiang.Along the Songhua and Ussuri Rivers and their tributaries, there are flood lands, lakes and lake depressions with mainly aquatic grass, wet meadows and shrubs.

The Phosphorus Model in EcoHAT-P
Total phosphorus (TP) is a measure of all the forms of phosphorus, such as active organic P, dissolved inorganic P, and particulate inorganic P. The increase in phosphorus in the soil is mainly from fertilization and plant litter decomposition.The ways in which phosphorus output from the soil occurs include absorption by plants, surface runoff through rainfall or artificial drainage, soil erosion and leaching.The phosphorus cycle is affected by many reactions and mechanisms.The simulation methods that are usually used, summarize the forms of organic phosphorus and inorganic phosphorus present, generally include the processes of decomposition, fixing/mineralization, adsorption/desorption and plant absorption.Some models increase the leaching and anthropogenic factors.
The phosphorus circular migration module is built up from processes in the nutrient-element cycles, vegetation absorption of nutrients, net primary productivity (NPP), the distribution of productive forces and the process of vegetation littering.We used Rs-DTVGM [32] hydrological models, and MUSLE [33,34] equations coupling with the phosphorus cycle migration module, to set up the distributed phosphorus migration (EcoHAT-P) model (Figure 4).
After we used multi-source remote sensing data, and hydrologic and meteorological observation data to establish an ecohydrological spatial database, a distributed phosphorus migration model driven by remote sensing data was realized.Finally, the content of dissolved phosphorus and adsorbed phosphorus in surface runoff was obtained.The sum of these forms of phosphorus is TP.This provided data for further study on the spatiotemporal changes in phosphorus.The model structure is shown in Figure 4. Total phosphorus (TP) is a measure of all the forms of phosphorus, such as active organic P, dissolved inorganic P, and particulate inorganic P. The increase in phosphorus in the soil is mainly from fertilization and plant litter decomposition.The ways in which phosphorus output from the soil occurs include absorption by plants, surface runoff through rainfall or artificial drainage, soil erosion and leaching.The phosphorus cycle is affected by many reactions and mechanisms.The simulation methods that are usually used, summarize the forms of organic phosphorus and inorganic phosphorus present, generally include the processes of decomposition, fixing/mineralization, adsorption/desorption and plant absorption.Some models increase the leaching and anthropogenic factors.
The phosphorus circular migration module is built up from processes in the nutrient-element cycles, vegetation absorption of nutrients, net primary productivity (NPP), the distribution of productive forces and the process of vegetation littering.We used Rs-DTVGM [32] hydrological models, and MUSLE [33,34] equations coupling with the phosphorus cycle migration module, to set up the distributed phosphorus migration (EcoHAT-P) model (Figure 4).
After we used multi-source remote sensing data, and hydrologic and meteorological observation data to establish an ecohydrological spatial database, a distributed phosphorus migration model driven by remote sensing data was realized.Finally, the content of dissolved phosphorus and adsorbed phosphorus in surface runoff was obtained.The sum of these forms of phosphorus is TP.This provided data for further study on the spatiotemporal changes in phosphorus.The model structure is shown in Figure 4.

Phosphorus Modules in EcoHAT
Most of EcoHAT-P is composed of two modules, the plant phosphorus cycle and the soil phosphorus cycle.Phosphorus biological cycles in vegetation are mainly associated with growth of vegetation.The vegetation growth process can be divided into four sub-models: vegetation net primary productivity (NPP), productivity distribution, nutrient absorption and vegetation litter.
The circular migration model of phosphorus in soil mainly refers to the SWAT model [35], and according to the need of the application, the SWAT model was simplified such that decomposition and mineralization considered only the leaf litter and upper soil fresh litter.Mobility and transfer of

Phosphorus Modules in EcoHAT
Most of EcoHAT-P is composed of two modules, the plant phosphorus cycle and the soil phosphorus cycle.Phosphorus biological cycles in vegetation are mainly associated with growth of vegetation.The vegetation growth process can be divided into four sub-models: vegetation net primary productivity (NPP), productivity distribution, nutrient absorption and vegetation litter.
The circular migration model of phosphorus in soil mainly refers to the SWAT model [35], and according to the need of the application, the SWAT model was simplified such that decomposition Water 2016, 8, 505 8 of 17 and mineralization considered only the leaf litter and upper soil fresh litter.Mobility and transfer of phosphorus in soil mainly refers to the phosphorus model with simplified soil and vegetation proposed by Jones [36], which elaborated details about the mobility and migration characteristics of phosphorus at the interface between soil and water.The main functions and parameters in the P model are listed in the Table 3.
Table 3. Functions and parameters in the P model.

No.
Model Name Equation Reference CASA model [37] 2 Production Distribution [38,39] 4 Vegetation Litter Fertilization The dissolved soil P from fertilization (kg P•ha −1 ); fert minP : The percentage of mineral P in fertilizer (%); orgP frsh,fert : The fresh organic P in soil come from fertilization (kg P•ha −1 ); fert orgP : The percentage of organic P in fertilization (%); orgP hum,fert : The humus organic soil P come from fertilization (kg P•ha −1 ); fert: The amount of fertilization (kg P•ha −1 ); P mina : The phosphorus mineralized from the humus active organic P pool (kg P•ha −1 ); orgP act : The mount of phosphorus in the active organic pool (kg P•ha −1 ); β min : The rate coefficient for mineralization of the humus active organic nutrients; γ tmp : The nutrient cycling temperature factor; γ sw : The nutrient cycling water factor; P sol/act,ly : Amount of phosphorus transferred between the active and stable mineral pools (kg P•ha −1 ); P slution,ly : The amount of phosphorus in solution (kg P•ha −1 ); minP act,ly : The amount of phosphorus in the active mineral pool (kg P•ha −1 ); pai: The availability index.

Validation of the EcoHAT-P Model
Through the random method, we got 41 sampling points, obtained their geography coordinates, and overlaid the coordinates of sampling points on the simulation images to get the simulated results.In Table 4, the simulated mean value of total phosphorus (TP) is 0.878 g/kg, and the observed mean value is 0.881 g/kg.Other statistical standards, which include the maximum values, standard deviations, and standard error of mean, were also used to show the distribution characteristics between the two data groups.We obtained observed and simulated values and used the data to represent the conditions of this area.In Figure 5, the test results show that R 2 of the observed and the simulated data is 0.661, demonstrating that the P model is reliable.Through the random method, we got 41 sampling points, obtained their geography coordinates, and overlaid the coordinates of sampling points on the simulation images to get the simulated results.In Table 4, the simulated mean value of total phosphorus (TP) is 0.878 g/kg, and the observed mean value is 0.881 g/kg.Other statistical standards, which include the maximum values, standard deviations, and standard error of mean, were also used to show the distribution characteristics between the two data groups.We obtained observed and simulated values and used the data to represent the conditions of this area.In Figure 5, the test results show that R 2 of the observed and the simulated data is 0.661, demonstrating that the P model is reliable.

Total Phosphorus during the Period
In detecting abrupt changes with a non-parameter Mann-Kendall test method, the results showed that the yearly total phosphorus load abruptly changed in 2005.The study period was divided into two intervals: 2000-2005, and 2006-2010.Figure 6 shows that the total phosphorus load in the Sanjiang Plain has a trend of overall increase.The P-load presents a relatively low growth rate during 2000-2005 (0.0036 g/(kg•a)), then increased during 2006-2010 (0.0122 g/(kg•a)).During 2006-2010, the growth rate of the total phosphorus load presented a trend of rapid growth for the Sanjiang Plain.This is because much wetland (~1800 km 2 ), pasture (~1500 km 2 ), and woodland (~1200 km 2 ) were converted to paddy field (~3500 km 2 ) and dry land (~1000 km 2 ), resulting in a rise in the amount of fertilizer application, and thus total phosphorus.The p value was 0.1569, which is larger than 0.01.The overall correlation between NPS mean total phosphorus and precipitation was not significant.

Total Phosphorus during the Period
In detecting abrupt changes with a non-parameter Mann-Kendall test method, the results showed that the yearly total phosphorus load abruptly changed in 2005.The study period was divided into two intervals: 2000-2005, and 2006-2010.Figure 6 shows that the total phosphorus load in the Sanjiang Plain has a trend of overall increase.The P-load presents a relatively low growth rate during 2000-2005 (0.0036 g/(kg•a)), then increased during 2006-2010 (0.0122 g/(kg•a)).During 2006-2010, the growth rate of the total phosphorus load presented a trend of rapid growth for the Sanjiang Plain.This is because much wetland (~1800 km 2 ), pasture (~1500 km 2 ), and woodland (~1200 km 2 ) were converted to paddy field (~3500 km 2 ) and dry land (~1000 km 2 ), resulting in a rise in the amount of fertilizer application, and thus total phosphorus.The p value was 0.1569, which is larger than 0.01.The overall correlation between NPS mean total phosphorus and precipitation was not significant.

Total Phosphorus Loads from Different Land Uses
A statistical test (one-way ANOVA) was applied to evaluate the significance of differences in TP load intensity among different land use types.The P value was 0.0012, showing that there was significant difference in TP load intensity among different land use types (p < 0.05).The total phosphorus load intensity is heterogeneous in different land use types.The average value of the

Total Phosphorus Loads from Different Land Uses
A statistical test (one-way ANOVA) was applied to evaluate the significance of differences in TP load intensity among different land use types.The P value was 0.0012, showing that there was significant difference in TP load intensity among different land use types (p < 0.05).The total phosphorus load intensity is heterogeneous in different land use types.The average value of the total phosphorus load intensity in different land use types can be ranked in descending order as paddy field, dry land, wetlands, grassland, and forestland (Figure 7a).

Total Phosphorus Loads from Different Land Uses
A statistical test (one-way ANOVA) was applied to evaluate the significance of differences in TP load intensity among different land use types.The P value was 0.0012, showing that there was significant difference in TP load intensity among different land use types (p < 0.05).The total phosphorus load intensity is heterogeneous in different land use types.The average value of the total phosphorus load intensity in different land use types can be ranked in descending order as paddy field, dry land, wetlands, grassland, and forestland (Figure 7a).
The total phosphorus load for dry land and paddy field accounted for 64.5%, suggesting that the dominant human agricultural activities have played a major role in the occurrence of phosphorus pollution in the study area.The total phosphorus load of forestland accounted for 25% of that in the study area, but the high forestland total phosphorus load is caused by the decomposition and mineralization of forest litter (Figure 7b).The total phosphorus load for dry land and paddy field accounted for 64.5%, suggesting that the dominant human agricultural activities have played a major role in the occurrence of phosphorus pollution in the study area.The total phosphorus load of forestland accounted for 25% of that in the study area, but the high forestland total phosphorus load is caused by the decomposition and mineralization of forest litter (Figure 7b).
The total phosphorus load intensity of paddy field increased rapidly, ranging from 9.72821 kg/ha in 2000 to 11.58 kg/ha in 2010, a 0.1683 kg/ha growth per year.The total phosphorus load intensity of dry land increased rapidly, ranging from 8.72307 kg/ha in 2000 to 10.606 kg/ha in 2010, a 0.1712 kg/ha growth per year.The total phosphorus load intensity of wetlands increased rapidly, ranging from 8.9662 kg/ha in 2000 to 9.678 kg/ha in 2010, a 0.0647 kg/ha growth per year.The total phosphorus load intensity of grassland decreased rapidly, ranging from 8.32733 kg/ha in 2000 to 7.7688 kg/ha in 2010, a decline of 0.0508 kg/ha per year.The total phosphorus load intensity of forestland decreased rapidly, ranging from 7.98461 kg/ha in 2000 to 7.3731 kg/ha in 2010, a decline of 0.0556 kg/ha per year.
Based on the data, the total phosphorus load intensity of paddy field and dry land increased year on year, while the other three classes (forestland, pasture and wetland) all showed a declining trend.The increasing trend of the total phosphorus load intensity in the paddy field and dry land is mainly due to the rise in the usage of phosphate fertilizer (Figure 8a).
The total phosphorus load for paddy field and dry land are increasing year on year, while forestland and pasture decline steadily.It appears that there is a smooth trend in wetland after the year 2004.This is because the policies of "Reverting Farmland to Wetland" and "Reverting Farmland to Pasture" on the Sanjiang Plain have been very well implemented.They have reduced the rate of the reduction in the area of wetland and pasture, while the loss of forestland area was still in substantial decline.By contrast, the area of paddy field and dry land showed steady growth (Figure 8b).
The total phosphorus load for paddy field and dry land are increasing year on year, while forestland and pasture decline steadily.It appears that there is a smooth trend in wetland after the year 2004.This is because the policies of "Reverting Farmland to Wetland" and "Reverting Farmland to Pasture" on the Sanjiang Plain have been very well implemented.They have reduced the rate of the reduction in the area of wetland and pasture, while the loss of forestland area was still in substantial decline.By contrast, the area of paddy field and dry land showed steady growth (Figure 8b).

Total Phosphorus Loads from Different Soil Types
The total phosphorus loads of four soil types (meadow, dark brown, planosol and bog) accounted for 90.84% of the total phosphorus load of all soil types.The largest share of total phosphorus in paddy fields is meadow soil, whereas paddy soil is the main soil component in paddy fields.Meadow soil is the main soil type of dry land (Figure 9a,b).
The one-way ANOVA test was applied to evaluate the significance of differences in TP load intensity among different soil types.The p value was 0.00005, showing that there was significant difference in TP load intensity among different soil types (p < 0.05).The total phosphorus load intensity is heterogeneous among different soil types.The average value of the total phosphorus load intensity in the main soil types can be ranked in descending order as paddy soil, bog soil, planosol, meadow soil, black soil, and dark brown earth between 2001 and 2010 (Figure 9c,d).

Total Phosphorus Loads from Different Soil Types
The total phosphorus loads of four soil types (meadow, dark brown, planosol and bog) accounted for 90.84% of the total phosphorus load of all soil types.The largest share of total phosphorus in paddy fields is meadow soil, whereas paddy soil is the main soil component in paddy fields.Meadow soil is the main soil type of dry land (Figure 9a,b).The one-way ANOVA test was applied to evaluate the significance of differences in TP load intensity among different soil types.The p value was 0.00005, showing that there was significant difference in TP load intensity among different soil types (p < 0.05).The total phosphorus load intensity is heterogeneous among different soil types.The average value of the total phosphorus load intensity in the main soil types can be ranked in descending order as paddy soil, bog soil, planosol, meadow soil, black soil, and dark brown earth between 2001 and 2010 (Figure 9c,d

Effect of Land Use Changes on Total Phosphorus
The transformation of soil types is relatively slow in nature, while the transformation of land use types is more frequent compared with soil types.Since Sanjiang Plain is an important agricultural area, the changes of land use are driven by economic interests (e.g., the economic benefit of rice is high, so much dry land are transformed into paddy field) or by policy (e.g., cultivated fields are transformed into the wetland or pasture).To ensure high yields, fertilizers have generally been applied at very high rates by farmers in this area [41].We believe that the increased application of phosphate fertilizer is one of the most important reasons for the significant change in P-loading by the end of the year 2005.Even the total phosphorus load strength is relatively stable in different land types, the area changes of some land use types will lead to total phosphorus load changes of different land use types, resulting in change in the total phosphorus load in the study area.In this study, the TP load intensities were 0.778 g/kg (Dryland), 0.814 g/kg (Paddy field), 0.620 g/kg (Forest), and 0.759 g/kg (Wetland).The TP load intensities in Ouyang's study are 0.878 ± 0.101 g/kg (Dryland), 0.878 ± 0.101 g/kg (Paddy field), 0.841 ± 0.369 g/kg (Forest), and 1.357 ± 0.780 g/kg (Wetland).The content of total phosphorus load Water 2016, 8, 505 13 of 17 intensity is similar to the results of the study by Ouyang et al. [42].However, it is much more efficient to use the mechanism model EcoHAT-P to estimate the spatial and temporal variations of NPS phosphorus pollution for a large region than the previous regional experimental sampling methods [41,43].Thus, the model is expected to greatly reduce the management cost in large agriculture regions.

Response of NPS Phosphorus Load with Temperature
The study area is located at high latitude, as is apparent from Figure 11a.The trend of the average temperature changed in the study area was very similar to the trend of average soil temperature during 2005-2010 (a unimodal distribution).The lowest and highest temperatures were in January and July, respectively, and from November to March of the next year, the air temperature and soil temperatures were below 0 • C. The trend of total phosphorus load appears to be a skewed distribution, with the most severe change mainly in May, June and July.This arises because these months coincide with the cropping season, when the phosphorus is applied (Figure 11b).It can be seen from Figure 12a that the main land use types that exhibited month to month fluctuations were dry land and paddy field.It can be seen from Figure 12b that the main soil types that exhibited month to month fluctuations were meadow soil, planosol and bog soils.
the end of the year 2005.Even the total phosphorus load strength is relatively stable in different land types, the area changes of some land use types will lead to total phosphorus load changes of different land use types, resulting in change in the total phosphorus load in the study area.In this study, the TP load intensities were 0.778 g/kg (Dryland), 0.814 g/kg (Paddy field), 0.620 g/kg (Forest), and 0.759 g/kg (Wetland).The TP load intensities in Ouyang's study are 0.878 ± 0.101 g/kg (Dryland), 0.878 ± 0.101 g/kg (Paddy field), 0.841 ± 0.369 g/kg (Forest), and 1.357 ± 0.780 g/kg (Wetland).The content of total phosphorus load intensity is similar to the results of the study by Ouyang et al. [42].However, it is much more efficient to use the mechanism model EcoHAT-P to estimate the spatial and temporal variations of NPS phosphorus pollution for a large region than the previous regional experimental sampling methods [41,43].Thus, the model is expected to greatly reduce the management cost in large agriculture regions.

Response of NPS Phosphorus Load with Temperature
The study area is located at high latitude, as is apparent from Figure 11a.The trend of the average temperature changed in the study area was very similar to the trend of average soil temperature during 2005-2010 (a unimodal distribution).The lowest and highest temperatures were in January and July, respectively, and from November to March of the next year, the air temperature and soil temperatures were below 0 °C.The trend of total phosphorus load appears to be a skewed distribution, with the most severe change mainly in May, June and July.This arises because these months coincide with the cropping season, when the phosphorus is applied (Figure 11b).It can be seen from Figure 12a that the main land use types that exhibited month to month fluctuations were dry land and paddy field.It can be seen from Figure 12b that the main soil types that exhibited month to month fluctuations were meadow soil, planosol and bog soils.The study area is located at higher latitude in a freeze-thaw region.This provides a possible condition by which to simulate the effect of NPS phosphorus pollution on change in freeze-thaw effects in the watershed.
The variation of NPS phosphorus load with changes in air and ground temperatures is a The study area is located at higher latitude in a freeze-thaw region.This provides a possible condition by which to simulate the effect of NPS phosphorus pollution on change in freeze-thaw effects in the watershed.
The variation of NPS phosphorus load with changes in air and ground temperatures is a complex process.In freeze-thaw cycles, biological activity is also closely related with temperature.Normally, with increasing temperature, the rate of biological processes would also increase.The temperature also controls water movement in the soil, and the movement of soil water strongly influences the migration and transformation of NPS phosphorus.In this study, there was a minimum phosphorus value in the April simulation, and in April thawing of the soil occurs when the soil temperature is above the freezing point.The peak value was in May (Figure 11b), which was due to peak fertilization for agriculture.The plot of the NPS phosphorus load had a skewed peak pattern, which is further validation that on-site sampling data are needed for further study.

Uncertainty for Agricultural NPS Phosphorus Pollution Control
The accuracy of the spatial data is an important factor that impacts on the accuracy of the simulation of phosphorus migration on the Sanjiang Plain.In this study, the main spatial data were derived from multiple-sourced remote sensing products provided by a common platform or model-based remote sensing inversion.Each product itself contains some errors, and remote sensing products from different data sources also have different resolutions in time and space.We used a simple resampling method to transform their scales, which inevitably leads to some errors, increasing the uncertainty of the model.The cumulative errors produce impacts on the final simulation results.
The spraying of chemical fertilizers and pesticides is an important source of agricultural NPS phosphorus pollution.Only phosphorus fertilizer application is considered in this process, and the study does not reflect the impact of spraying pesticides on the NPS phosphorus pollution load.In addition, it does not consider the impacts from farmland ridges and ditches for phosphorus retention.This may also cause uncertainty in the NPS results regarding the phosphorus pollution load.

Suggestions for Sustainable Agricultural Land Management in the Future
As can be seen from Figure 5, there was a rapid growth trend in the NPS phosphorus pollution load in Sanjiang Plain from 2000 to 2010.This is explained by the increasing amount of phosphate fertilizer used in the cultivated fields.In order to ensure that the NPS phosphorus pollution can be effectively controlled, farmers need to undertake management measures according to the total phosphorus load intensity in different soil types or land use types in the study area (Figure 8a,b).Such actions will reduce the waste of phosphate fertilizer, and avoid the NPS phosphorus pollution caused by them.In addition, farmers need to consider controlling the use of pesticides.In this study, the relationship between the NPS phosphorus load and the temperature is relatively coarse.We can conduct a further study to determine the relationship between the two by increasing the numbers of soil temperature monitoring sites.Then, we will be able to provide data and theoretical support for agricultural NPS pollution control in the area, which is located in the mid-high latitude areas, such as Sanjiang Plain.

Conclusions
We applied a distributed NPS model driven by remote sensing data to simulate the phosphorus NPS pollution load on the Sanjiang Plain from 2000 to 2010.The NPS phosphorus pollution load was analysed according to land use types and soil types.We found that the NPS phosphorus pollution load on the Sanjiang Plain shows a trend of increasing year on year, and that the NPS phosphorus pollution load is heterogeneous for different land use types and soil types.The dry land and paddy field account for the majority of total phosphorus load in the study area.This is mainly caused by extensive use of phosphate fertilizers on cultivated land.It can be seen that intensive human agricultural activities play a major role in the phosphorus pollution of the study area.If we do not effectively control the agricultural use of chemical fertilizers, there is a significant risk that the NPS phosphorus pollution load will continue to rise on the Sanjiang Plain.

Figure 1 .
Figure 1.Geographical location and DEM (Digital Elevation Model) of the Sanjiang Plain in Northeast China.

Figure 2 .
Figure 2. Land use distributions and conversion from 2000 to 2010.

Figure 1 . 17 Figure 1 .
Figure 1.Geographical location and DEM (Digital Elevation Model) of the Sanjiang Plain in Northeast China.

Figure 2 .
Figure 2. Land use distributions and conversion from 2000 to 2010.Figure 2. Land use distributions and conversion from 2000 to 2010.

Figure 2 .
Figure 2. Land use distributions and conversion from 2000 to 2010.Figure 2. Land use distributions and conversion from 2000 to 2010.

Figure 3 .
Figure 3. Distribution of soil types in Sanjiang Plain.

Figure 3 .
Figure 3. Distribution of soil types in Sanjiang Plain.

1 .
The Phosphorus Model in EcoHAT-P

Figure 4 .
Figure 4. Conceptual framework of the soil P model.

Figure 4 .
Figure 4. Conceptual framework of the soil P model.

Figure 5 .
Figure 5. Observed and simulated results for the random sample points.

Figure 5 .
Figure 5. Observed and simulated results for the random sample points.

Figure 6 .
Figure 6.Annual watershed precipitation and NPS mean total phosphorus (TP) yield from 2000 to 2010.

Figure 6 .
Figure 6.Annual watershed precipitation and NPS mean total phosphorus (TP) yield from 2000 to 2010.

Figure 7 .
Figure 7. Load intensity (a); and load (b) of total phosphorus (TP) from each land use between 2001 and 2010.

Figure 7 .
Figure 7. Load intensity (a); and load (b) of total phosphorus (TP) from each land use between 2001 and 2010.

Figure 8 .
Figure 8. Load intensity (a); and load (b) of total phosphorus (TP) for each land use between 2000 and 2010.

Figure 8 .
Figure 8. Load intensity (a); and load (b) of total phosphorus (TP) for each land use between 2000 and 2010.

Water 2016, 8 , 505 12 of 17 Figure 9 .
Figure 9. Proportion of soil types in different land use types (a); proportion of land use types in different soil types (b); load (c); and load intensity (d) of total phosphorus (TP) for each soil type between 2001 and 2010.

Figure 9 .
Figure 9. Proportion of soil types in different land use types (a); proportion of land use types in different soil types (b); load (c); and load intensity (d) of total phosphorus (TP) for each soil type between 2001 and 2010.
).The total phosphorus load intensity of paddy soil increased rapidly, ranging from 11.2404 kg/ha in 2001 to 13.0769 kg/ha in 2010, a 0.18365 kg/ha growth per year.The total phosphorus load intensity of bog soil increased rapidly, ranging from 10.5554 kg/ha in 2001 to 12.0433 kg/ha in 2010, a 0.14879 kg/ha growth per year.The total phosphorus load intensity of planosol increased rapidly, ranging from 8.67532 kg/ha in 2001 to 10.677 kg/ha in 2010, a 0.200168 kg/ha growth per year.The total phosphorus load intensity of meadow soil increased rapidly, ranging from 7.93381 kg/ha in 2001 to 9.47273 kg/ha in 2010, a 0.153892 kg/ha growth per year.The total phosphorus load intensity of black soil increased rapidly, ranging from 7.01221 kg/ha in 2001 to 8.77915 kg/ha in 2010, a 0.176694 kg/ha growth per year.The total phosphorus load intensity of dark brown soil decreased rapidly, ranging from 7.65015 kg/ha in 2001 to 7.28174 kg/ha in 2010, a decline of 0.036841 kg/ha per year.Based on the data, the total phosphorus load and total phosphorus load intensity are rising in most of the major soil types, with only dark brown soil showing a declining trend (Figure10).

Figure 9 .
Figure 9. Proportion of soil types in different land use types (a); proportion of land use types in different soil types (b); load (c); and load intensity (d) of total phosphorus (TP) for each soil type between 2001 and 2010.

Figure 10 .
Figure 10.Load (a); and load intensity (b) of total phosphorus (TP) from each soil type between 2001 and 2010.

Figure 10 .
Figure 10.Load (a); and load intensity (b) of total phosphorus (TP) from each soil type between 2001 and 2010.

Figure 12 .
Figure 12.Monthly total phosphorus (TP) load from each: land use type (a); and soil type (b).

Table 1 .
Description of data available for the Sanjiang Plain.

Table 2 .
Estimated areas of each land use and land use changes from 2000 to 2010.

Table 2 .
Estimated areas of each land use and land use changes from 2000 to 2010.

Table 4 .
Statistical analysis of the observed value and the simulated value.

Table 4 .
Statistical analysis of the observed value and the simulated value.