Modelling Changes in Soil Phosphorus When Phosphorus Fertiliser Is Reduced or Ceases

In order for land managers and policy makers to manage excessive soil phosphorus (P) concentrations and reduce the risk of this particular source of P from impacting water bodies, models of soil P decline under various scenarios are needed. We modelled the decrease in calcium chloride-extractable P (CaCl2-P), and sodium bicarbonate-extractable P (Olsen-P and Colwell-P) using data from six Australian grazed pasture soils with contrasting P sorption properties, over a period of 4.5 years. Each soil had four initial soil P concentrations (Pinit), each of which received four on-going rates of P fertiliser (Pfert). The model predicts the final P concentration (Pfinal) by taking into account the P concentration previously measured (CaCl2-P, Olsen-P or Colwell-P), Pfert applied since measurement, and time since previous measurement: Final P concentration = (previously measured P concentration + ep x P fertiliser applied) exp (-dp x years since previous P concentration measurement). Where ep is the increase in soil P for each unit of applied P and dp is the decay constant representing how quickly the soil P decreased. The greatest decreases in proportion to Pinit occurred for CaCl2-P, followed by Olsen-P, and then Colwell-P. The model tended to fit the dataset well for Olsen-P and Colwell-P, with mean overestimation (modelled Pfinal concentration greater than actual Pfinal) of the Pfinal concentrations of 6.1 (32%) and 4.3 mg/kg (10%), respectively. Although there was less CaCl2-P data, the model successfully described it, with a mean overestimation of Pfinal CaCl2 of 3.1 mg/kg (26%). The overestimation of Pfinal CaCl2 was possibly due to the high CaCl2-P concentrations of the low P buffering index soils. The model predicted an average of 32 years (ranging from 26 to 49 years) for Olsen-P concentrations of between 55 and 96 mg/kg to decrease to an agronomic optimum of 17 mg/kg. Agronomic optimum was not a reliable indicator of environmental risk as some soils did not exceed the CaCl2-P environmental threshold until Olsen-P concentrations were twice the agronomic optimum, whereas low P sorbing soils tended to exceed the threshold before reaching agronomic optimum. Further work with more soils is required to examine the influence of soil properties – such as P sorption – on decreases in soil P.


INTRODUCTION
Excessive soil phosphorus (P) concentrations are a concern in parts of Europe, the United States of America and Australia (Paulter and Sims, 2000;Sharpley et al., 2001;Gourley et al., 2015). Soils with more P than is needed for agronomic purposes can result in unnecessarily high concentrations of P in runoff (Sharpley and Rekolainen, 1997;Burkitt et al., 2010) which can, in addition to other landscape and land management factors, contribute to eutrophication of waterways (Carpenter et al., 1998). These excessive soil P concentrations can be the result of P inputs (fertiliser, livestock feeds and manures) exceeding P removal in produce. Some land managers still apply P fertiliser when soil P concentrations are above agronomic optimum, despite being aware of the environmental risks, especially in intensive pasture industries (Nicon Rural Services, 2010). This practice is used as a form of risk management to ensure that soil P does not suddenly decrease below the agronomic optimum and result in a subsequent decrease in pasture production and consequent economic loss. However, this approach results in inefficient use of P fertiliser, a finite resource which is in increasing world demand (Smil, 2000;Condron et al., 2005).
In response to increasing P concentrations in some water resources, policies have been implemented in some countries to prevent further increases in soil P concentrations that are already above agronomic optimum (European Parliament, 2000). Models predicting increases in plant available P concentrations are well established, however, few models predict declines. In order for land managers and policy makers to better understand how to manage soil P so as to decrease excessive concentrations without risking production losses, models of the changes in soil Pparticularly decreases -under various scenarios are needed.
A number of researchers have described an exponential decrease in extractable soil P after P fertiliser application has ceased (McCollum, 1991;Paris et al., 2004;van der Salm et al., 2009). A number of mechanistic (Russell, 1977) and descriptive (Cox et al., 1981;Schulte et al., 2010;Dodd et al., 2012) models have been devised to model soil P decline. Schulte et al. (2010) and Dodd et al. (2012) developed models across a range of initial P (P init ) concentrations and soils, with Schulte et al. (2010) also presenting a number of negative P-balance scenarios. However, these models did not include the influence of different ongoing P fertiliser (P fert ) rates on decreases in extractable soil P, a practice which commonly occurs for intensive pasture soils and agricultural enterprises more generally. Coad et al. (2014) described changes in extractable soil P in a range of field soils which each had a wide range of initial soil P concentrations and on-going fertiliser P rates.
In this current paper, we use a Bayesian statistical framework to describe the changes in environmental P (CaCl 2 -P) (Houba et al., 2000) and agronomic P (Olsen-P and Colwell-P) (Olsen et al., 1954;Colwell, 1963) using datasets from the field study of Coad et al. (2014). We use a Bayesian approach to allow for uncertainty in measurement and parameters to be incorporated in a natural way through the appropriate prior specification. Bayesian methods also allow complex models to be more easily constructed (Brooks, 2003). Our aim was to develop a model to predict the final extractable soil P concentration (P final ) after a specified time period for various soil and P fertiliser scenarios. This model could be used by policy makers and land managers to assist their decision making and to set realistic timeframes to return soils that have excessive soil P to agronomic optimum concentrations and/or concentrations that pose low environmental risk.

Site Management and Soil Characterisation
The materials and methods employed to obtain the dataset used to develop the Bayesian model are published in Coad et al. (2014). Briefly, six field sites were selected that had a wide range of P sorption capacities (Table 1). Subsequently a wide range of P rates were applied as a single initial application to generate varying initial soil P concentrations. Then P was applied every 6 months at rates of zero, 0.5, 1, and 2 times maintenance rates for up to 4.5 years. The experimental design was fully factorial with 3 replicates at each site. The sites were sampled regularly for up to 4.5 years. Agronomic (Colwell and Olsen) measures of P and an environmental indicator (0.01 M CaCl 2− P) of P status were measured on soil samples. We also measured both unadjusted P buffering index (PBI) (Burkitt et al., 2008) -commonly used in Australia -and the amount of P sorbed after addition of 1000 mg P/kg (Ps 1000 ) for the initial and final soil samples for all treatments. Both PBI and Ps 1000 were undertaken using the PBI method (Burkitt et al., 2008;Rayment and Lyons, 2011). Briefly, 1 g of soil was shaken for 17 h in 0.01 M CaCl 2 with an addition of 1000 mg P/kg soil (Burkitt et al., 2008). The difference between the added P and the final P concentration in the equilibrating solution was the amount of P sorbed (Ps 1000 ), with PBI calculated as P sorbed (mg/kg)/[final P solution (mg/L) 0.41 ]. Phosphorus sorption capacity (PSC) was calculated for each soil according to Borling et al. (2001;

Bayesian Statistics
Bayesian inference is the process of fitting a probability model to a set of data and summarising the results via probability distributions on the parameters of the model. In a Bayesian model, the parameters are considered to be uncertain and are therefore represented by probability distributions, known as prior distributions, or priors. These are combined with the model likelihood to obtain posterior distributions that provide their most probable values given the data and modelling assumptions (Berger, 1993). The posterior distributions are summarised using summary statistics. Using a Bayesian framework allows the inclusion of uncertainty in parameters and models, the probabilistic representation of the model outcomes, as well as conveniently constructing very complex models.
Consequently we developed a model with an exponential decay function. The model was based on (1) P concentration previously measured (P init concentration), and (2) the amount of P added as on-going P fert , with P concentration for each plot being repeatedly measured over time. Although P init and P fert were used in constructing data summaries, they were not independent of each other within the model due to the overlapping treatment structure i.e., the 16 treatments comprised the combination of the four P ferts with the four P init concentrations. Time was measured in 12 month intervals for CaCl 2 -P and 6 month intervals for Olsen-P and Colwell-P. Thus, the CaCl 2 -P parameter estimates were based on less than half the data of Olsen-P and Colwell-P. All treatments were measured for Olsen-P and Colwell-P, while CaCl 2 -P was measured from nine selected treatments; the two lowest and the highest P init concentrations receiving either zero or one of the two highest P fert rates. The model was fitted using Markov Chain Monte Carlo techniques by means of JAGS software, version 3.1, using 500 000 iterations, a burn-in of 250 000, and thinning every 250 iterations. The model used Eq. 1 to reflect the changes in soil P that were offset by P fert application. In Eq. 1 Y p,t was the soil P concentration at time (t) for each treatment plot (p); Y p,t−1 represented the previous soil P concentration (P init ); A p,t−1 was the amount of P applied since the previous soil P concentration measurement (divided by 1000 for numerical convenience); S represented the time in years since previous measurement; and e p and d p were parameters which were estimated and were required to be positive. e p represented the increase in extractable soil P for each unit of applied P, and d p represented a decay constant that reflects how quickly the extractable soil P decreases i.e., moves from the labile to unavailable pools, is removed in product, or otherwise lost from the soils.
Final Pconcentration = (previously measured P concentration + e p x P fertiliser applied) exp (−d p x years since previous P concentration measurement) The model is universal for all three P extractions (mg/kg), P fertiliser applied is expressed as kg P/ha, and estimated soil specific e and d values are provided (see Table 4). The observed P final concentrations for each of the measures were denoted by P p,t , which was assumed to be normally distributed, with means given by the modelled value Y p,t , and with reciprocal variance τ (Eq. 3).
Prior distributions were assigned to the parameters of the Bayesian model ( Table 2). Vague priors were used since they make minimal assumptions, resulting in the inferences being TABLE 2 | Assignment of priors to parameters e p and d p for a Bayesian model of soil P concentration change across six soils, each with a range of P init concentrations and P ferts .

P fert constant
Decay constant The observation precision, τ, in Eq. 3 and all precisions in this table (τ es , τ ep , τ e 0 , τ dp , τ ds , τ d 0 ) were assigned vague gamma priors with shape and rate parameters set to 0.01; e.g., τ es ∼ Ga (0.01, 0.01). The initial predicted value; Y p , 1 was unlike later predictions in that it was assumed to not be subject to earlier P applications, but informed only by the first observed P value, and was assigned a uniform distribution so that any value between the upper and lower limits was equally probable. informed by the data rather than the prior parameters. The parameters e p and d p were assumed to be normally distributed, each with soil specific means, e s and d s , where the subscripts p and s are indices for plot and soil. The e s and d s were assumed to be normally distributed with common global means, e 0 and d 0 , respectively. Thus, plot parameters could vary around their soil means, which themselves were distributed around a global mean. The normal distributions were truncated on the left at zero so that the parameters remained positive ( Table 2). Using the model, the time required for P init to decrease to agronomic optimum -defined as 17 mg/kg for Olsen-P (Gourley et al., 2007) and as 55,40,23,29,55, and 29 mg/kg for Soils 1 to 6, respectively, for Colwell-P (Gourley et al., 2007) and an environmental threshold (defined as <0.25 mg/kg for CaCl 2 -P), was calculated. We chose the CaCl 2 -P threshold of 0.25 mg/kg (equivalent to solution concentration of 0.05 mg/L based on a 5:1 solution:soil extraction ratio) because the phosphorus trigger value for slightly disturbed lowland streams in Tasmania and Victoria at the time of the study was 0.05 mg/L (ANZECC, 2000) and CaCl 2 extractable P (mg/L) has been shown to approximate runoff P concentrations in a number of studies (Dougherty et al., 2008;Burkitt et al., 2010).

Predicted Versus Observed Final CaCl 2 -, Olsen-and Colwell-P Concentrations
Posterior predictions and observed soil CaCl 2 -P, Olsen-P and Colwell-P concentrations (Figure 1) were highly correlated (r 2 > 0.85), showing the model successfully fitted the dataset well for all three P extractions based on 4.5 years of data for 6 soil which ranged from extremely low to high P sorption capacity. Final Olsen-P and Colwell-P concentrations, which were initially approximately twice agronomic optimum or lower, were the most accurately predicted. The precision of the model for P concentrations which were initially near agronomic optimum should be reassuring to land managers and policy makers as it limits the chances of soil P concentrations decreasing below agronomic optimum faster than predicted. However, the large range in extractable soil P concentrations of the dataset resulted in modelled values diverging from observed values as soil P concentrations increased (Figure 1). This divergence meant that the model overestimated the P final concentration i.e., the greater the P init concentration, the greater the predicted P final concentration was compared to the actual P final concentration (Figure 1). This may have been the result of the smaller population of data at high P init concentrations, for example, Soil 3 did not contain high Olsen and Colwell P init concentrations.
Since the dataset used to develop this model was sourced from grazed pasture field soils, some decreases in extractable soil P would be the result of P being exported in grazing animal products (estimated to be between 2.9 and 12.3 kg P/ha.yr depending on site), and environmental losses. These estimated exports changed according to soil and P init concentration, as would any environmental losses which may occur. The current model includes these exports as contributing to the decreases in extractable soil P concentrations. Further population of the model with data which compares decreases in extractable soil P and various rates of P export from a soil is required, and may lead to incorporation of a P export term. Nevertheless, the majority of the decrease in extractable soil P may be driven by readily extractable P being sorbed by the soil (Barrow, 1980) and conversion to forms not extractable by CaCl 2 or bicarbonate, for example organic P (Po) or recalcitrant inorganic-P (Pi).
It could be postulated that our model was predicting decreases in the effectiveness of the capital fertiliser applied at trial site establishment. As a result, the P final concentrations we observed may be expected to be lower than those in soils where similarly high P init concentrations were established over a longer period. However, if that were the case, our model would more closely represent the higher P final concentrations since the model overestimated P final , primarily for high P init soils. Irrespective of this, greater decreases in soil extractable P have been reported for high P init concentrations, whether built-up over a short or long period (McCollum, 1991). McCollum (1991) found that the loss constants for Mehlich-1 P init , when the P init was achieved either through a single application of P fertiliser or progressively over 8 years of P fertiliser applications, had similar loss constants of −0.097 (r 2 = 0.89) and −0.073 (r 2 = 0.81), respectively, for initial Mehlich-1 P of 100 g m −3 ; −0.068 (r 2 = 0.73) and −0.080 (r 2 = 0.72) for initial values of 40-60 g m −3 ; and −0.047 (r 2 = 0.52) and −0.048 (r 2 = 0.41) for initial values of 20-30 g m −3 . When averaged across all soils P init concentrations and P fert , our model overestimated final Olsen-P and Colwell-P concentrations by 6.1 mg/kg (32%) and 4.3 mg/kg (10%), respectively ( Table 3). The overestimation for final Colwell-P was within analytical error, suggesting that the model was most successful at predicting final Colwell-P concentrations. Still, the overestimation for Olsen-P was relatively small considering the range in P init concentrations and soils. The Olsen-P method involves a short extraction (0.5 h) and generally extracts soluble orthophosphate P, including P that can be rapidly desorbed from soil particles (Coventry et al., 2001). The longer extraction time (16 h) for the Colwell-P method measures orthophosphate as well as some molybdate-reactive non-orthophosphate P compounds (Coventry et al., 2001). Thus, in addition to the background P exports, the model reflects the conversion of P from readily available to less available forms, and to lesser degree, the rapid and reversible sorption reactions which would have influenced the decreases in CaCl 2 -P and Olsen-P concentrations. Similar to Olsen-P and Colwell-P, the divergence of the predicted CaCl 2 -P concentrations from the observed concentrations was greater as P init concentrations increased (Figure 1). The greater a soil's CaCl 2 -P concentration, the greater the amount of P available to be converted to less available P through processes such as sorption and precipitation, to be removed in animal products, and to be lost to the surrounding environment. The model overestimated average final CaCl 2 -P concentrations by 3.07 mg/kg (26%) ( Table 3). The P final concentrations for Soils 3 and 6, with the highest CaCl 2 -P concentrations and lowest PBIs, were overestimated by a mean of 9.1 mg/kg (140%) ( Table 3). Whilst these may be large overestimations when compared to the environmental CaCl 2 -P threshold i.e., ≤0.25 mg/kg (0.05 mg/L), they were not when compared to the initial CaCl 2 -P concentration range of the dataset ( Table 5). We hypothesise that the extremely high CaCl 2 -P concentrations measured from Soils 3 and 6, and the fewer data available for CaCl 2 -P, contributed to the overestimation. The extreme CaCl 2 -P concentrations of Soils 3 and 6 appear to be due to these soil's low PBI. Soils 3 and 6 contained the lowest oxalate-extractable iron and aluminium contents, lowest clay and highest sand contents, and some of the highest organic carbon contents, of all the soils examined in the study (Table 1).
Phosphorus movement from the surface 10 cm, and conversion to forms not CaCl 2 -P extractable, may have had a greater contribution to P decreases compared to sorption processes for Soils 3 and 6. Indeed, the greater a CaCl 2 -P concentration, the greater the risk of environmental loss . This is supported by the modelled P final concentrations being much closer to the observed P final concentrations for the very low PBI Soil 4 which had lower initial CaCl 2 -P concentrations than the extremely low and very low PBI Soils 3 and 6, respectively ( Table 3). Final CaCl 2 -P concentrations for Soils 1, 2, and 5, which had much lower initial CaCl 2 -P concentrations and higher PBIs than Soils 3, 4, and 6, were underestimated by an average of −0.13 mg/kg (−53%).

Influence of P fert on the Modelled P final Concentration
It is important to examine the influence of P fert on the modelled P final concentration because it is a variable that the land manager can control. Whilst some land managers would rather not omit P fert inputs altogether, soil P concentrations may still decrease if the P fert applied is less than that required to maintain soil P concentrations. There was a range in e s values (increase in extractable soil P concentration for each unit of applied P) across the soils, where smaller e s values indicated soils which were less responsive to P fert when presented in isolation from the soil effect (Table 4). Changes in P concentrations of Soils 3 and 6 were the least sensitive to P fert when extracted with CaCl 2 ( Table 4); contrasting with current assumptions that relatively low rates of P fertiliser are required to increase fertility of low PBI soils (Burkitt et al., 2001). This is likely to largely be due to the on-going P fert moving beyond the surface 10 cm of these very low PBI, coarse textured soils which receive high annual rainfall (∼1100 mm). Indeed, difficulty was encountered in increasing agronomic extractable soil P concentrations for the extremely low PBI Soil 3. Despite application of 500 kg P/ha, split over 12 months, the Olsen-P of Soil 3 only increased from 13 to 19 mg/kg. Alternatively, the high initial CaCl 2 -P concentrations of these soils, and conversion of the applied P fert to forms unextractable by CaCl 2 , such as Po through incorporation into the plant biomass and microbial population, may have also contributed to this low sensitivity to the applied P fert . Olsen-P and Colwell-P concentrations of Soils 3 and 6 were more sensitive to P fert than Soils 1 and 2 which had higher PBI (Table 4). This suggests that Soils 3 and 6 had some capacity to weakly sorb P which was not CaCl 2 -P extractable, and/or desorb molybdatereactive non-orthophosphate P compounds.

Influence of Soil on the Modelled P final Concentration
Despite soil PBI being an important factor in determining the appropriate P fertiliser rate to increase soil P concentrations, few studies have modelled P decreases across a range of soils. e represents the increase in extractable soil P for each unit of applied P, and d represents a decay constant that reflects how quickly the extractable soil P decreases. The larger the number, the more sensitive the extraction is to changes in soil P. Kerridge et al. (1990) reported that the negligible P sorption capacity of the soils they studied meant that a soil term did not need to be included in the model they used. However, the dataset that was used to construct the current model, showed that extractable P concentration decreases did vary according to soil. Olsen-P and Colwell-P concentrations decreased fastest (highest d s values) for Soils 3 and 6 (Table 4 and Figure 2), despite having some of the lowest P init concentrations. It is thought that vertical movement of P beyond the soil sampling depth, and conversion to forms not directly extractable, were the two most likely pathways responsible for the decreases in extractable P for Soils 3 and 6. However, despite absolute decreases in CaCl 2 -P concentrations of Soils 3 and 6 being the largest, they had the lowest ds, suggesting the CaCl 2 -P concentrations of the two lowest PBI soils decreased more gradually than those of higher PBI soils ( Table 4, and Figure 2). Since the CaCl 2 -P concentrations of Soils 3 and 6 were similar to the Olsen-P concentrations, it is possible that their high organic carbon content may have led to greater hydrolysis of Po to readily available P. The acidic molybdate-solution may cause hydrolysis of Po so that the molybdate-blue colour reaction of Murphy and Riley (1962) does not exclusively represent Pi (Haygarth and Sharpley, 2000). Thus, it is possible that the gradual decreases in CaCl 2 -P concentrations were partly due to decreases in Po.
We have successfully modelled decreases across a range of soil types, and identified varying decreases according to P extraction, soil type, and PBI. However, we acknowledge that these findings are based on a limited data set of 6 soils sampled over a period of 4.5 years and we suggest more data on P decreases across a wider range of soil types are needed before a term based on PBI can be incorporated into our model. If such a soil term could be incorporated, we suggest the accuracy of the predictions would increase further.

Influence of P fert and Soil on Modelled P Decrease
The magnitude of the effect of P fert was less than soil since the decreases in P were more variable between soils than between rates of P fert ( Table 4). This may in part be due to different animal product exports (2.9 to 12.3 kg P/ha.yr according to soil and soil P concentration), and potential for surface runoff and leaching of P according to soil, from the trial sites which populated our model. Schulte et al. (2010) reported that their large P deficit due to P export in fodder removed from their trial sites was the driver of P decreases (contributing 63%) for eight moderate to low PBI soils. Despite this, and the small range in PBI of the soils studied by Schulte et al. (2010) they reported a relationship between Morgan-P decreases and parent material of the individual soils that was marginally insignificant (P = 0.052). The greater influence of soil than P fert in the current study was thought to be due to the large range of soils examined and supports the need to investigate P decreases across a more diverse set of soils which cover a wider range of P sorption and their properties.
Since there was one experimental site per soil, each of the site-level parameters corresponds to a soil type with associated properties such as PBI. However, a highly significant (P < 0.001; r 2 = 0.95) linear relationship existed between Ps 1000 and the d p parameter for CaCl 2 -P. Due to the PBI scale not being linear; i.e., a PBI of 400 does not have double the sorption capacity of a PBI of 200, there was a highly significant (P < 0.001) logarithmic relationship between PBI and the d p parameter for CaCl 2 -P (r 2 = 0.97). There was also a significant (P < 0.05) linear relationship between Ps 1000 and d p for Olsen-P and Colwell-P (r 2 = 0.24 and 0.38, respectively). However, the polynomial relationship for d p of Olsen-P was stronger (r 2 = 0.40) than the linear relationship, as was the exponential relationship for Colwell-P d p (r 2 = 0.57). Although the parameters for this model are soil specific, information on their magnitude and interrelationship will assist in estimating decreases for other soils. For example, changes in P concentrations of low PBI soils which are coarse textured, contain high organic matter contents, and have low clay and oxalate-aluminium and oxalate-iron contents, would be expected to be different to medium and high PBI soils.

Time Until Environmental and Agronomic Optimum P Concentrations Are Reached
Based on the current data set, the model estimated that for the majority of the soils, it would take up to 2.7 years for CaCl 2 -P concentrations of ≤1.71 mg/kg (CaCl 2 -P concentrations of Soils 1, 2, 4 and 5) to decrease to an environmentally acceptable level of ≤0.25 mg/kg (0.05 mg/L) ( Table 5). This is a shorter period of time than that reported by Dodd et al. (2012) of 22-44 years for surface 7.5 cm water extractable-P (WEP) concentrations of 0.039-0.141 mg/L to decrease to an environmentally acceptable level of 0.02 mg/L (0.10 mg/kg). When Dodd et al. (2012) reported WEP and CaCl 2 -P concentrations were plotted, there was no relationship between the two measures. This explains why Dodd et al. (2012) did not report exponential decreases in CaCl 2 -P for three of their four sites, while they did report exponential decreases in WEP for three sites. Dodd et al. (2012) estimated that WEP concentrations of the two soils with the greatest P retentions (PR) would take twice as long to reach the target concentration than a soil with a lower PR. We converted the PR (%) reported by Dodd et al. (2012) to approximate PBI values using data from Burkitt et al. (2002) to conclude that the PBI range was approximately 90-600 for their lowest P fertiliser treated soils. Based on this conversion, the PBI's of two of Dodd et al. (2012) soils were less than 150, as were three of our soils. When we modelled our two lowest PBI soils the model estimated it would take on average 20 years, ranging from 3 to 35 years, for CaCl 2 -P concentrations to decrease from between 19 and 28 mg/kg for Soil 3, and on average 34 years, ranging from 17 to 41 years, to decrease from between 10 and 29 mg/kg for Soil 6; to ≤0.25 mg/kg ( Table 5). For our CaCl 2 -P concentrations to decrease to Dodd et al. (2012) WEP threshold of 0.10 mg/kg, our model predicted it would take 37 years, ranging from 4 to 53 years for Soil 3, and 47 years, ranging from 21 to 51 years, for Soil 6. While the latter mean estimates were similar to those reported by Dodd et al. (2012) for WEP, in contrast to Dodd et al. (2012) our lowest PBI soils took longer to decrease than higher PBI soils. This may be due to our initial CaCl 2 -P concentrations, which were greatest for the lowest PBI soils, being much greater than the respective CaCl 2 -P and WEP concentrations of Dodd et al. (2012). This suggests that the extreme CaCl 2 -P concentrations of our two lowest PBI soils, and high CaCl 2 -P concentrations of our high PBI soils, decreased more rapidly than those reported by Dodd et al. (2012) supporting the exponential decrease we modelled for CaCl 2 -P.
Our findings are in line with simple model calculations by Schoumans and Groenendijk (2000) which showed that it would take many decades before high soil P concentrations in noncalcareous low P sorbing sandy soils would be reduced to environmentally acceptable concentrations. It is assumed that Po at least partly buffers the decreasing extractable P concentrations. Considering Pi is the dominant P form in highly P fertilised soils (Dougherty et al., 2006) that are able to hold P, it is assumed that in low PBI soils which only hold small amounts of Pi, Po supply via microbial processes can buffer loss of Pi that occur via desorption.
An important finding of the model is that even if some soils are at or close to agronomic optimum, environmental losses can still occur. For the extremely low PBI Soil 3 and very low PBI Soil 6, CaCl 2 -P concentrations were above the environmental threshold when Olsen-P and Colwell-P concentrations were at or just above agronomic optima. These findings suggest the use of such soils for agricultural purposes could pose a long-term elevated risk to the environment as optimum agronomic production is unlikely to be achieved without a high risk of P loss. Similarly, for the low PBI Soil 4, Olsen-P concentrations below agronomic optima resulted in CaCl 2 -P concentrations above the environmental threshold. Since the CaCl 2 -P concentrations of Soil 4 were not as great as those of Soils 3 and 6, the time taken to reach the environmental threshold was less than 6 years. However, not all soils will pose an environmental risk when soil P concentrations are at agronomic optimum. From Table 5, CaCl 2 -P concentrations of Soils 1 and 2 were not above the environmental threshold until Olsen-P concentrations were twice, and Colwell-P concentrations three to four times the agronomic optima. For Soil 5, Olsen-P and Colwell-P concentrations at or near the agronomic optima were associated with CaCl 2 -P concentrations close to the environmental threshold.
The current model predicts it would take, on average, 14 years, ranging from 11 to 20 years for Olsen-P concentrations of between 34 and 44 mg/kg to decrease to agronomic optimum of 17 mg/kg. An initial Olsen-P concentration of between 55 and 96 mg/kg would take on average 32 years, ranging from 26 to 49 years, depending on the soil ( Table 5). In comparison, Dodd et al. (2012) reported a shorter period of up to 7 years for initial Olsen-P concentrations (0-7.5 cm) to decrease from 51 to 25 mg/kg. This is equivalent to an initial Olsen-P of ∼43 mg/kg decreasing to 22 mg/kg if sampled to 10 cm (Coad et al., 2010) which our model suggests would take 25 years when using global d and e values. The time reported by Dodd et al. (2012) for Environmental threshold for CaCl 2 -P was defined as <0.25 mg/kg, while agronomic optima for Olsen-P was defined as 17 mg/kg, and Colwell-P as 55, 40, 23, 29, 55, and 29 mg/kg for Soils 1 to 6, respectively. n/a means not applicable since P init was already below agronomic optimum, or was not analysed for CaCl 2 -P.
Olsen-P concentrations to decrease to agronomic optimum was shorter than the time they reported for WEP concentrations to decrease to the environmental threshold, despite WEP decreasing at a greater rate. This suggests that the majority of the soils they studied may have soil properties similar to Soils 3, 4 and 6 of the current study, in that they cannot be maintained at agronomic optimum without environmental loss. Colwell-P agronomic optimum varies according to soil and as such, it would take on average 9 (mean range of 7 to 11 years), 47 (mean range of 39 to 59 years), and 84 years (mean range of 69 to 117 years) for initial Colwell-P concentrations of 30-90, 102-181, and 274-461 mg/kg to decrease to their respective agronomic optima for Soils 1 to 6 ( Table 5).
From Table 5, and as mentioned above, there is a large range in predicted decreases for our model. Schulte et al. (2010) reported an uncertainty in their model predictions of 3 to >20 years for Morgan-P to decrease from Index 4 to 3. Dodd et al. (2012) reported an estimate of 23 years could be 21 or 31 years for WEP. Our greater range in predicted decrease period may be due to the greater range of soils used. In addition, our model extrapolated our dataset more than Dodd et al. (2012) who had field data for 7, 16, 21, and 26 years, depending on site. However, the advantage of our model, as opposed to a purely descriptive empirical one, is that the predictions beyond the observed data range are reasonable, as long as the model assumptions are valid. We suggest that the large range in P init concentrations of the dataset used for the current model contributed to these wide ranges, since some of our low P init concentrated soils had low prediction ranges, similar to those reported by Dodd et al. (2012), with the range widening as P init concentrations increased. In addition, the on-going P fert would have extended some of our periods of decrease.

CONCLUSION
The model developed in this study has provided critical information for our long term management of soil P concentrations, as we found that decreases in soil P are generally slow (14 years (ranging from 11 to 20 years) for Olsen-P concentrations of between 34 and 44 mg/kg to decrease to an agronomic optimum of 17 mg/kg). As a result of our findings, we suggest that policy makers need realistic expectations about achievable P decreases, and land managers can be reassured of the implications of withholding P fertiliser, as soil P concentrations are unlikely to decline rapidly on these Australian pasture soils. Our findings also highlight the risk of using agronomic soil test concentration as a reliable guide to environmental risk, as some soils in the study did not exceed the environmental threshold until Olsen-P concentrations were twice, and Colwell-P concentrations three to four times the agronomic optima, whereas soils with extremely low and very low PBI values exceeded the threshold, despite Olsen-P and Colwell-P concentrations being at or just above agronomic optima.
To achieve more rapid decreases in soil P concentrations, increases in P exports in plant and animal products would need to occur in addition to withholding P fertiliser application. To enable modelling of scenarios where various P exports occur, a P export term needs to be incorporated into future versions of the model. Within the 6 soils examined in the current study, results have shown that soil with inherently very low and extremely low PBIs, require different management than higher PBI soils, in order to limit P losses i.e., application of different quantities and forms of P. However, further work on a greater range of soils types is needed to understand the influence of soil properties such as PBI on extractable soil P decreases, in order to refine this model.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
JT completed all field sampling, laboratory and non-Bayesian data analysis, and presentation, and drafted the manuscript. RC undertook all the Bayesian modelling and assisted JT with the interpretation of these data. LB planned the field experiment and made a substantial intellectual contribution to the research, and the interpretation of data presented in this manuscript. WD made a substantial intellectual contribution in terms of the interpretation and presentation of results presented in this manuscript. All authors contributed to the article and approved the submitted version.