Nutrient Load Mitigation with Wintertime Cover as Estimated by the INCA Model

Increased nutrient loading causes deterioration of receiving surface waters in areas of intensive agriculture. While nitrate and particulate phosphorus load can be efficiently controlled by reducing tillage frequency and increasing vegetation cover, many field studies have shown simultaneously increased loading of bioavailable phosphorus. In the latest phase of the Rural Programme of EU agri-environmental measures, the highest potential to reduce the nutrient loading to receiving waters were the maximum limits for fertilization of arable crops and retaining plant cover on fields with, e.g., no-till methods and uncultivated nature management fields. Due to the latter two measures, the area of vegetation cover has increased since 1995, suggesting clear effects on nutrient loading in the catchment scale as well. We modeled the effectiveness of agri-environmental measures to reduce phosphorus and nitrogen loads to waters and additionally tested the performance of the dynamic, process-based INCA-P (Integrated Nutrients in Catchments—Phosphorus) model to simulate P dynamics in an agricultural catchment. We concluded that INCA-P was able to simulate both fast (immediate) and slow (non-immediate) processes that influence P loading from catchments. Based on our model simulations, it was also evident that no-till methods had increased bioavailable P load to receiving waters, even though total P and total N loading were reduced.


Introduction
Increased nutrient loading of both phosphorus (P) and nitrogen (N) cause deterioration of receiving surface waters in areas of intensive agriculture. In Finland, agriculture is most intensive and is the dominant land cover on clayey soils situated in the southwestern part of the country from where runoff enters the Baltic Sea. In this region, a major proportion (60-90%) of total phosphorus (TP) load is transported as particulate phosphorus (PP) with eroded soil particles [1], and cultivated fields contribute 66-100% of the total suspended sediments delivered from a small agricultural catchment [2]. Moreover, past and present P excess fertilization is reflected in high soil test p values (STP) of the fields in the area, which is again connected with high dissolved reactive P (DRP) concentrations in runoff from fields [3,4].
Increased P loads, the major nutrient controlling eutrophication in many aquatic systems (e.g., [5]), have been observed to cause eutrophication of surface waters. The Baltic Sea may be seasonally or spatially N limited [6]. Further, P is not always the limiting nutrient in freshwaters. In Finland, some smaller lakes are observed to be partly N limited [7], calling for measures to reduce the loading of both P and N from agriculture.
Nitrogen (N) has partly different sources and pathways than P. Most commonly, it is transported as nitrate (NO 3 − ), which is highly leachable and can also percolate into ground water. Nitrate can reach receiving waters because of agricultural activities, such as N fertilization and tillage, enhancing the decay of crop residues and soil organic matter (SOM). Under anaerobic conditions, NO 3 − may be denitrified to N 2 . For cutting nutrient load from agriculture in Finland, the main tool since 1995 has been the Rural Programme of EU, whereby more than 90% of the farmers have joined and have been subsidized for using agri-environmental measures (AEM). In the latest phase of AEM (2014-2020), measures with the highest potential for reducing the nutrient loading on waters were upper limits for P and N fertilization of arable crops, winter surface cover of fields, and green set-aside as uncultivated nature management fields [8,9].
The most common measures to mitigate PP load to receiving waters by preventing erosion and suspended sediment transport include autumn ploughing or no-till that leaves an undisturbed soil surface over winter. For example, after changing over from intensive autumn tillage (ploughing) to permanent vegetation cover on clayey soils, erosion and PP loss decreased by 70-65% and 70-54%, respectively [10,11]. Moreover, N loads often decrease due to less intensive tillage [12].
While PP load can be clearly controlled by soil surface cover over winter, many field studies have shown simultaneously increased loads of DRP [11,[13][14][15][16]. One explanation is that the shallower tillage depth results in both stratification of P in the upper soil layers, (e.g., [17][18][19]) and, in Finnish conditions, higher risk of surface runoff [11,20]. The estimated P balance (P fertilization − P uptake in yield) can be 3 kg ha −1 higher in no-till fields than in cultivated fields, due to lower crop yield, which together with plant residues that are left on the soil surface gradually increase soil P status in the top layer and thus DRP load potential from the fields [13]. When tests have shown only around 5% (0-13%) of PP bioavailable in river waters but DRP totally bioavailable [21], their magnitudes and proportions are important to the actual ecological consequences of the different measures in place in an agricultural catchment.
In North America, Lake Erie is an example of severe water quality problems due to harmful algal blooms that have been connected to controversial effects of surface cover. Even though agri-environmental measures (no-till, buffer strips, etc.) have been able to decrease considerably suspended sediment and PP load from catchments to the lake, DRP loads have been increasing. Algal blooms are directly linked with increasing riverine input of DRP due to agricultural management and stratification of P at the soil surface [17,[22][23][24].
To estimate the efficiency of agricultural water protection measures, e.g., the success of AEM; water quantity and quality are frequently monitored in small catchments (e.g., [25]). Especially high frequency observations together with statistical tools can improve understanding of the processes that influence water quality [26]; however, that dense data is not always available at most locations, and/or time series are often too short.
Mathematical modeling can be useful to consider the combined effect of different pressures and water quality models are developed to determine the source, transformation, transport, and delivery of nutrients and suspended sediments in catchments and into water bodies. These models can be used to predict water quality in areas where monitoring is not feasible, or to help estimate the changes in water quality due to management actions. Compared to empirical models, process-based models can be more accurate to simulate conditions outside the current observations (e.g., [27]). Catchment-scale models, like INCA (Integrated Nutrients in Catchments) [28] and SWAT (Soil and Water Assessment Tool) [29], are especially valuable in assessing both water quality and quantity.
As catchments are heterogeneous systems and hydro-biogeochemical processes are non-linear, the use of dynamic, process-based models has also raised several questions concerning parameterization and calibration, such as equifinality and over-parameterization (e.g., [30]), and whether grid-based observations are valid in describing heterogeneous catchments [31][32][33][34]. Different sensitivity and uncertainty analyses are primarily concerned with the question of how model outputs are affected by variability in model parameters and input values and provide useful information when these components are not completely known. Sensitivity analysis [35] identifies which parameters predominately control model behavior, whereas uncertainty analysis relates prediction errors to uncertainty in the model structure and parameters. Good modeling practice requires that the modelers provide an evaluation of the confidence in the model, possibly assessing the uncertainties associated with the modeling process and with the outcome of the model itself (e.g., [36,37]).
Our main aim was to evaluate the effectiveness of AEM introduced by the Rural Programme to reduce P and N loads to receiving water courses using a chained set of catchment-scale models. Our motivation to concentrate on fields and agri-environmental measures is based on our previous study, where changes in general land use, animal density, point sources, or climate did not explain changes in water quality. Our second aim was to use sensitivity and uncertainty analysis methods to evaluate the performance of the dynamic, process-based INCA-P model for simulating P dynamics in an agricultural catchment against a large set of observational data at different levels.

Study Area
The Yläneenjoki catchment (233 km 2 ) has been one of the follow-up study areas of the effects of AEM (earlier FAEP) [8]. The catchment is located in the upper reaches of the Eurajoki River basin in southwestern Finland, as shown in Figure 1. The Eurajoki River discharges into the Bothnian Sea that is mostly P but intermittently N limited [6]. The mean annual precipitation in the area is 700 mm and the mean annual temperature is +4 • C (data from Finnish Meteorological Institute). With a mean discharge of 1.94 m 3 s −1 in 1991-2010 [38], the Yläneenjoki River flows to the Pyhäjärvi Lake, which is an important lake for water supply, commercial fishing, and recreation [39]. Eutrophication of the lake is ongoing. The Yläneenjoki River itself is estimated to achieve good ecological status by 2027 according to the river basin plan of the Water Framework Directive [40].
Main soil types in the catchment are clay (Eutric Cambisol 2, Vertic Cambisol) in the cultivated area and till and rock (Lithic Leptosol 1, Haplic Podzol 1, 2) in the forested area with some organic soil types (Fibric/Terric Histosol 1) [41]. Arable fields cover 28% of the catchment area, and the main crops are spring cereals with some sugar beet (Beta vulgaris L.) cultivation. The plant production is accompanied by animal husbandry, mainly pig and poultry, with animal density of 0.53 AU ha −1 in 2005 [42]. In the catchment, most of the farmers are committed to AEM included in the Rural Programme [42].

Catchment-Scale Data
Meteorological data, including daily measurements of precipitation and temperature, were obtained from the Finnish Meteorological Institute (FMI), Jokioinen Observatory station (WGS84 60.81, 23.5).
Maps of soil types, land cover, and the digital elevation model all had a grid size of 25 × 25 m. Soil data were obtained from [41] and land cover data were obtained from the CORINE 2006 GIS (geographical information system). Land cover data were supplemented by more detailed field parcel data of the Natural Resources Institute Finland.
The farmers were interviewed about their annual cultivation methods since 1995 up to 2005 [42][43][44]. Data between 2008 and 2016 originated from annual statistics of Natural Resources Institute Finland. Along with AEM, there has been an increase in the share of no-till methods and decrease in autumn ploughing as the primary tillage method of the spring cereal fields, as shown in Figure 2, while the share of reduced tillage (e.g., stubble cultivation) has remained under 20% of the field area throughout the study period, and the main changes have happened between no-till and ploughing. Measurements of river discharge (daily values) and suspended sediment, TP, total dissolved phosphorus (TDP), and DRP concentrations (grab samples: on average 33 per year in the Yläneenjoki River at Vanhakartano station) were available from national monitoring programs [45]. The main observation site (WGS84 60.87, 22.40) is located close to  Measurements of river discharge (daily values) and suspended sediment, TP, total dissolved phosphorus (TDP), and DRP concentrations (grab samples: on average 33 per year in the Yläneenjoki River at Vanhakartano station) were available from national monitoring programs [45]. The main observation site (WGS84 60.87, 22.40) is located close to the river outlet, as shown in Figure 1. In the upper reaches there was one automatic measurement station (WGS84 60.88, 22.58), which collected daily data of discharge and turbidity over five periods between 15 March 2006 and9 November 2007. In the analysis of TP, polyphosphates and organophosphorus compounds were first converted to orthophosphate by acid peroxodisulphate digestion under pressure at 120 • C. Orthophosphate ions reacted in an acidic solution containing molybdate and antimony ions to form a phosphomolybdate complex. Reduction of the complex with ascorbic acid formed a colored molybdenum blue complex, the absorption of which was measured at 880 nm by a spectrophotometer. TDP was similarly analyzed after first filtering the sample through a Nuclepore polycarbonate membrane (0.4 µm). Phosphate P was analyzed after filtering the sample and used as a measure of DRP. Particulate phosphorus (PP) was calculated as difference of TP and TDP. From TP, around 69% was PP, 19% DRP, and 12% dissolved unreactive phosphorus (estimated as the difference between TDP and DRP). Turbidity was used as an estimate of suspended sediment concentration at the upper measurement location due to their good correlation (y = 0.86 × x + 2.38, r 2 = 0.94). TP concentration correlated with turbidity (R 2 = 0.68) and TDP with DRP (R 2 = 0.78). TN determination was initiated by digestion with peroxodisulphate, followed by reduction of NO 3 − with a Cd amalgam and determination of NO 2 by the azo color method.

Field Scale Data
Field scale data with runoff from 2002-2007 (green set-aside) and 2008-2017 (crop cultivation) were obtained from the Kotkanoja experimental field of Natural Resources Institute situated on a clay soil (Protovertic Luvisol, Cleyic, Cutanic) in Jokioinen. In the study years, spring cereals were cultivated in the field with either continuous no-till (two 0.5 ha plots) or autumn ploughing (two 0.5 ha plots) as primary tillage and winter cover options. The field had a gentle slope of 2% (1-4%) and a well-functioning subsurface drainage system (the averaged share of sub-surface flow of total measured runoff, sub-surface plus surface runoff, was 63 and 75% in no-till and ploughing, respectively). Quantities of both the surface runoff and the sub-surface drainage flow were recorded throughout the years by tipping buckets and the flow-weighted concentrations of TP, DRP, and TN were measured using similar methods as described above with the exception of using a filter pore size of 0.2 µm for DRP analyses (partly published for the above time interval, [11]).

Models and Model Set-Ups
We used a modification of the Universal Soil Loss Equation (USLE) [46] to calculate the relative erosion sensitivity for each of the fields in the catchment based on land use, soil type, and topography. In the modified version RUSLE [47], spatial differences in climate are not considered, and the range of the factors has been adjusted for Finnish conditions. The relative erosion sensitivity (K) is calculated as a product of four constant factors, determined by soil type (E), slope (R), land use (M), and distance to water (V) as

PERSiST
We used a flexible rainfall-runoff modeling toolkit, PERSiST (the Precipitation, Evapotranspiration, and Runoff Simulator for Solute Transport), with the INCA family of models. PERSiST is designed for simulating present-day hydrology, projecting possible future effects of climate or land use change on runoff and catchment water storage [48]. PERSiST has limited data requirements and is calibrated using observed time series of precipitation, air temperature and runoff, at one or more points in a river network.

INCA-N and INCA-P
INCA-N [27,49] and INCA-P [50] are dynamic mass-balance models, which attempt to track the temporal variations in the hydrological flow paths and nutrient transformations and stores, in both the land and in-stream components of a river system. INCA models provide daily and annual land-use specific nutrient loads for all stores and transformation processes, and include concentrations in the soil water, ground water and direct runoff. Decay of the soil organic N pool is described by first-order kinetics, and the decay rate is limited by temperature and moisture. The change in mass of P adsorbed to or desorbed from the soil is calculated using a formulation of the Freundlich isotherm [51], where the change in mass adsorbed is proportional to the difference between TDP concentration in soil water and the equilibrium TDP concentration at which no adsorption or desorption occurs (EPC 0,soil ). It is a dynamic variable that depends on the mass of P in the labile soil store. The other soil store is the inactive P store. This allows for the simulation of longer term soil P dynamics and fits with monitoring data showing EPC 0,soil increasing as soil P enrichment increases (e.g., [52]).

PEST
The Model-Independent Parameter Estimation system (PEST; [53]) was used to make a parameter sensitivity and uncertainty analysis. PEST is a non-linear parameter estimation tool for automatic model calibration and predictive analysis. It aims to minimize the weighted sum of squared differences between the model simulation and an observed set of data. The optimization problem is iteratively solved by linearizing the relationship between model output and parameters using a Taylor series expansion. The parameter vector is updated at each step using the Gauss-Marquardt-Levenberg algorithm [54,55]. The derivatives of the model outputs with respect to its parameters provide a measure of parameter sensitivity. Uncertainty analysis is conducted by keeping the objective function with an aim to minimize or maximize a certain set of observations.

Model Set-Up
We used the branching version 1.4.11 of the INCA models in this study. To reduce the number of parameters and thus the risk for over parameterization, and to enable multiple model runs, the calibration was simplified as much as possible, while still aiming to keep the characters of the system.
The INCA-P model was calibrated against daily observations of discharge and monthly to biweekly observations of suspended sediments, and TDP and TP concentrations observed close to outlet. There were also available short series of automated discharge and turbidity measurements in the upper sub-catchment. The calibration period was 2003-2008 and the validation period 2009-2011. The INCA-N model, in turn, has previously been calibrated against NO 3 -N and NH 4 -N concentrations in the same river [56].
The Yläneenjoki catchment was divided into two parts according to the erosion sensitivity map based on the modified USLE. In the less sensitive part of the study area, soil types are typically not prone to erosion and fields are not directly connected to the river. In the highly sensitive area, soil types are prone to erosion and fields are located on steep slopes and are typically also connected to surface water courses [57]. These two areas have different parameter values for erosion and transport capacity processes. In the Yläneenjoki catchment, 40% of the field area was allocated in the less sensitive area. In addition, the upper reach where the automatic measurement station was located was split to form a separate sub-catchment. Moreover, the river was divided into 3 sections according to channel width.
The average fertilizer use, crop yield, P uptake, and initial soil-test P (STP) were derived from farmer interviews (e.g., [42]). Land use was divided into three classes. The most erosion sensitive class included fields where spring cereals were cultivated on clay soils. Grass fields are less sensitive to erosion due to a longer growing period and denser vegetation cover. Sugar beet was cultivated on coarser soil types than spring cereals, and the fields do not have permanent or long-lasting vegetation cover. In addition, P balances that affect STP tend to be high for the sugar beet fields. P loads from forests were calibrated against observed loads [58] and P equilibrium [59]. River water and sediment equilibrium is based on the literature [60], but no observations of P exchange between sediment and water in the river were available.
In the PEST analysis conducted against suspended sediment and TP concentrations, three groups of parameters were used. Furthermore, the parameter ranges, initial parameter values, and parameter increments must be provided for PEST. Parameters with known values, e.g., catchment area and areal percentages of different crops, belonged to the first group. These values were not allowed to change. The second group included relatively well-defined parameters for which the range of values was set narrow and to have known limits. The third group includes parameters for which the values were not well-defined, and the range was set to either +/−25% or +/−50%.
Kotkanoja research field measurements were used to calibrate sub models and equations in the INCA models. Nutrient concentrations in single land use classes (crop-soil type-tillage method combination) were calibrated against concentrations measured in the research field. Equilibrium coefficients were calibrated against those obtained for different soil profile layers in the research field. In addition, the calibrated equilibrium coefficients were compared to measured equilibrium coefficients in corresponding laboratory trials [61].
We used different goodness-of-fit indices. Pearson's coefficient of determination (R 2 ) is an index of the degree of linear relationship between observed and simulated data, so that value 0 indicates no linear relationship. The Nash-Sutcliffe efficiency (NSE) is commonly used in hydrology. It is a normalized statistic that determines the relative magnitude of the residual variance compared to the measured data variance [62]. Values lower than 0.0 indicate that the mean observed value is a better predictor than the simulated value. RMSE (Root Mean Square Error) is one of the commonly used error index statistics, and it is commonly accepted that the lower the RMSE, the better the model performance. Percent bias (PBIAS) measures the average tendency of the simulated data to be larger or smaller than observations. The optimal value of PBIAS is 0.0, with low-magnitude values indicating accurate model simulation.

Scenarios
We created three scenarios that were all run for a 22-year-long period using the same initial conditions and climatic data but different agricultural land uses:

1.
Autumn ploughing as the primary tillage for all spring cereal fields in the Yläneenjoki catchment; 2.
Primary tillage as actually performed in the catchment in 1994-1995; 3.
Primary tillage as actually performed in the catchment in 2008-2010.
Scenario 1 was set as the baseline, and the two other tillage scenarios were compared to the baseline.

Calibration and Uncertainty Bounds
Split-sample calibration and validation periods cover years 1995-2007, when exact cultivation practices in the area were available. The starting point for sensitivity and uncertainty analysis was in turn the manually calibrated model for the period 2003-2007, when automatic data was available at the upper reaches. Manual calibration was further finetuned by automatic PEST calibration. The validation period was 1995-2002. Goodnessof-fit values R 2 , RMSE, percentage bias (PBIAS), and Nash-Sutcliffe efficiency (NSE) for discharge are presented in Table 1. For phosphorus validation, measures with R 2 ≤ 0.40 were considered as not satisfactory, whereas for nitrogen 0.30 < R 2 ≤ 0.60 were satisfactory and 0.60 < R 2 ≤ 0.70 were good. Values for PBIAS ≥ |30| were considered not satisfactory, |20| ≤ PBIAS < |30| satisfactory, |15| ≤ PBIAS < |20| good, and PBIAS < |15| very good in their accuracy according to published guidelines [63,64].
Discharge calibration was good based on NSE. Nitrate calibration was satisfactory, and TP and DRP calibration were not satisfactory. During the validation period, goodness-of-fit values of discharge and nutrients decreased considerably along the years. On the other hand, PBIAS of discharge was <10 during both calibration and validation periods, which resulted in "very good" calibration. In addition, NO 3 − simulations were very good, and TP and DRP simulations satisfactory. Simulated TP and DRP concentrations reached the same level and seasonality as the observations, as shown in Figures 3 and 4. There was greater variation in the observed than in the simulated concentrations, but that variation fell between 95% uncertainty bounds, as shown in Figure 5.

Plot and Field Measurements
In the INCA-P soil compartment, the labile P pool corresponds to the easily soluble P pool from which soluble P leaching occurs. Soil-test P values can thus be used as indicators of the size of that pool. The simulated influence of soil P balance (the difference between the input of P in fertilizers and output in crop uptake) on the labile P pool was compared with the results of an empirical model [61], which related the P balance surplus (or deficit) in a farm to the edge-of-field loads of algal-available P. Based on long-term fertilizer trials, the model first estimated the change in soil-test P of the topsoil with the

Plot and Field Measurements
In the INCA-P soil compartment, the labile P pool corresponds to the easily soluble P pool from which soluble P leaching occurs. Soil-test P values can thus be used as indicators of the size of that pool. The simulated influence of soil P balance (the difference between the input of P in fertilizers and output in crop uptake) on the labile P pool was compared with the results of an empirical model [61], which related the P balance surplus (or deficit) in a farm to the edge-of-field loads of algal-available P. Based on long-term fertilizer trials, the model first estimated the change in soil-test P of the topsoil with the aid of the soil-surface P balance. Soil-test P was then used to approximate the concentration of DRP in surface runoff and drainage flow, as adjusted for different P application types. The INCA-P model produced slightly faster change in soil labile P, as indicated by STP, than the empirical model. Even though in the Yläneenjoki system the simulated DRP concentration was lower, the slope of the change was the same, as shown in Figure 6. aid of the soil-surface P balance. Soil-test P was then used to approximate the concentration of DRP in surface runoff and drainage flow, as adjusted for different P application types. The INCA-P model produced slightly faster change in soil labile P, as indicated by STP, than the empirical model. Even though in the Yläneenjoki system the simulated DRP concentration was lower, the slope of the change was the same, as shown in Figure 6.     Figure 6. The effect of P balance on soil-test P (STP, based on the acid ammonium acetate method, pH 4.65, used in Finland, [62]) values (a) and DRP concentration (b) in river water.
Simulated soil water TP, DRP, and NO3 − loads corresponded to the calculated loads based on measured drainage water and surface runoff concentrations, as shown in  Simulated equilibrium isotherms for P were close to what [61] measured in agricultural soils (Figure 9), as shown in Figure 7, and in forest soils, the isotherm values were close to what [59] measured in Podzol soils in Myrtillus type forests. Further, equilibrium concentrations of P in river bottom sediments and in water were in line with measurements of [60]. Cross stands for mean and dots for outliers. Simulated equilibrium isotherms for P were close to what [61] measured in agricultural soils (Figure 9), as shown in Figure 7, and in forest soils, the isotherm values were close to what [59] measured in Podzol soils in Myrtillus type forests. Further, equilibrium concentrations of P in river bottom sediments and in water were in line with measurements of [60]. Cross stands for mean and dots for outliers. Figure 9. Equilibrium isotherms of P in soil as modeled and measured in the Kotkanoja research field at different soil layers.

Parameter Sensitivity
The most sensitive parameters were those describing flow velocity in the river (flow) and catchment and soil sensitivity to erosion, as shown in Figure 10. The parameters describing fertilization, vegetation growth, or P processes in the soil were not among the 20 Figure 9. Equilibrium isotherms of P in soil as modeled and measured in the Kotkanoja research field at different soil layers.

Parameter Sensitivity
The most sensitive parameters were those describing flow velocity in the river (flow), and catchment and soil sensitivity to erosion, as shown in Figure 10. The parameters describing fertilization, vegetation growth, or P processes in the soil were not among the 20 most sensitive parameters. Figure 9. Equilibrium isotherms of P in soil as modeled and measured in the Kotkanoja research field at different soil layers.

Parameter Sensitivity
The most sensitive parameters were those describing flow velocity in the river (flow), and catchment and soil sensitivity to erosion, as shown in Figure 10. The parameters describing fertilization, vegetation growth, or P processes in the soil were not among the 20 most sensitive parameters.

Scenario Results
According to the scenarios, change in tillage methods from autumn ploughing towards no-till methods decrease nitrate and suspended sediment loading from the Yläneenjoki catchment but increase DRP loading compared to the baseline with autumn ploughing only, as shown in Table 2. Similar change is seen in land use specific loads, as shown in Table 3, which are in the range of literature values for the corresponding cultivation methods [12].

Model Development
We used a set of goodness-of-fit measures to evaluate the success of the model performance, as they measure different features in the fit between simulated and observed values. Pearson's coefficient of determination (R 2 ) has been widely used for model evaluation; however, it is oversensitive to high extreme values (outliers) and insensitive to additive and proportional differences between model predictions and measured data [65]. NSE is very commonly used in evaluating the fit of a hydrograph, which provides extensive information on reported values. PBIAS has the ability to clearly indicate poor model performance [66]. Even though R 2 values in our study were low, BIAS gave relatively good fit between simulated and observed values. That was confirmed visually in a figure, where observed values lay between simulated uncertainty bounds.
Catchment models offer useful methods for exploring problems in space and time, but the quantity and quality of data are not sufficient to parameterize the models uniquely and the modeling addresses a wide range of uncertainties [34]. In addition, several empirical studies have shown that many models and many parameter combinations give equally good fits to data, indicating that it is impossible to find an optimal model or an optimal parameter set in hydrological modeling-a problem termed equifinality [67]. Thus, the use of experimental knowledge of nutrient processes, soft data, is recommended during calibration as it helps to constrain the model based on the available knowledge of the catchment processes [68]. In this study, we used data from experimental fields to calibrate P processes and their outcomes to observed levels.
Various methods from expert assessment to sensitivity analysis can be applied to evaluate uncertainty of deterministic model outputs, concluding that the best method for uncertainty evaluation is determined by the definition of a model, and the amount of available information [69]. A sensitivity analysis framework was recommended as a good starting point to consider uncertainty of spatially distributed environmental models in general [70]. In our study, the model outcome was most sensitive to parameters that defined river flow and catchment erosion processes. These are the processes that give the immediate response to simulated discharge and TP concentrations. The values of these parameters are easy to either measure or calibrate against the measurements. The model was not sensitive to parameters defining DRP processes, reflecting that these processes do not have immediate effect on P concentrations. Thus, change in environmental conditions/changes that influence these processes are difficult to catch in short-term simulations, or see in calibration.
Most of the P is transported to receiving waters attached to erosion material, and thus most of the P transport models are concentrated on describing these processes accurately. When modeling the fate of DRP, equilibrium isotherms both describe the dynamic, interactive P transport between the soil labile P pool and soil water, and can be defined in the laboratory for different land cover types. As equilibrium equations add complexity to modeling process, most of the models use more simple descriptions, either by using empirical equations that are based on soil physical properties, or by limiting the active land cover types (e.g., [30,71,72]). The interpretation of the results should thus be proportionate to how well the model describes the relevant processes.
The Yläneenjoki catchment is an erosion-sensitive area due to the high elevation differences and erosion-prone clay soil types [73]. As such, it is a typical example of an area where TP and suspended sediment load correlate. Even though field percentage is not exceptionally high, around 70% of TP load comes from fields [74]. The share of forests is presumably reflected in water quality, so that TDP concentration is higher than DRP concentration (assumedly P in organic matter). There have not been major changes or trends in anthropogenic pressures on forests, so we assume that changes in water quality are related to agriculture. In the previous study [57], general changes in land use, number of domestic animals, or climate did not explain changes in water quality.
We concluded that INCA-P was able to describe both fast, immediate processes but also slow, non-immediate processes that influence P loading from catchments. We recommend paying special attention when calibrating slow processes controlling P pools in the soil, because they are easily missed when focusing on fast, immediate processes in suspended sediment and PP loading.

Nutrient Load Mitigation by Winter Cover
Agri-environmental measures (AEM) of the Rural Programme have shown success in improving nutrient balances by introducing maximum fertilization levels, and in increasing the area of reduced tillage and wintertime vegetation cover [8]. However, while previous studies have shown that lower fertilization has the potential to reduce nutrient loads (e.g., [75,76]), the observed decrease in riverine nutrient loads has been modest [74]. In controlling the phosphorus loads, AEM have so far succeeded only in catchments with erosion-sensitive soils where the decrease in suspended sediment and PP loads exceeded the increase in DRP load [56].
According to this study in the Yläneenjoki catchment, the percentage increase in simulated algal available DRP load was higher than the percentage decrease in TP load when the land use in 2008-2010 statistics was compared to the baseline (autumn ploughing in spring cereal fields). Thus, according to [20], this change in tillage was enough to increase bioavailable P load to the Pyhäjärvi Lake up to 10%. No-till methods were relatively popular already in 1994, before launching AEM, and that change was enough to increase bioavailable P load by 0-2%. Statistics from different years are not directly comparable as they originate from different sources and have different resolutions. In any case, we can conclude, based on our model simulations, that no-till methods probably have a contribution to increased bioavailable P load to the Pyhäjärvi Lake.
This study reflects the contradictory nature of vegetation cover and especially notill methods. When planning catchment-scale water protection methods, the tradeoff is between N, P, and suspended sediment loads, at least when there is no such favorable system change in, e.g., soil structure due to the introduced methods that would clearly influence the result in the long term, like the reduced proportion of surface runoff or the lower need of P fertilizer input for good crop growth. In our research field data [11], the immediate increasing effect induced by green set-aside (that is very similar to the management of the uncultivated nature management field of AEM) on DRP loss was clearly less, while the decreasing effect on TP loss was even more than for no-till.
If the focus is on N-limited parts of the Baltic Sea, no-till methods are one possible measure to reduce N loading. If the focus is on P-limited lakes or sea areas, other methods to decrease DRP loads should also be considered. Even though nutrient balances in Finland have decreased [8], input of P in mineral fertilizers and manure is still above the need of plants, if soil reserves are accounted for [77], and P fertilization needs to be further optimized to speed up the decline of STP values on fields where they are unnecessary high for good crop growth. By this means, it is possible to achieve a permanent reduction in the P loading potential of fields. Moreover, some researchers (e.g., [17]) have even recommended that no-till fields should be occasionally tilled or ploughed conventionally to cut down the stratification of P in the top layers. Obviously, the controversial effects of wintertime cover on nutrient loading to surface waters call for advanced decision tools for site-specific and optimized soil management planning.

Conclusions
An equilibrium isotherm describes the dynamic, interactive P transport between the soil labile P pool and soil water in the INCA-P model. We calibrated successfully these equations against a long-term field trial. The model outcome was still the most sensitive to parameters that defined river flow and catchment erosion processes. These are the processes that give the immediate response to simulated discharge and TP concentrations.
We concluded that INCA-P was able to describe both fast, immediate processes but also slow, non-immediate processes that influence P loading from catchments.
According to this study in the Yläneenjoki catchment, the percentage increase in simulated algal available DRP load was higher than the percentage decrease in TP load, when the land use in 2008-2010 statistics was compared to the baseline (autumn ploughing in spring cereal fields). This change in tillage was enough to increase bioavailable P load to the Pyhäjärvi Lake up to 10%. We can conclude based on our model simulations, that no-till methods probably have a contribution to increased bioavailable P load to the Pyhäjärvi Lake.  Data Availability Statement: Publicly available discharge and river water quality datasets were analyzed in this study. This data can be found here: https://www.syke.fi/avointieto (accessed on 5 September 2018). Statistical data is available at https://www.luke.fi/avoin-tieto/tilastopalvelu/ (accessed on 5 September 2018). The field scale data presented in this study are available on request from the corresponding author. The data are not publicly available due to content of personal information.

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