Phosphorus and Nitrogen Yield Response Models for Dynamic Bio-Economic Optimization : An Empirical Approach

Nitrogen (N) and phosphorus (P) are both essential plant nutrients. However, their joint response to plant growth is seldom described by models. This study provides an approach for modeling the joint impact of inorganic N and P fertilization on crop production, considering the P supplied by the soil, which was approximated using the soil test P (STP). We developed yield response models for Finnish spring barley crops (Hordeum vulgare L.) for clay and coarse-textured soils by using existing extensive experimental datasets and nonlinear estimation techniques. Model selection was based on iterative elimination from a wide diversity of plausible model formulations. The Cobb−Douglas type model specification, consisting of multiplicative elements, performed well against independent validation data, suggesting that the key relationships that determine crop responses are captured by the models. The estimated models were extended to dynamic economic optimization of fertilization inputs. According to the results, a fair STP level should be maintained on both coarse-textured soils (9.9 mg L−1 a−1) and clay soils (3.9 mg L−1 a−1). For coarse soils, a higher steady-state P fertilization rate is required (21.7 kg ha−1 a−1) compared with clay soils (6.75 kg ha−1 a−1). The steady-state N fertilization rate was slightly higher for clay soils (102.4 kg ha−1 a−1) than for coarse soils (95.8 kg ha−1 a−1). This study shows that the iterative elimination of plausible functional forms is a suitable method for reducing the effects of structural uncertainty on model output and optimal fertilization decisions.


Introduction
The motivation for using nutrient inputs in agriculture more efficiently has never been greater than it is currently.This motivation originates from increased demand for food, finite available arable land and fertilizer material resources, expected increase in energy costs, growing public concerns related to environmental effects of agriculture, and the need for improved productivity and profitability of cropping systems [1,2].The changing climate and long-term socioeconomic trends are anticipated to alter the future challenges in agricultural production in complex, yet uncertain and possibly controversial ways [3][4][5].Sustained supply and long-term efficiency in the use of inorganic fertilizers, namely nitrogen (N) and phosphorus (P), will be critical for meeting the increasing global food demand [6].
Bio-economic models can be used as tools to study economically and environmentally sound fertilizer inputs in crop production.Multiple empirical models have been developed to separately capture yield response to N inputs [7][8][9], and to P inputs [10][11][12], and as a response to plant-available soil P [13][14][15].Both the response to P fertilizer and soil P dynamics have been incorporated into an economic analysis of nutrient inputs [16,17].Lambert et al. [18] studied spatially and temporally optimal fertilizer inputs for soybean and corn rotation.To the best of our knowledge, this is the only prior study that simultaneously considered the N and P inputs with P carryover in an economic analysis of crop production.The joint impact of N and P inputs on crop yield was explicitly considered in dynamic crop growth simulation models, such as DAISY [19], HYPE [20], APSIM [21], and EPIC [22,23], and in decision support systems for agrotechnology transfer like DSSAT [24].Although these models provide a rich and detailed description of the processes driving the yield response, they are not directly suited to dynamic economic optimization due to the large number of state variables and the extensive data requirements [25].
The existing literature investigating the joint impacts and dynamic feedbacks associated with economically optimal N and P fertilization is scarce, despite its importance for crop production.The obvious reason for the lack of this type of study is the need for extensive and long-term data from field experiments that capture the weather-induced variation in the impacts of N and P inputs and soil quality on crop growth [26].Such field experiments are expensive to conduct [27].Yield response modeling is a data-intensive process, as empirical data is needed for both building and validating the model [28].In Finland, a long tradition of N and P fertilizer field experiments exists, motivated by the challenging agro-climatic production conditions due to the Nordic location and natural scarcity of P in the soil [9,11,12,[29][30][31].As a result of these experiments, extensive empirical data for main cereals and different soil textures have been accumulated over five decades.However, identifying the effects of the individual nutrients from the compound fertilizer experimental data is difficult due to the complex interactions between N and P, and the general multicollinearity problem [32].
The objective of this study was to estimate and validate yield response models suitable for simultaneous optimization of inorganic P and N fertilization for spring barley (Hordeum vulgare L.).The models were developed for the following two soil texture groups: clay soils and coarse-textured mineral soils.The plant-available phosphorus dynamics (soil test phosphorus: STP) were included in the modeling to allow for the examination of temporally optimal fertilization paths.The models were designed for optimizing long-term farm level N and P fertilization decisions.The relative importance of different sources of uncertainty, including structural and parametric uncertainty, were evaluated.The robustness of various model specifications for economic optimization was also assessed.

Materials and Methods
In this study, an iterative approach to model building was applied [33][34][35].In the first phase (Figure 1), the conceptual model was developed and the available and relevant data for model estimation and validation were gathered.The second and the third phases included model estimation, verification, and validation.In the fourth phase, the yield response model was coupled with the transition function describing STP dynamics.In the fifth phase, the dynamic system model was applied for economic optimization.Finally, the sixth phase included uncertainty and sensitivity analyses.[36].STP, soil test phosphorus.

Phase I: Conceptual Model and Data
Our conceptual model was built on the following findings from the literature: (1) the amounts of nutrients in mineral fertilizers that should be applied depend on the nutrient requirement of a crop and the nutrient supply from the soil [2]; (2) inputs of inorganic P and N fertilizers have a direct positive effect on annual yield up to a certain level, as the availability of N and P are typically the most limiting factors for crop growth [37][38][39]; (3) P fertilization has a direct positive effect on STP due to the accumulation of P in the soil if the application exceeds the P removal by the harvest [40][41][42]; (4) N fertilization has an indirect negative effect on P balance by positively and indirectly affecting crop P uptake, meaning that the P is removed annually from the system in the harvested yield; (5) a positive association exists between STP and crop yield [14,15,43,44]; and (6) P fertilizer responses are lower for higher STP levels as P is available for plants from the existing reserve in soil.Consequently, yield responses to P fertilization are probable for STP classes ranging from poor to fair, whereas the responses become improbable for STP classes ranging from satisfactory to excessive [12,45,46].The applied ranges for the STP classes on coarse soils are here as follows: poor, 0-3; rather poor, 3-6; fair, 6-12; satisfactory, 12-20; good, 20-33; high, 33-50; excessive, >50 mg l −1 a −1 .On clay soils, the corresponding ranges are 0-2, 2-3.5, 3.5-7, 7-14, 14-23, 23-40, and >40 mg l −1 a −1 [47].STP was extracted from soil with acid ammonium acetate, pH 4.65 [48].
These findings were used to develop the conceptual farm-level dynamic decision model, where the planning horizon of the agricultural production, occurring on a representative field hectare, was divided into T stages at time t = 1, 2, . . ., T − 1 (Figure 2).The length of a stage is one year.In each stage, the producer decides the amount of N and P fertilizers used.Each stage t is characterized by a STP level, the state variable of the system.The system changes from one state at stage t to another at the next stage, t + 1, as described by the transition function ϑ.The transition function is a function of the current state of the system (STP t ), the fertilization decisions made at that stage, N t and P t , and the resulting output, i.e., the current crop yield y t : STP t+1 = ϑ(STP t , P t , y t (STP t , P t , N t )).At each stage, the decisions made and the state of the stage result in a stage return, R t = r t (y t (STP t , P t , N t ), P t , N t , p), where r t is a stage return function, which is a function of the current output, input decisions, and fertilizer and crop prices (denoted by a vector p).[49]), where t denotes a stage (one year), STP is the soil test P, which is the state of the system, P and N denote annual input decisions, R denotes the annual state return, r is the stage return function, y is the annual crop yield, p denotes the prices of inputs and output, and ϑ denotes a transition function that determines the evolution of the system from one stage to the next.
The data applied for the estimation of the yield response models included six separate datasets collected by Valkama et al. [9,12] that described the measurements from Finnish P and N fertilizer field experiments of cereals for clay and coarse-textured soils.The standardization of the separate datasets is shown in Appendix A. Additional datasets for specifying the impact of STP on crop yield without added P fertilization were obtained from Saarela et al. [11], since the information regarding this relationship was limited in the data obtained from Valkama et al. [12].Consequently, we were able to model the linkage between P response and STP, and the linkage between yield without added P and STP.However, the N fertilizer experiment datasets contained no information regarding the amount of N in the soil.Instead, the yields without added N were reported in the N datasets.Furthermore, the association between N response and the yields without added N was significantly negative (correlation = −0.63(p-value = 2.35 × 10 −8 ) for coarse soils and correlation = −0.49(p-value = 0.00027) for clay soils).The datasets were initially unbalanced panel datasets, implying that the datasets contained multiple observations for each cross-section unit in the sample, but the time dimension was not the same for each individual experiment.However, the observations in the applied datasets were averaged over the experimental years.Therefore, the observations from longer experiments were more reliable, since these data were less affected by random variation attributable to weather-driven events [50].Notably, the observations within the datasets were serially correlated, since multiple measurements were obtained from the same site.Moreover, for model validation, two independent datasets were gathered, including 28 short-term and long-term NPK fertilizer field experiments conducted in Finland between 1964 and 1988 at 10 sites on clay and coarse-textured soils.All the eight applied datasets are provided in the supplemental material.

Phase II: Model Development
In the applied datasets, the yield responses to individual nutrients were described as a relative increase in yield compared to the yield without the added nutrient.The Cobb−Douglas production function [51] is commonly used to describe such relationships.Our model specification resembled the Cobb−Douglas Mitschelich−Baule type production function [52] in that it included three multiplicative elements: where where y i denotes the crop yield (kg ha −1 ), i ∈ {1, . . . ,n} is a particular fertilization experiment, y P0 i denotes the yield without added P (kg ha −1 ), ω P i ∈ (0, ∞) is a scaling factor for P fertilization, ω N i ∈ (0, ∞) is a scaling factor for N fertilization, f l is a given nonlinear function specified by a certain structure for a model element l ∈ [1, 3], θ l,j indicates an unknown parameter with subscript j denoting a parameter number, ε l,i indicates the model residual error, and y N0 i denotes the yield without added N (kg ha −1 ).It is assumed that Equation (1) satisfies the typical properties of a production function [53].
The first element, defined by Equation ( 2), determines a yield without added P (kg ha −1 ) as a function of plant available soil P. The second element, defined by Equation (3), determines the yield response (%) to P fertilizer as a scaling factor that scales the first element up depending on the P application rate (kg P ha −1 ) and the STP level.The third element, defined by Equation ( 4), determines the yield response (%) to N fertilizer (kg N ha −1 ) as a scaling factor that scales the yield up or down depending on whether the N application rate is lower or higher than the average N rate in P experiments.The third model element is a decreasing function of yield without added N (kg ha −1 ).In model applications, y N0 was treated as a parameter rather than a variable due to the non-availability of a suitable transition function describing the N dynamics in the soil.The absence of a model component describing the N stock development is recognized as a model limitation.We hypothesized that a yield without added N must be lower than that without added P because typically N has a stronger effect on yields compared with P or STP [30,54].Therefore, y N0 was fixed to a "reference-level" of 2/3y P0 (2148 kg ha −1 and 2466 kg ha −1 for coarse-textured and clay soils, respectively).We used the weighted nonlinear least squares method for separately estimating the parameters for the particular functional forms used for f l in Equations ( 2)-( 4).We weighted the observations with the duration of the respective experiment.The heteroscedasticity and serial correlation of the data were considered by applying heteroscedasticity-and-autocorrelation-consistent (HAC) standard errors [55].For the N fertilization data on clay soils, however, the application of HAC standard errors was insufficient to prevent the estimates from being almost perfectly correlated.Therefore, in this case, we applied nonlinear mixed-effect (nlme) modeling, where the unobserved site-specific characteristics were treated as random effects and the data were grouped by the experimental plots.Estimations were based on the maximum likelihood in this approach.The impacts of the random effects were assumed to be asymptotically normally distributed and zero on average.The residual error included two mutually independent and additively separable parts.These are a random effect for group k = 1, . . ., K, and the random residual for individual response i of group k is as follows: u .Consequently, Equation (4) takes the following form for clay soils: ω N i = f 3 N i , y N0 i ; θ 3,j + z 3,i,k + u 3,i .Since the average expected yield is considered within the optimization and predictions, the site-specific random effects were set to their average expected level.More details about the estimation process can be obtained in Appendix A.
Typically, multiple models can be fitted successfully into data [8,56].Therefore, we examined several variations and combinations of commonly applied functional forms of f l in model elements in Equations ( 2)-(4) [35,57].As a result of this process, we obtained a set of candidate models.We chose four or five models out of the numerous iteratively obtained alternatives for each model element in Equation (1) for both soil textures.We ranked the estimated models using the modified second order variant of Akaike information criteria (AIC) [36].In addition, we calculated the relative likelihoods to quantitatively compare the models.We also performed a statistical goodness-of-fit analysis, as AIC is unable to indicate how well the models fit the data used for estimation [58].In particular, the minimum requirement for a model performance is that the Nash−Sutcliffe efficiency (NSE), defined as the ratio between the residual variance and the measured data variance [59], must be positive.Following the recommendations of Moriasi et al. [60], we used NSE as an overall indicator of model efficiency instead of the index of agreement (IA) by Willmott [61], because IA clearly overestimated the model performance.

Phase III: Verification and Validation
We explored the model behavior via various simulations.We also used graphical comparison, statistical goodness-of-fit analysis, and standard statistical tests to compare the simulated versus observed data.We tested the normality of the residuals between observations and predictions using the Shapiro−Wilk normality test.We also estimated the following simple regression model: where y is the observed yield and ŷ is the simulated model prediction, to evaluate the linear relationship between predictions and observations.In this study, we tested the following composite hypothesis: H0: α = 0 and β = 1 [62,63].Due to the serial correlation, following the suggestions of Marcus and Elias [63], HAC standard errors were again applied.Finally, we performed a paired t-test where the null hypothesis was that the mean of the observations is equal to the mean of the predictions, to quantify the differences between the observations and the predictions [62].

Phase IV: Formulation of the Dynamic System Model
To obtain the dynamic decision model (Figure 2), we coupled the estimated yield response models with the transition functions describing the soil P dynamics.We applied the following model [64]: where t indicates time (year), ∆STP t = STP t+1 − STP t defines a transition from state STP t to the next state STP t+1 , δ q with q ∈ [1, 4] as parameters, and P bal,t indicates the P balance of the soil at time t.
The model of Uusitalo et al. [64] was developed for similar conditions to those estimated in this study.
As such, the model required no calibration.We used the following model for the P balance [17]: where β 1 and β 2 are parameters, and the term (β 1 log(STP t ) + β 2 ) determines the P concentration of the crop yield.Thus, a P balance describes the difference between annually added P and removed P with the harvested yield, which is an increasing function of both P t and N t .

Phase V: Model Application for Economic Optimization
We formalized the economic optimization problem as a recursive finite horizon dynamic programming problem [65], where the objective is to maximize the net present value (NPV) of the discounted sum of the annual returns over the planning horizon by adjusting the control variables: subject to STP t+1 = ϑ(STP t , P t , y(STP t , P t , N t )) P t , N t , SPT t+1 ≥ 0 (10) STP 1 given (14) where β = (1 + ρ) −1 is a discount factor, ρ is a constant discount rate, r t is a state return function, p y is a constant price for yield, p P and p N denote constant prices for P and N, respectively, y t denotes an average expected yield at time t (Equation ( 13) is an expectation of Equation ( 1) for a given t), ∅ T denotes a scrap value [66], which approximates the following value of land after the production period: β T+t r T+t ≈ β T r T ρ , and STP 1 indicates the initial soil P level.This type of problem formulation is appropriate if we assume that the producers are risk-neutral, profit maximizers, and price takers.Furthermore, we assumed that farmers will own their farm for the decades ahead.Typically, farmers are interested in maximizing expected long-term profit if the farmer is a landowner and thus have a long-term planning horizon [67].We evaluated the suitability of the certain functional form for economic optimization by examining whether the following requirements were satisfied: the optimal control trajectory must converge to a finite steady state where the optimal solution for each stage remains approximately the same, and the optimal steady-state fertilizer rates, STP levels, and the yields do not exceed the corresponding sample maximums.

Phase VI: Uncertainty and Sensitivity Analysis
We evaluated the sensitivity of the optimization results to the structural uncertainty by solving Equations ( 8)-( 14) for different functions in constraint Equation ( 13) [8,57].We also calculated the differences in NPV by simulating the best models for input vectors obtained from other model candidates (Appendix A, Table A1) to approximate the potential economic losses resulting from incorrect decisions when selecting a yield response model [8].The model combinations applied in structural uncertainty analysis are described in Appendix A. Moreover, we examined the effects of parameter uncertainty on model outputs using Monte Carlo analysis.In the analysis, the probability distribution of all model parameters were assumed to asymptotically follow a normal distribution.We generated 10,000 values for each model parameter and computed the corresponding yields to obtain one realization for the yield distribution, which we described with the mean and the variance.To obtain the sampling distributions for the descriptive statistics of the yield distributions, we generated 10,000 realizations for the given statistic by looping the procedure.We also performed a sensitivity analysis to determine the sensitivity of the optimization results to the exogenous model parameters.We applied partial sensitivity analysis where the effect of each parameter on economic optimization results was examined while all the other parameters were held constant.

Model Estimation
After estimating numerous different model formulations for each model element in Equations ( 2)-( 4), by fitting various types of nonlinear functional forms for each dataset, and ranking the obtained model candidates using AIC (Appendix A, Table A1), we obtained the explicit form for the Cobb−Douglas type yield response model for coarse soils and for clay soils, respectively, as follows: where y is the average expected yield, STP is the soil test phosphorus, P is the phosphorus fertilizer, N is the nitrogen fertilizer, and θ l,j is a parameter with l ∈ [1, 3], indicating a model element, and j indicates a parameter number in the given model element.Parameter ω P,min denotes the minimum level of the scaling factor defining the P response.The numbers below Equations ( 15) and ( 16) refer to the model elements given in Equations ( 2)-( 4).Hence, the first model element in Equations ( 15) and ( 16) defines the average expected yield without added P as a function of STP (recall Equation (2): The second model element in Equations ( 15) and ( 16) defines the average expected yield response to P fertilization as a function of P and STP (recall Equation (3): ).The third model element in Equations ( 15) and ( 16) defines the average expected yield response to N fertilization as a function of N and yield without added N (recall Equation ( 4): Table 1 shows the estimated parameters and the summary statistics for each model element of Equations ( 15) and ( 16).All the parameters were statistically significant at the 5% significance level.The standard errors are heteroscedasticity-and-autocorrelation-consistent (HAC) for the coarse soils.In the case of first model element for the clay soils, the parameters were estimated via bootstrapping since the data sample was small.For the second model element, the standard errors were not HAC since the distribution of the residual errors was not significantly different from a normal one.Nonlinear mixed-effect (nlme) estimation was applied to the third model element.
For both soil textures, the ability of the models to capture the observed variability in the data was highest for the N response models, indicating that the addition of the variable yield, without added N, corrects most of the omitted variable bias related to the omission of the soil N. In particular, the addition of site-specific random effects mitigated the omitted variable bias.The high predictive accuracy of the N response model for clay soil resulted from the fact that most of the variation could be explained by the random effects.The average value of the random effects was 0.0015 across the 15 groups and the associated standard deviation was 0.13, whereas the standard deviation of the residual errors was 0.088.However, the nlme approach was not applied to the other models since it reduces the generalizability of the models.This became apparent by the poorer fit to the validation data.Notably, the use of independent datasets for validation is known to be a convenient way to avoid the problem of overfitting [68].The Nash-Sutcliffe efficiency (NSE) was also highest for the N response models for both soil textures (Table 2).In contrast, the discrepancy between observations and the model predictions was greatest for the model element describing the relationship between STP and yield without added P for both soil textures.The high prediction error was related to the relatively high variation in the dependent variable (yield) within the respective datasets, indicating that many other unobserved factors determining the variability in yields without added P, possibly more important than STP, exist.

Verification and Validation Results
Figure 3a,b show the simulated yield N-P-response surfaces for the fair STP class.The response surfaces are concave with respect to both nutrients.Moreover, the response surfaces for clay and coarse-textured soils were of the same order of magnitude, although for coarse soils, the observed P response (8.4% on average) was greater than that for clay soils (4.6% on average).For clay soils, in the other hand, the observed N response (74% on average) was greater than that for coarse soils (53% on average) (see footnotes for Tables S2 and S5 and for S3 and S6). Figure 3a,b also show that P fertilization does not affect the yield increases caused by N fertilization and vice versa.In this study, the yield responses of individual nutrients were estimated from the separate datasets; therefore, the interactions that are typically captured by the interaction term in polynomial function form could not be explicitly captured.Figure 3c,d illustrate the yield response curves for various STP levels and show a non-monotonic yield response to STP for a certain range of the domain for substantial P fertilization applications.The non-monotonicity results from the yield-increasing effect of P fertilization (modeled by the second model element) is considerable for low STP levels, whereas the effect vanishes for higher STP levels (Appendix A, Figure A1c,d).The yield in absolute terms, however, was an increasing function of STP (the relationship was modeled by the first model element).Therefore, the yield increased again for even higher STP levels.Thus, the interaction of the first and second model elements and the shapes of the respective response curves and surfaces (Appendix A, Figure A1a-d) were a result of the non-monotonicity observed in Figure 3c,d.
For both soil textures, the correlation between predicted and observed yield was reasonably high, suggesting that the model adequately represents the experimental results (Table 3).In addition, the relatively high NSE indicates that the variation in the observed yields was captured well, and the model simulation can be judged as "good" for coarse soils and "satisfactory" for clay soils [60].The statistics measuring the prediction accuracy indicated that the model for coarse soils had a smaller prediction error than that for clay soils, which is mainly due to the greater amount of variation in yields in the clay soil data compared to that in the coarse soil data.For both soils, the distribution of the residual errors was not significantly different from a normal distribution (Table 4).The estimation results of the simple regression described in Equation (5) showed that, for both soil textures, the α parameter was nonsignificant at the 5% significance level, whereas the estimate for the β-parameter was almost one and significant at the 5% significance level, indicating accurate predictive performance.The results of the paired t-test indicate that the mean of the observations was not significantly different from the mean of the predictions for either models, indicating that the absence of the interaction effect of P and N fertilization may be of minor importance for the performance of the models.This conclusion holds particularly well for the coarse-soils model.For clay soils, more observations were overestimated than underestimated, demonstrated by the relatively high mean bias, which was expected because the average yield in the validation dataset for clay soils was only 3307 kg/ha (Table S8).However, the systematic bias was minor (Figure 4e).The normality of the residual errors appears in Figure 4c,f.Moreover, the distribution of the residual errors was more normal in clay soils, as suggested by the Shapiro−Wilk normality test.

Economically Optimal Nutrient Inputs and Uncertainty Analysis
For comparison, Table 5 shows the steady-state solutions obtained by solving the economic optimization problem by dynamic programming and by simulating the model using the present maximum permissible N and P rates denoted by the voluntary Finnish agri-environmental scheme.Applying rates equal to or below these maximum permissible rates is necessary for receiving the payments that provide compensation for additional or opportunity costs or both resulting from obeying environmentally friendly farming practices.The comparison shows that, for both soil textures, the economically optimal steady-state N rates were close to the maximum permissible N rates when the exemplary prices are valid.Instead, the economically optimal steady-state P rate was 36% higher on coarse soils and 58% lower on clay soils than for the maximum permissible P rates.Consequently, the steady-state STP level for economically optimal fertilizer rates was slightly higher for coarse soils and slightly lower for clay soils than the steady-state STP level for the fertilizer rates set at the maximum permissible level.Despite these differences, losses in annual profits and NPV were marginal in both cases.However, the chosen prices and the discount rate affected these results.The examination of the sensitivity of the optimization solution to the input costs, output price, and discount rate is shown in Appendix A.
The statistics in Table 6 show that the average effect of the structural uncertainty is minor for both soil textures.The average economic loss of applying the "wrong" model, measured by the difference in NPV for the best model and any other model candidate, indicated by ∆NPV (€ ha −1 ), is also minor.The relative (mean-scaled) effect of the structural uncertainty on model output was greater for coarse soils.Regardless, the relative cost of structural uncertainty was greater for clay soils.
For coarse soils, most of the structural uncertainty effect was associated with the choice of the third model element (Figure 5a,b), which indicates that we were able to successfully fit function specifications, providing different response surfaces to the initial dataset.This can also be seen in Table 6, which shows that the relative effect (SD/Average) of the structural uncertainty was greater for optimal N fertilization (0.139) than for optimal P fertilization (0.07).Instead, for clay soils, the effect of structural uncertainty was mostly explained by the choice of the second model element (Figure 5c,d).Notably, Table 6 shows that, for clay soils, the relative effect of the structural uncertainty was greater for optimal P fertilization (0.41) than for optimal N fertilization (0.053).For both soil textures, the choice of the first model element had only a minor effect on the steady-state yield and NPV.However, the choice of the functional form affected the determination of the optimal steady-state yields (Figure 5a,c), and the incorrect choice may have economic consequences for the producer (Figure 5b,d).The exemplary fixed prices are 0.91€ kg −1 , 1.99€ kg −1 , and 0.12€ kg −1 , for N, P, and barley, respectively.The discount rate is 3.5% and the initial STP levels are 7.5 mg L −1 and 4.5 mg L −1 for coarse and clay soils, respectively.‡ STP indicates soil test P and NPV indicates net present value.§ The present maximum permissible N rate is 100 kg ha −1 and that of P is 34 kg ha −1 for the poor STP class, 26 kg ha −1 for the rather poor STP class, 16 kg ha −1 for the fair STP class, 10 kg ha −1 for the satisfactory STP class, 5 kg ha −1 for the good STP class, and 0 kg ha −1 for the high and excessive STP classes [69].The effect of parametric uncertainty on model output was greater for clay soils than for coarse soils (Table 7).For both soil textures, most parametric uncertainty stemmed from the third model element, the N response model.Moreover, the effect of parameter uncertainty on model output was greater than the effect of structural uncertainty.For coarse and clay soils the SD/Average ratios associated with the structural uncertainty were 0.037 and 0.027 for coarse and clay soils, respectively, whereas the corresponding ratios associated with the parametric uncertainty were 0.106 and 0.196 for coarse and clay soils, respectively., the values for N, P, and STP were fixed to their economically optimal steady state levels as shown in Table 5.The asymptotic 95% CI for the average yield and its standard deviations illustrate the uncertainty of the statistic.‡ The first model element refers to Equation ( 2), the second model element refers to Equation (3), and the third model element refers to Equation (4).
Lastly, we explored the effect of the initial STP level on optimization results using sensitivity analysis.The compared strategies were those where fertilization rates are adjusted to the initial STP levels by optimization and where the present maximum permissible rates of the voluntary agro-environmental scheme for each STP class are applied.When the initial STP level was below the steady-state level, applying a P rate above the steady-state level and a N rate below the steady-state level at the beginning of the planning horizon is optimal to rapidly achieve the steady state for STP (Figure 6).This resulted from the fact that, in the applied transition function for STP, the effect of N fertilization on STP development is negative: Instead, when we differentiated the transition function with respect to P fertilization, we obtained the following condition: In numerical applications, the effect of P on the development of STP is always positive (∂∆STP t /∂P t > 0), since β 1 = 0.000186, β 2 = 0.003 [17], and 0 < ∂y t /∂P t < 1. Conversely, when the initial STP level is above the steady-state level, there is no need for P fertilization for a long time.However, the STP level eventually declines below the steady-state level if no P fertilization is applied.To avoid this outcome, the optimal P fertilization rate thus increases steadily from zero to its steady-state level.Correspondingly, the optimal N rate is above the steady-state level at the beginning of the planning horizon and lowers to its steady-state level steadily.The convergence to the steady state was more gradual for clay soils than for coarse soils, since the N input is relatively more important for clay soils (Figure 6c,f).When the present maximum permissible rates were applied, the path for STP development was similar to the path for optimally adjusted fertilization rates (Figure 6a,d).Additionally, the P fertilization paths for the two strategies resembled one another (Figure 6b,e).For coarse soils, the optimally adjusted P fertilization rate for the poor STP class was considerably higher than that of the maximum permissible rate strategy.For clay soils, the P fertilization rates for all STP classes were higher for the maximum permissible rate strategy than those for the optimally adjusted rates strategy.The differences between the management strategies were greatest for the first 10 years of the planning horizon.The difference between the N fertilization paths was minor in absolute terms between the two management strategies for both soil textures, although the optimally adjusted N fertilization rates adjusted slowly to the steady state for clay soils (Figure 6c,f).
The initial STP had a relatively minor effect on the NPV for clay soils, with the difference between very low and high STP being approximately 1000€, compared to coarse soils where the difference was approximately 3000€ (Figure 7).Moreover, the application of the present maximum permissible rates led only to small economic losses.
Figure 7.The relationship between the initial soil test P (STP) and net present value (NPV) for the cases where the fertilization paths were determined by economic optimization and by applying the present maximum permissible rates (100 kg N ha −1 a −1 ) and declining maximum permissible P rates for the increasing STP classes.

Discussion
An advantage of using the existing extensive experimental data, available separately for N and P responses, is that we were able to avoid the multicollinearity problem.Another strength is that applying existing data is cheaper than conducting new field experiments [70].A disadvantage is that possible interactions in the joint effect of N and P fertilizer inputs were neglected due to the use of data from separate experiments for the two nutrients [38,39,71].Conversely, since the estimated models performed relatively well against the independent NPK fertilizer experimental data used for validation, the neglected interaction might be of less importance.Furthermore, the structure of the model and the modeling process described in this study can, after calibration and validation, generally be applied for any region, given that the required data are available.
Our results about the effect of the initial STP level on the optimal fertilization management showed that applying P on soils characterized by high STP levels according to the present interpretation is not practical, whereas for soils characterized by low STP levels, substantial P rates are economically justified.Similar findings were reported by Cox [72] and Valkama et al. [12].Maintaining fair or satisfactory STP levels has also been shown to require surplus P fertilization [15,29,73].In some cases, accumulating and maintaining a certain level of STP to obtain optimal crop yields may be profitable [44,74].According to our results, maintaining a fair STP level on both coarse-textured and clay soils was profitable.However, for coarse soils, this required higher P fertilizer inputs than clay soils, and the actual STP value was higher.This result shows the general effect of soil texture on the relation of P availability and STP for the examined specific region, and the generalization of the result requires exploration in other climatic regions.Thus, on coarse soils, the optimal fertilization strategy is more or less a so-called "build-up and maintenance" approach with a focus on the available amount of a nutrient in the soil.On clay soils, the optimal fertilization strategy follows more a "sufficiency" approach with a focus on the plant needs [2,44].For N fertilization, the rate of convergence to the steady state was rapid for coarse soils and gradual for clay soils, reflecting the relative importance of the nutrients.Somewhat similar findings were obtained by Lambert et al. [18].Thus, based on our results, the optimal steady-state fertilization rates and STP levels depend on the soil texture, among many other factors not included in the models.This result is consistent with suggestions provided by Syers et al. [44] and Reetz [2].
The developed models can be used to estimate the compensation farmers should be given for lowering the fertilization rates below the economically optimal rates.Our optimization results indicate that, although the present maximum permissible fertilization rates of the voluntary agro-environment scheme were somewhat different compared with the economically optimal fertilization rates determined by the models, the economic losses resulting from applying the maximum permissible rates were low.These results indicate that the present maximum fertilization rates applied in Finland successfully consider the effect of soil P dynamic feedback on optimal fertilization rates.For such cases, the compensation that farmers receive from applying fertilization rates that are different than the economically optimal rates could be minimal.
Previous studies showed that imperfect information regarding the model structure is often the main source of uncertainty associated with the model output [56,75].Our results of the uncertainty analysis showed the opposite: the effect of the parametric uncertainty on model output was greater for both soil textures than the effect of the structural uncertainty.However, we examined only the structural uncertainty for the mathematical formulations and ignored the structural uncertainty caused by the choice of processes included in the model.The effect of functional forms on optimal fertilization rates have often been shown to be significant [8,76,77].Our results for the effect of structural uncertainty on optimization results showed that the effect of functional form choice on optimal fertilizer rates was minor.Thus, via iterative elimination of plausible functional forms, by dropping the less suitable functional forms from further consideration in multiple phases, we efficiently reduced the effects of structural uncertainty on model output and optimal input decisions.
The results of the parametric uncertainty analysis showed that the parameter uncertainty had a greater impact on model output for clay soil models compared with coarse soil models.The impact of parametric uncertainty for clay soil models was associated with the greater variation in yield response to fertilization inputs within the initial data.This positive association between noisy data and model uncertainty has been recognized in previous literature [78][79][80].
Modeling exercises help detect areas where knowledge and data are scarce [81].We identified several issues that need additional research.First, given the large variation in observed yields without added P, particularly for clay soils, more long-term data on the relationship between yield without added P and STP, and the critical factors affecting P supply from soil to plants that may not be covered by STP measurements only, are needed [82,83].Secondly, more data and research are needed about the dynamics between soil N and soil organic matter (SOM), and the effect of these variables on crop yield.The weak association between SOM and a yield response was ignored in this study [9].Instead, the N status in the soil is quite volatile [84][85][86].Due to the data limitations and to consider the absence of soil N, our approach was to include a yield without added N into a model as an independent variable.Consequently, the N response models performed well against the initial N experimental data.Nevertheless, we recognized the absence of the N transition function from the economic model as a model limitation.To fully capture the dynamic behavior of the system, the soil N should be included in the economic model as a state variable.Hence, a natural extension for this study is the estimation of the transition function for N and the coupling of the transition function with the economic model.This extension, however, can be included in the model only once the data limitations are overcome.Moreover, the economic model could be expanded by introducing stochastic weather-related variables.Agriculture is strongly influenced by seasonal, weather-related, and climate-related uncertainty [87][88][89].Annual variation in crop yields is often explained better with annual weather-related variation in production conditions than with variation in fertilizer application rates [12,30,90].Nevertheless, when considering the long-term management of a farm, focusing on average expected yield is reasonable.Last, the explicit examination of the leaching potential remains a subject for further studies, since this would require coupling the system models with appropriate leaching functions for N and P.

Conclusions
The validation results suggest that the key relationships determining the crop responses were captured by the Cobb−Douglas type model specification consisting of multiplicative elements.The estimated models can be extended to assess the impacts of alternative management prescriptions and agricultural policies on the economic surplus of crop cultivation.This information may help improve long-term fertilization management plans.On coarse-textured soils, the optimal fertilization strategy can be described as build-up and maintenance, whereas on clay soils the optimal fertilization strategy is in line with the sufficiency rule of thumb.However, to fully capture the nutrient dynamics in bio-economic modeling, further research is required.Modeling approaches based on iterative elimination of plausible functional forms were efficient at reducing the effects of structural uncertainty on model output and on optimal input decisions.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2073-4395/8/4/41/s1,Table S1: Dataset for the first model element for coarse soils, Table S2: Dataset for the second model element for coarse soils.STP was measured at the beginning of experiments from the control (without added P), Table S3: Dataset for the third model element for coarse soils, Table S4: Dataset for the first model element for clay soils, Table S5: Dataset for the second model element for clay soils, Table S6: Dataset for the third model element for clay soils, Table S7: Validation dataset for coarse soils, Table S8, Validation dataset for clay soils.
STP level ranged from the poor to fair STP class.Instead, for clay soils, the yield without added P increased smoothly for higher STP levels.Figure A1c,d illustrate the associated predictive response surfaces for estimated best P response models.The yield response to P fertilization was observed only for STP classes ranging from poor to fair, whereas the response vanished rapidly for higher STP levels.The maximum P response was considerably higher for coarse soils than clay soils, indicating that the P fertilizer is a relatively more important input for coarse soils.Furthermore, the surface converged faster to its maximum plateau on clay soils when the P rate increased.Figure A1e,f illustrate the N response surfaces for the estimated models predicting lower yield responses to N fertilization for higher yields without added N. The maximum response to N fertilization was higher for clay soils, indicating that N fertilizer is a relatively more important input for clay soils.In addition, the coarse soil models converged to the maximum plateau for lower N rates.The illustrated models were the best model candidates.RSS is the sum of squared residuals, AICc is the second order variant of Akaike's information criteria, ∆ is the AIC difference, w is the Akaike's weight, and w 1 /w j is the evidence ratio between the best model (1) and another model j, and the parameter ω P,min is the minimum observation for yield response (0.96).‡ The 1. Model element refers to Equation (2) in the main text describing the association between soil test P (STP) and yield without added P (y P0i = f 1 STP i ; θ 1,j + ε 1,i ), 2. Model element refers to Equation (3) in the main text describing the yield response to P fertilization (ω P i = f 2 P i , STP i ; θ 2,j + ε 2,i ), and the 3. Model element refers to Equation (4) in the main text describing the yield response to N fertilization (ω N i = f 3 N i , y N0i ; θ 3,j + ε 3,i ).§ y P0 indicates expected average yield without added P, STP indicates soil test P, ω P indicates the scaling factor for the average expected P response, ω N indicates the scaling factor for the average expected N response, and y N0 indicates the yield without added N. effect on optimal fertilization rates.Furthermore, the effect of barley yield price on yield was greater for clay soils than for coarse soils, indicating greater elasticity in the yield price for clay soils.
Conversely, the effect of a discount rate was relatively minor compared with the other examined parameters.The higher the discount rate, the lower the optimal P fertilization rate for both soil textures.Instead, the optimal N fertilization rate was a decreasing function of the discount rate on coarse soils, whereas the optimal N fertilization rate was an increasing function of the discount rate on clay soils.This results from the fact that N fertilization was a more important factor for production on clay soils.Therefore, on clay soils, compensating for the marginal decrease in P applications caused by the increased discount rate by increasing the N fertilization rate is optimal.For coarse soils, the P fertilizer was a more important factor for production, so the reduction in P fertilization cannot be compensated by increasing applications of N.

Figure 1 .
Figure 1.Schematic diagram of the modeling process.AIC stands for Akaike information criteria [36].STP, soil test phosphorus.

Figure 2 .
Figure 2. The conceptual farm-level dynamic decision model (adapted from Hardaker et al.[49]), where t denotes a stage (one year), STP is the soil test P, which is the state of the system, P and N denote annual input decisions, R denotes the annual state return, r is the stage return function, y is the annual crop yield, p denotes the prices of inputs and output, and ϑ denotes a transition function that determines the evolution of the system from one stage to the next.

Figure 3 .
Figure 3. Above the simulated yield N-P-response surfaces for the fair STP class for (a) coarse and (b) clay soils.Below the yield as a function of STP for various P and N rates for (c) coarse and (d) clay soils.

Figure 4 .
Figure 4. Graphical model validation results for (a-c) coarse-textured soils and (d-f) clay soils.

Figure 5 .
Figure 5. Steady-state yields for (a) coarse-textured soils and (c) clay soils determined by the different model candidates, and the associated economic losses resulting from incorrect model specification measured in differences in net present value (NPV) for (b) coarse and (d) clay soils.The dashed line indicates the level associated with the best model candidate.

Figure 6 .
Figure 6.The effect the initial STP level on paths for STP development and for P and N fertilization for two fertilization management strategies.

Table 1 .
Parameters and the associated statistics for the models describing the yield response to soil test P (STP), P and N fertilization on coarse and clay soils as per Equations (15) and (16) † .

Table 2 .
Evaluation of the goodness-of-fit between the yield response models and the initial datasets † .Pearson's product-moment correlation coefficient between measured and simulated values is indicated by r, NSE is Nash−Sutcliffe efficiency, RMSE is root mean squared error, MB is the mean bias, and MAE is the mean absolute error.The statistics are a percentage of the observed mean values, in brackets.‡ Model Element (1) refers to Equation (2), Model Element (2) refers to Equation (3), and Model Element (3) refers to Equation (4).

Table 3 .
Evaluation of the goodness-of-fit between the model simulations and the observations in the validation datasets † .The mean-normalized statistics are shown in brackets.Pearson's product-moment correlation coefficient between measured and simulated values is indicated by r, NSE is Nash−Sutcliffe efficiency, RMSE is root mean squared error, MAE is the mean absolute error, and MB is the mean bias.

Table 4 .
Summary of the statistical tests for the model validation results.

Table 5 .
Steady states obtained by applying economically optimal fertilization rates and the maximum permissible fertilizer rates of the voluntary agro-environmental scheme in Finland † .

Table 6 .
Effect of the structural uncertainty on the robustness of the economic optimization results † .The number of the model candidates was 30 for coarse soils and 32 for clay soils.‡ STP indicates soil test P, NPV indicates net present value, and ∆NPV indicates the difference in NPV between the best model candidate and any other model candidate.

Table 7 .
Summary of the Monte Carlo analysis results used for examining the effect of the parametric uncertainty on model output (the barley yield), associated with individual model elements and the entire model for coarse and clay soils † .

Table A2 .
The sensitivity of the steady-state N and P fertilizer rates and yields to the economic parameters of the model.