Organic phosphorus cycling may control grassland 1 responses to nitrogen deposition : a long-term field 2 manipulation and modelling study 3

Abstract. Ecosystems limited in phosphorous (P) are widespread, yet there is limited
understanding of how these ecosystems may respond to anthropogenic
deposition of nitrogen (N) and the interconnected effects on the
biogeochemical cycling of carbon (C), N, and P. Here, we investigate the
consequences of enhanced N addition for the C–N–P pools of two P-limited
grasslands, one acidic and one limestone, occurring on contrasting soils, and we
explore their responses to a long-term nutrient-manipulation experiment. We
do this by combining data with an integrated C–N–P cycling model (N14CP). We
explore the role of P-access mechanisms by allowing these to vary in the
modelling framework and comparing model plant–soil C–N–P outputs to
empirical data. Combinations of organic P access and inorganic P
availability most closely representing empirical data were used to simulate
the grasslands and quantify their temporal response to nutrient
manipulation. The model suggested that access to organic P is a key
determinant of grassland nutrient limitation and responses to experimental N
and P manipulation. A high rate of organic P access allowed the acidic
grassland to overcome N-induced P limitation, increasing biomass C input to
soil and promoting soil
organic carbon (SOC) sequestration in response to N addition. Conversely,
poor accessibility of organic P for the limestone grassland meant N
provision exacerbated P limitation and reduced biomass input to the soil,
reducing soil carbon storage. Plant acquisition of organic P may therefore
play an important role in reducing P limitation and determining responses
to anthropogenic changes in nutrient availability. We conclude that
grasslands differing in their access to organic P may respond to N
deposition in contrasting ways, and where access is limited, soil organic
carbon stocks could decline.



Introduction 51
widespread, than N limitation [Fay et al., 2015;Du et al., 2020;Hou et al., 2020]. In a meta-analysis of 75 grassland nutrient addition experiments spanning five continents, Fay et al. [2015] found that 76 aboveground annual net primary productivity was limited by nutrients in 31 out of 42 sites, most 77 commonly through co-limitation of N and P [Fay et al., 2015]. Similarly, P additions in 652 field 78 experiments increased aboveground plant productivity by an average of 34.9% [Hou et al., 2020], 79 while it is estimated that P limitation, alone or through co-limitation with N, could constrain up to 80 82% of the natural terrestrial surface's productivity [Du et al., 2020]. 81 Furthermore, P limitation may be exacerbated by N deposition [Johnson et al., 1999;Phoenix et al., 82 2004], or become increasingly prevalent as previously N-limited ecosystems transition to N-sufficient 83 states [Goll et al., 2012]. For example, in parts of the Peak District National Park, UK, N deposition has 84 exceeded 3 g m -2 yr -1 , with further experimental additions of 3.5 g m -2 yr -1 leading to decreases rather 85 than increases in productivity of limestone grasslands [Carroll et al., 2003]. This makes P limitation 86 critical to understand in the context of global carbon and nutrient cycles. By definition, N deposition 87 should impact P-limited ecosystems differently to N-limited ones, yet there is little understanding of 88 how N deposition impacts these systems. 89 While N deposition may worsen P limitation in some instances, plant strategies for P acquisition may 90 require substantial investments of N, suggesting that increased N supply may facilitate enhanced P 91 uptake [Vance et al., 2003;Long et al., 2016;Chen et al., 2020]. Indeed, previous work from long-term 92 experimental grasslands has shown strong effects of N deposition on plant enzyme production 93 [Johnson et al. 1999;Phoenix et al. 2004], whereby the production of additional extracellular 94 phosphatase enzymes was stimulated. While it is not clear if this response is driven by exacerbated P-95 limitation resulting from N deposition or extra N availability making elevated enzyme production 96 possible, such changes in plant physiology may promote cleaving of P from organic soil pools. Over 97 time, the accumulation of plant-available P from organic sources may provide a mechanism by which 98 plants exposed to high levels of N deposition may overcome P limitation [Chen et al. 2020]. 99 By using the integrated C-N-P cycle model N14CP, Janes-Bassett et al. [2020] suggest that the role of 100 organic P cycling in models may be poorly represented, as the model failed to simulate empirical yield 101 data in agricultural soils with low P fertiliser input. Organic P access is therefore likely an important 102 means of nutrient acquisition for plants in high N and low P soils [Chen et al. 2020], yet our 103 understanding of organic P cycling in semi-natural ecosystems is fairly limited [Janes-Bassett et al. 104 2020]. Such interdependencies of the C, N and P cycles make understanding an ecosystem's response 105 to perturbations in any one nutrient cycle challenging, particularly when ecosystems are not solely 106 limited in N. This highlights the need for integrated understanding of plant-soil nutrient cycling across 107 the C, N and P cycles, and in ecosystems that are not solely N-limited. 108 Process-based models have a role to play in addressing this, as they allow us to test our mechanistic 109 understanding and decouple the effects of multiple drivers. There has been increasing interest in 110 linking C with N and P cycles in terrestrial ecosystem models [Wang et al., 2010;Achat et al., 2016;111 Jiang et al., 2019] as the magnitude of the effects that anthropogenic nutrient change can have on 112 biogeochemical cycling are realised [Yuan et al., 2018]. Yet, few modelling studies have explicitly 113 examined the effects of P limitation, or the role of organic P access in determining nutrient limitation, 114 likely mirroring the relatively fewer empirical studies of these systems. 115 By combining process-based models with empirical data from long-term nutrient-manipulation 116 experiments, we may simultaneously improve our understanding of empirical nutrient limitation, the 117 role(s) of organic P acquisition, and their interactions with anthropogenic nutrient pollution. In 118 particular, this approach offers a valuable opportunity for understanding ecosystem responses to 119 environmental changes that may only manifest after extended periods of time, such as with changes 120 in soil organic C, N and P pools, which typically occur on decadal timescales [Davies et al., 2016a, Janes-121 Bassett et al., 2020]. Here, we combine new data from a long-term nutrient manipulation experiment 122 on two P-limited upland grasslands (acidic and limestone) occurring on contrasting soils, with the 123 mechanistic C-N-P plant-soil biogeochemical model; N14CP [Davies et al., 2016b]. 124 We use these experimental data to explore the role of organic P access in determining ecosystem 125 nutrient limitation and grassland responses to long-term nutrient manipulations. Specifically, we aim 126 to explore how variation in P acquisition parameters, that control access to organic and inorganic 127 sources of P in the model, may help account for differing responses of empirical grassland C, N and P 128 pools to N and P additions. Second, we explore the effects of long-term anthropogenic N deposition 129 and experimental N and P additions on plant and soil variables of the simulated acidic and limestone 130 grasslands. This will help improve our understanding of organic P process attribution within the model 131 and may suggest how similarly nutrient limited grasslands could respond to similar conditions. 132 We hypothesise that 1) access to organic P will be an important determinant of ecosystem nutrient 133 limitation, 2) increased organic P availability may alleviate P limitation resulting from N deposition and 134 3) grasslands capable of accessing sufficient P from organic forms may overcome P limitation resulting 135 from N deposition and nutrient treatments, whereas grasslands lacking such accessibility will not. 136 137 138 139 2. Methods 140

Field experiment description 141
The empirical data is from Wardlow Hay Cop (henceforth referred to as Wardlow), a long-term 142 experimental grassland site in the Peak District National Park (UK) [Morecroft et al., 1994]. Details of 143 empirical data collection are available in supplementary section 1. There are two distinct grassland 144 communities occurring in close proximity; acidic (National vegetation classification U4e) and 145 limestone (NVC CG2d) semi-natural grasslands (Table S2). Both grasslands share a carboniferous 146 limestone hill but the limestone grassland sits atop a thin humic ranker [Horswill et al., 2008] and 147 occurs predominantly on the hill brow. In contrast, the acidic grassland occurs in the trough of the 148 hill, allowing the accumulation of wind-blown loess and the formation of a deeper soil profile of a 149 palaeo-argillic brown earth [Horswill et al., 2008]. 150 Despite contrasting soil types, both the acidic and limestone grasslands are largely P-limited 151 [Morecroft et al. 1994;Carroll et al. 2003], though occasional N and P co-limitation can occur 152 [Phoenix et al. 2003] and more recently, positive growth responses in solely N-treated plots have 153 been observed, in line with the latest understanding that long-term N loading may increase P supply 154 by increasing phosphatase enzyme activity [Johnson et al. 1999 (distilled water control -0N), 3.5 (low nitrogen -LN) and 14 g N m -2 yr -1 (high nitrogen -HN). The P 163 treatment is applied at a rate of 3.5 g P m -2 yr -1 (phosphorus -P). 164 Data collected from the Wardlow grasslands for the purpose of this work are; aboveground biomass 165 C, SOC, and total N, which is assumed to be equivalent to modelled SON. This new data is combined 166 with total P data that was collected by Horswill et al. at the site [Horswill et al., 2008]. Summaries of 167 these data are available within the supplementary material (Table S4)  The N14CP ecosystem model is an integrated C-N-P biogeochemical cycle model that simulates net 173 primary productivity (NPP), C, N and P flows and stocks between and within plant biomass and soils, 174 and their associated fluxes to the atmosphere and leachates [Davies et al., 2016b]. N14CP was 175 originally developed and tested on 88 northern Europe plot-scale studies, including grasslands, 176 where C, N and P data were available. All but one of the tested ecosystems exhibited N limitation 177 [Davies et al., 2016b]. It has also been extensively and successfully blind-tested against SOC [Tipping 178 et al., 2017] and NPP data from unimproved grassland sites across the UK [Tipping et al., 2019]. 179 However, N14CP has not been extensively tested against sites known to exhibit P limitation, 180 especially where these are explicitly manipulated by long term experimental treatments. While the 181 importance of modelled weatherable P (P Weath0 ) and historic N deposition on N-limited C, N and P 182 have been investigated [Davies et al. 2016b], the potential influence of organic P on ecosystem 183 nutrient limitation and responses to nutrient perturbations have yet to be explored. 184 Here, we modify N14CP to add experimental N and P additions to simulate a long-term nutrient 185 manipulation experiment similar to that at the limestone and acidic grasslands at Wardlow, and we 186 use empirical data from Wardlow to explore the role of organic P cleaving in determining ecosystem 187 state. A full model description can be found in Davies et al., [2016b], however, a summary of the 188 most relevant features is given here for convenience. 189

Net primary productivity and nutrient limitations 190
Plant biomass is simulated in the model as two sets of pools of coarse and fine tissues representing 191 both above and belowground plant C, N and P, with belowground biomass for each plant functional 192 type (PFT) represented by a root fraction. NPP adds to these on a quarterly basis with growth 193 occurring in quarters 2 and 3 (spring and summer). In N14CP, NPP depends on a single limiting 194 factor, in accordance with Liebig's law of the minimum. The factors that can limit growth in the 195 model include available N and P, temperature or precipitation, the latter two being provided as 196 input driver data (see section 2.3.2). 197 First, the potential maximum NPP limited by climate is calculated using regression techniques, as in 198 Tipping et al. [2014]. The corresponding plant demand for N and P to achieve this potential NPP is 199 values for a given plant functional type and is downregulated by anthropogenic N deposition, but 213 not soil N content more generally, as it is assumed that atmospherically deposited N is readily 214 available to N-fixers. Nitrogen fixation in the model is also related to P availability. The degree to 215 which P availabilty limits this maximum rate of fixation is determined by a constant; K Nfix [Davies et 216 al. 2016b]. The initial rate of N fixation is based on literature values for a given PFT but is 217 downregulated by anthropogenic N deposition and related to P availability. The degree to which P 218 availabilty limits this maximum rate of fixation is determined by a constant; K Nfix [Davies et al. 219 2016b]. This means that while modelled NPP is limited by availability of a single nutrient, co-220 limitation may occur through P limitation of N fixation [Danger et al. 2008]. to its rapid rate of turnover whereas coarse litter contributes N and P through decomposition and 230 does not join the SOM pool. Plant available P also comes from SOM and coarse litter decomposition, 231 direct treatment, desorption of inorganic P from soil surfaces, and sometimes cleaving of organic P 232 [Davies et al., 2016b]. The sorbed inorganic P pool builds over time with inputs of weathered P and 233 sorption of any excess plant available inorganic P, and desorption occurs as a first order process. 234 Phosphorus enters the plant-soil system by weathering of parent material, the initial value of which 235 (P Weath0 within the model) can be set to a default value, or made site-specific by calibrating this initial 236 condition to soil observational data (as in methods section 2.3.3). From this initial pool, annual 237 releases of weathered P are determined by first-order rate constants that are temperature 238 dependent, with the assumption that no weathering occurs below 0 degrees Celsius. This weathered 239 P can then contribute toward plant-available P in soil water or be sorbed to soil surfaces. In principle, 240 P can be added in small quantities by atmospheric deposition [Ridame and Guieu, 2002] or by local 241 redistribution [Tipping et al., 2014]. In principle, P can be added in small quantities by atmospheric 242 deposition [Ridame and Guieu, 2002] but for the purpose of this work, P deposition is set to zero in 243 the model. While the contribution of P through atmospheric deposition is increasingly realised 244 [Aciego et al. 2017], we cannot account for the losses of P that may also occur through landscape 245 redistribution [Tipping et al. 2014].For the purpose of this study, P deposition is set to zero as its net 246 contribution to the total P pool in comparison to weathering is assumed to be minimal. 247 The size of the available P pool is determined by summing: P retained within plant biomass prior to 248 litterfall, inorganic P from decomposition, dissolved organic P and P cleaved from SOP by plants. 249 Accessibility of each P form is determined by a hierarchal relationship in the order mentioned above, 250 whereby plants and microbes access the most readily available P sources first and only move onto 251 the next once it has been exhausted. 252 When N is in sufficient supply and more bioavailable P forms have been exhausted from the total 253 available pool, simulated plants can access P from SOM via an implicit representation of extracellular 254 P-cleaving enzymes with a parameter termed P Cleave . While empirical data quantifying this parameter 255 is scarce, N14CP constrains P Cleave by utilising a maximum SOM C:P ratio; [C:P] fixlim , that ensures SOM 256 stoichiometry is not unrealistically disrupted by excessive removal of organic P (Equation 1). 257 The functioning of the P Cleave parameter, including its stoichiometric constraint, remains the same in 261 this work but we have introduced a modifier to adjust the rate at which plants can access this P 262 Equation 1 source. This parameter; P CleaveMax , represents the maximum amount (g m -2 season -1 ) of cleaved P that 263 plants can acquire from the available P pool to satiate P demand. 264 A fraction of plant biomass is converted to litter in each quarterly time step and contributes a 265 proportion of its C, N and P content to SOM, which is sectioned intro three pools (fast, slow and 266 passive) depending on turnover rate [Davies et al., 2016b]. Soil organic P (SOP) is simulated 267 alongside SOC and SON using C:N:P stoichiometries of coarse and fine plant biomass. Decomposition 268 of SOP, and its contribution to the available P pool, is subject to the same turnover rate constants as 269 for SOC and SON. 270 Carbon is lost as CO 2 following temperature-dependent decomposition and as dissolved organic . Red lines highlight modifications to N14CP for the purpose of this work, including adding experimental nutrients and allowing uptake of cleaved P to be more flexible. Solid lines indicate input to another pool and a dashed line indicates either a feedback or interaction with another pool. In the model, N can enter the available pool via atmospheric deposition, nutrient treatments, biological fixation, and decomposition of coarse litter and SOM. For P, the two main sources are the inorganic sorbed pool and from the turnover of SOM. The former is derived initially from the weatherable supply of P, defined by its initial condition (P Weath0 ). P can also be added to this pool experimentally as with N. The dashed line going from available N and P to N fixation represents the downregulation of N fixation by N deposition and the dependency of N fixation on P availability. The cleaving of organic P from SOM and its incorporation into the plant-available nutrient pool, is represented by the dashed red line and its uptake by plants, determined by P CleaveMax, shown with a solid red line.

Simulating the field manipulation experiment with the model 288
We use data from the Wardlow limestone and acidic grasslands to explore the potential role organic 289 P access may have in determining grassland nutrient limitation when exposed to long-term N 290 deposition and more recently, experimental nutrient manipulation. We use environmental input 291 data collated from Wardlow to drive model processes. Empirical data regarding contemporary soil C, 292 N and P for the contrasting grasslands is used to calibrate the initial size of the weatherable P pool 293 within the model, and to allow access to organic cleaved P to vary to account for patterns in the 294 data. We do not aim to perfectly replicate the Wardlow grasslands but rather use the unique 295 opportunity that Wardlow provides to test our understanding of such P-limited ecosystems and how 296 our conceptualisation of P access mechanisms within the model may affect them. In addition, we 297 can use the model-simulated grasslands to investigate the potential effects of long-term N 298 deposition and nutrient manipulation on ecosystems which may differ in their relative availability of 299 different P forms. 300 (10,000 BP in the model). This is to capture the length of time required for soil formation following 312 deglaciation in north west Europe and is not an attempt to truly model this long term period. 313 Instead, it allows us to form initial conditions for modern day simulations that takes in what we 314 know about the site's history and forcings. 315 To use this spin up phase and simulate contemporary soil C, N and P stocks, we use a variety of input 316 driver data. Inputs nearer the present are more accurately defined based on site-scale 317 measurements and assumptions are made regarding past conditions. This approach of spinning up 318 to present-day observations avoids the assumption that ecosystems are in a state of equilibrium, 319 which is likely inaccurate for ecosystems exposed to long-term anthropogenic changes in C, N and P 320 availability. Input driver data include PFT plant functional type history, climatic data and N 321 deposition data. A summary of the data used for model input is provided in supplementary Table S3. 322 To simulate the sites' PFT plant functional type history, we used data on Holocene pollen 323 stratigraphy of the White Peak region of Derbyshire [Taylor et al. 1994], which captures important 324 information regarding Wardlow's land-use history for the entire duration of the model spin up 325 phase. 326 Input drivers are provided as annual time series to drive the model and as the acidic and limestone 327 sites are co-located, these input timeseries are shared for both grasslands. It is assumed in the 328 model that anthropogenic N deposition was negligible prior to 1800 and the onset of the industrial 329 revolution. After 1800, N deposition is assumed to have increased similarly across Europe [Schopp et 330 al., 2003]. In N14CP, this trend is linearly extrapolated from the first year of data (1880) back to 1800 331 [Tipping et al., 2012]. Data regarding N deposition that is specific to Wardlow was incorporated To provide climate forcing data, daily minimum, mean and maximum temperature and mean 335 precipitation records beginning in 1960 were extracted from the UKPC09 Met office CEDA database 336 (Table S3). The data nearest to Wardlow was calculated by triangulating latitude and longitude data 337 and using Pythagoras' theorem to determine the shortest distance. These data were converted into 338 mean quarterly temperature and precipitation. Prior to this, temperature was assumed to follow 339 trends described in Davies et al. [2016b] and mean quarterly precipitation was derived from Met 340 Office rainfall data between 1960 to 2016 and held constant. 341 342 2.3.3. Model parameters for the acidic and limestone grasslands 343 The N14CP model has been previously calibrated and tested against a wide range of site data to 344 provide a general parameter set that is applicable to temperate semi-natural ecosystems, without 345 extensive site-specific calibration [Davies et al., 2016b]. The majority of those parameters are used 346 here for both grasslands. However, two parameters relating to P sources and processes were 347 allowed to vary between the sites: the initial condition for the weatherable P pool, P Weath0 ; and the 348 rate of plant access to organic P sources, P CleaveMax (Figure 1). We allowed P Weath0 to vary for each 349 grassland as variation in a number of factors including lithology and topography mean that we 350 should expect the flux of weathered P entering the plant-soil system to vary on a site-by-site basis 351 [Davies et al. 2016b]. Indeed, we should expect that P Weath0 differs between the acid and limestone 352 grasslands, as despite their proximity, they have differing lithology. Davies et al. [2016b], show that 353 variation in this initial condition considerably helps explain variance in contemporary SOC, SON and 354 SOP stocks between sites. However, it is difficult to set this parameter directly using empirical data, 355 as information on lithology and P release is limited at the site scale. 356 As this is the first time that N14CP has been knowingly applied to ecosystems of a largely P-limited 357 nature, we also allowed the maximum rate at which plants could access cleaved P (P CleaveMax ) to vary, 358 to investigate how plant P acquisition might change when more readily accessible P forms become scarcer. Empirical quantification of organic P access is poor [Janes-Bassett et al. 2020], hence we use 360 a similar data-driven calibration for P CleaveMax as we do for P Weath0 . 361 We ran a series of simulations systematically varying P Weath0 and P CleaveMax and comparing the results 362 to observations., Wwe simulated the two grasslands and their treatment blocks with a set of 200 363 parameter combinations. This captured all combinations of 20 values of P Weath0 between 50 and 1000 364 g m -2 and 10 values of P CleaveMax between 0 to 1 g m -2 per growing season using a log 10 spacing to focus 365 on the lower range of P CleaveMax values. The P Weath0 range was set to capture the lower end of P Weath0 366 estimates described in Davies et al. [2016b], which were more likely to be appropriate for these P-367 poor sites. We explored a range of values for P CleaveMax , from zero where no access to organic sources 368 is allowed, to 1 g m -2 per growing season -a rate in the order of magnitude of a fertilizer application. 369 The model outputs were compared to measured, SOC, SON and total P (Table S4)  Plant biomass C data were excluded from the cost function to allow for blind testing of the model's 382 performance against empirical observations. As the variable most responsive to nutrient additions, 383 both in terms of rapidity and magnitude of the response, we deemed these the most rigorous data 384 to use for separate testing. We included soil C, N and P data from all nutrient treatments rather than 385 just the control to ensure that the selected parameter combination could better account for 386 patterns in empirical data. For instance, we know that empirical N treatments can increase plant and 387 soil enzyme activity in both Wardlow grasslands, [Johnson et al. 1999 Below, we first present data regarding the results of the calibration of P Weath0 and P CleaveMax for each 402 grassland, and how simulated grassland C, N and P using these parameter combinations compares to 403 the empirical data (section 3.1, Figure 2). Raw empirical data is available in table S1 in section 2 of 404 the supplementary material. Second, we explore how the limiting nutrient of the modelled 405 grasslands has changed through time in response to N deposition and experimental treatment 406 (section 3.2, Figure 3). Third, we explore how C, N and P pools in the simulated grasslands have 407 responded to N deposition and nutrient treatment within the model, and include empirical data to 408 contextualise changes (section 3.3, Figure 4). Finally, we present the C, N and P budgets for both 409 modelled grasslands to examine changes in C, N and P pools more closely, in order to better our 410 mechanistic understanding of changes in nutrient flows within the model (section 3.3, Figure 5). 411 The model calibration selected parameter values for P Weath0 and P CleaveMax that indicate contrasting 415 use of P sources by the two simulated grasslands, with the acidic grassland capable of acquiring 416 more P from organic sources, having a P CleaveMax value of 0.32162 g m -2 season -1 compared to the 417 limestone, with a value 10 times smaller at 0.0316 g m -2 season -1 . Conversely, inorganic P availability 418 was greater in the limestone grassland due to the larger weatherable pool of P, P Weath0 , at 300 g m -2 419 compared to 150 g m -2 in the acidic. 420 The selected parameter combinations resulted in the model simulating the acidic grassland as N-421 limited and the limestone as P-limited, with reasonable congruence between observed and 422 modelled data. The outputs for the calibrated model are shown in Figure 2 against the observations 423 for above-ground biomass C, soil organic C, and N for both the acidic and limestone grasslands (Fig  424   2). Raw data used for Figure 2 are provided in supplementary tables S5 and S6. 425 Overall, N14CP more accurately simulated the magnitude of limestone grassland C, N and P pools 426 than the acidic, and it generally captured the pattern of responses to nutrient treatment, albeit this 427 is not always supported by high r 2 values. The model estimates of above ground biomass C are 428 broadly aligned with the observations: capturing variation between the grasslands and treatments 429 (r 2 = 0.58), and on average overestimating the magnitude by 12.9% (SE ± 11.9) and 12.1% (SE ± 9.4) 430 for the acidic and limestone grasslands respectively (Fig 2a). 431 Soil organic C on average was slightly overestimated (7.1% with SE ± 3.3) for the limestone grassland 432 (Fig 2b), with a larger average overestimate for the acidic grassland (39.9% with SE ± 6.8). However, 433 in this latter case the variation between treatments was better captured. Despite a low r 2 value for 434 SOC (0.01), the model broadly captured the patterns we observe in the empirical data, with N 435 addition increasing SOC in the acidic and P addition increasing SOC in the limestone. However, the 436 intermediate increase in SOC with P in the acidic grassland is not captured by the model, nor is the 437 magnitude of the negative effect of LN treatment on limestone SOC. 438 Simulated magnitudes of SON are well-aligned with observations for the acidic grassland, with an 439 average error of 2.3% (SE ± 3.2), whilst SON for the limestone grassland was on average 440 underestimated by 17.8% (SE ± 3.6) (Fig 2c). The variation between treatments was better captured 441 for acidic than limestone SON but was overall reasonable (r 2 = 0.39). 442 Finally, the model overestimated total soil P (defined in the model as organic P plus sorbed P) by an 443 average of 6.0% (SE ± 4.3) for the limestone but underestimated by 54.7% (SE ± 8.0) in the acidic 444 grassland, which was the least accurately predicted variable out of those investigated (Fig 2d). With 445 only two empirical data points for TP across only two nutrient treatments, it is difficult to discern the 446 relationship between treatments and TP so an r 2 value is of little relevance here. 447 448 Figure 2: A comparison of the observed values of a) aboveground biomass carbon, b) soil organic carbon, c) soil organic nitrogen and d) total soil phosphorus from both grasslands, with simulated values from the model. The blue line represents a 1 to 1 relationship and the closer the data points are to the line, the smaller the discrepancy between observed and modelled data. All data are in grams per metre squared and all treatments for which data were collected are presented. The horizontal error bars represent the standard error of the empirical data means. The r 2 value of regression models fitted to the data give an overall indication of the direction of response of each variable to nutrient addition, hence a low value is not necessarily indicative of poor model fit

The limiting nutrient through time 452
Modelled acid grassland NPP remained N-limited from 1800 through to 2020 under most nutrient 453 treatments (Fig 3). Nitrogen deposition increased the potential NPP through time and the grassland 454 moved toward co-limitation in the LN treatment (i.e. the N and P lines were closer) but remained N-455 limited (Fig 3b). In the HN treatment, the acidic grassland shifted to P limitation as N-limited NPP 456 surpasses P-limited NPP (Fig 3c). 457 The simulated limestone grassland was also initially N-limited, but was driven through a prolonged 458 (c. 100 year) state of apparent co-limitation until clearly reaching P-limitation in 1950, solely as a 459 result of N deposition (Fig 3). In the 0N treatment, the grassland remained P-limited but the 460 potential NPP values for N and P are similar, suggesting the grassland is close to co-limitation (Fig  461   3e). The LN and HN treatment amplified pre-existing P-limitation, lowering the potential NPP of the 462 grasslands (Fig 3f, g). With the addition of P in 1995, P limitation is alleviated, and the ecosystem 463 transitions to a more productive N-limited grassland (Figure 3h). 464 Another way to interpret the extent of nutrient limitation within N14CP with specific reference to P-465 demand, is to assess the rate of P cleaving through time. These data corroborate the N and P-limited 466 NPP data, showing that in the limestone grassland, the maximum amount of cleavable P is accessed 467 by plants in the 0N, LN and HN treatments from approximately 1900 through to the end of the 468 experimental period in 2020 ( Fig S1, Table S14), highlighting its consistent state of P limitation. 469 Conversely, while cleaved P is used in the 0N treatment in the acidic grassland, it occurs at 470 approximately one third of the total rate, hence the grassland is not entirely P-limited (Fig S1, Table  471 S10). The LN treatment increases the rate of access to cleaved P and HN causes it to reach its 472 maximum value, confirming the shift to P limitation suggested by the NPP data ( Fig S1, Table S10). 473 Soil organic P cleaving does not occur in the P-treated plots of either grassland. 474  (1995). The value of the lines represents the maximum amount of productivity attainable given the availability of N and P separately. Due to a Liebig's law of the minimum approach to plant growth, it is the lowest of the two lines that dictates the limiting nutrient of the grassland and represents actual modelled productivity. Where lines share a value, it can be considered in a state of N-P co-limitation.

Modelled trends and responses to nutrient additions 477 478
The model allows the temporal trends and responses to nutrient additions to be further explored. 479 Figure 4 provides the temporal responses for the treatments, and Figure 5 a full nutrient budget for 480 the year 2020. Full data for changes in soil C, N and P and plant biomass C pools since the onset of 481 large-scale N deposition (1800 within the model) for both grasslands are included in supplementary 482 Table S15. All data used for determining responses of biomass C and soil organic C, N and P pools to 483 experimental nutrient additions are in supplementary Tables S16 (acidic) and S17 (limestone). 484 485 3.3.1. Acidic grassland 486 The modelled time series suggest that in the 0N (control) treatment for the acidic grassland, 487 background levels of atmospheric N deposition between the period 1800-2020 resulted in an almost 488 four-fold increase in biomass C, a near-twofold increase in SOC and SON and increased the size of 489 the SOP pool by almost a fifth (Fig 4). 490 Since initiated in 1995, all C and N pools responded positively to N but not P treatments (Fig 5a, c,  491 Tables S7, S8). The LN and HN treatments further increased aboveground biomass C by 36.2% and 492 61.7% (Fig 4a) and increased the size of the total SOC pool by 11.5% and 20.6% respectively (Fig 4c). 493 Similarly, the total SON pool in the acidic grassland increased by 9.7% in the LN treatment and 36.6% 494 in the HN (Fig 4e). 495 Responses of the SOP pool are in contrast to those of the SOC and SON pools, with LN and HN 496 decreasing SOP by 4.4% and 9.1% respectively, while P addition substantially increased the size of 497 the SOP pool by 76.7% (Fig 4g). Nitrogen treatments facilitated access to SOP from both subsoil and 498 topsoil, increasing plant available P and facilitating its uptake into biomass material (Fig 5e, Table  499 S9). 500 501

Limestone grassland 502
Model simulations for the limestone grassland also suggest N deposition between 1800 and 2020 503 considerably increased aboveground biomass C, SOC and SON pools (Fig. 4), but to a lesser extent 504 than in the acidic grassland. Soil organic C and SON increased by almost half and biomass C more 505 than doubled. Soil organic P accumulated at a faster rate than in the acidic grassland, increasing by 506 about a third (Fig 4, Table S15). 507 Responses of the aboveground biomass C and SOC pools in the limestone grassland differ greatly to 508 those of the acidic, declining with N addition and increasing with P addition (Fig 4). This response 509 was ubiquitous to all C pools, with declines in subsoil, topsoil and biomass C (Fig 5b, Table S11). 510 Biomass C declined by 2.4% and 7.3% with LN and HN addition (Fig 4b) and SOC declined by 0.5% 511 and 1.4% with the same treatments (Fig 4d). Phosphorus addition increased biomass C and SOC by 512 22.0% and 6.1% respectively (Fig 4b, d). 513 Nitrogen treatments increased the size of subsoil, topsoil and available N pools, but led to small 514 declines in biomass N (Fig 5d, Table S12). The P treatment slightly reduced subsoil and topsoil SON 515 compared to the control yet increased available N and biomass N, to the extent where biomass N is 516 greater in the P than HN treatment (Fig 5d, Table S12). Total SON increased by 6.4% and 15.0% with 517 LN and HN respectively and declined by 0.2% with P treatment (Fig 4f). 518 The response of the limestone P pools mirrors that of carbon, with declines in subsoil SOP, topsoil 519 SOP, available P and biomass P with LN and HN addition (Fig 5f, Table S13). The limestone grassland 520 SOP pool declined by 0.2% with LN and 0.5% with HN addition, with an increase of 20.0% upon 521 addition of P (Fig 4h). The P treatment substantially increased total ecosystem P in the limestone 522 grassland, particularly in the topsoil sorbed pool (Fig 5f, Table S13). 523 524 525 Figure 4: Time series plots of aboveground biomass C, soil organic C, N and P for the acidic (panels a, c, e and g respectively) and limestone modelled grasslands (panels b, d, f and h respectively). The vertical dashed line represents the first year of nutrient addition (1995) and marks the beginning of the experimental period. The inset subplots focus on this experimental period (1995-2020) and highlight changes occurring as a result of nutrient additions rather than background N deposition. All nutrient treatments at Wardlow are represented in all panels though not all lines are visible if they do not differ from 0N. Both grassland share a y axis. Empirical data from figure 2 are plotted on the respective panels, with the exception of panels g and h, where empirical data is incompatible with modelled data (total P versus organic P).

Figure 5:
Modelled C, N and P budgets for the acidic (panels a, c and e) and limestone (panels b, d, f) grasslands for the year 2020. Modelled sizes of C and N pools are in grams per metre squared, and P pools are presented as log n grams per metre squared. Temporary pools such as available N and P and fixed N are not presented here to avoid 'double counting' in other pools and wood litter C, N and P are not presented due to their negligible sizes.

Simulating contrasting grasslands by varying plant access to P sources 530
This is the first instance in which N14CP, and to the best of our knowledge, any other integrated C-N-531 P cycle model, has explicitly modelled P-limited ecosystems and investigated their responses to N 532 deposition and additional nutrient treatments. By using empirical data from long-term experimental 533 grasslands to drive and calibrate N14CP, we could test the model's ability to simulate two 534 contrasting P-limited grasslands, and how organic P access may affect this ability. While the purpose 535 of this work was not to explicitly reproduce the Wardlow grasslands within N14CP, by comparing 536 data from Wardlow to the simulated grasslands, we can simultaneously develop our understanding 537 of the model's representation of under-studied P cycling processes and contextualise what this may 538 mean for empirical systems such as Wardlow. 539 The model suggests that the acidic grassland was characterised by high access to organic P, with 540 comparatively low inorganic P availability, whereas the limestone grassland was the opposite, with 541 low organic and high inorganic P availability. These simulated differences could reflect the relative 542 availability of different P sources at Wardlow. As the acidic grassland formed in a hillside depression, 543 loess has accumulated, thickening the soil profile and distancing the plant community from the 544 limestone bedrock. The plant rooting zone of the acidic grassland is therefore not in contact with the 545 bedrock, and roots almost exclusively occur in the presence of organic P sources which can be 546 cleaved and utilised by plants [Caldwell, 2005;Margalef et al., 2017]. Conversely, the limestone 547 grassland soil rarely exceeds 10 cm depth, and the rooting zone extends to the limestone beneath, 548 providing plants with greater access to weatherable calcium phosphate [Smits et al., 2012]. 549 Such parameter combinations allowed for reasonable congruence between empirical and simulated 550 data, with an average discrepancy of only 6.6% (SE ± 9.1) and 1.2% (SE ± 4.4) for the acidic and 551 limestone grasslands respectively across all variables (Table S5). However, model performance 552 differed greatly between the two grasslands. For instance, the model accurately captured the 553 magnitude of limestone C, N and P data and their expected P-limited responses to nutrient 554 treatment, but was less effective at simulating the acidic grassland. N14CP did not simulate an 555 increase in biomass C or SOC with P addition in the acidic grassland, instead simulating a solely N-556 limited grassland. While this may be expected of a model that employs a law-of-the-minimum 557 approach, N14CP has a number of mechanisms to account for N and P interdependence, meaning 558 that in principle, it is capable of simulating positive responses to LN, HN and P treatment, as 559 observed in the empirical data from 2017 (section 2.2.2). 560 The overestimation of acidic C pools and underestimation of total P suggests that the model is 561 simulating that too much organic P is being accessed by plants in response to N addition and 562 transferred into plant biomass pools (Fig 2d). Few parameter sets where simultaneously able to 563 simulate the magnitude of the empirical TP pool and the positive response of biomass to N addition 564 in the acidic grassland. This may also be due to limitations in the empirical P data, as P data used for 565 calibrating P cycling were available for only two nutrient treatments and represented total soil P, not 566 organic P. While we acknowledge the technical and theoretical issues associated with distinguishing 567 between organic and inorganic P pools [Lajtha et al. 1999; Barrow et al. 2020], such distinctions 568 would help in understanding this discrepancy and likely improve the model's ability to simulate P-569 limited systems, particularly when organic P availability may be important. 570 Additionally, N14CP's representation of organic P cleaving likely underestimates the ability of soil to 571 rapidly occlude and protect organic P that enters solution. For example, inositol phosphate, a major 572 constituent of organic P, has been found to be used extensively by plants grown in sand but is hardly 573 accessed by plants grown in soil [Adams and Pate 1992]. Such organic phosphates become strongly 574 bound to oxides in the soil, protecting them from attack by phosphatase enzymes [Barrow 2020]. 575 This may be particularly prevalent in the acidic grassland at Wardlow where N deposition has 576 resulted in acidification and base cation depletion [Horswill et al. 2008], potentially enhancing the 577 formation of iron and aluminium complexes and immobilising P [Kooijman et al., 1998]. 578 In addition to physico-chemical processes reducing P availability, in P-limited grassland soils, 579 microbial processes may be dominant drivers of ecosystem P fluxes [Bünemann et al. 2012]. For 580 instance, while mineralisation of organic P may increase inorganic P in solution [Schneider et al. 581 2017], this can be rapidly and almost completely immobilised by microbes, particularly when soil P 582 availability is low [Bünemann et al. 2012]. As the model lacks a mechanism for increasing access to 583 secondary mineral P forms comparable to organic P-cleaving, and microbial P immobilisation is 584 incompletely represented for P-limited conditions, it is possible that the uptake of organic P by the 585 acidic grassland in the model is exaggerated. 586 As the model lacks a mechanism for increasing access to secondary mineral P forms comparable to 587 organic P-cleaving, the uptake of organic P by the acidic grassland is very likely exaggerated. 588 The model's inability to simulate a positive response to both N and P addition in the acidic grassland 589 may be an unintended consequence of the downregulation of N fixation by N deposition included 590 within N14CP [Davies et al. 2016b]. While this representation is appropriate [Gundale et al. 2013], 591 when N deposition exceeds fixation (as at Wardlow), fixation is essentially nullified (as in Tables S8,  592 S12), meaning deposition becomes the sole source of N to the grassland. This in effect, removes the 593 dependence of N acquisition on P availability, and could make modelling behaviour akin to N-P co-594 limitation [Harpole et al., 2011] under high levels of N deposition challenging. This suggests that 595 current C-N-P cycle models that employ a Liebig's law of the minimum can provide a broad 596 representation of multiple variables by calibrating access to both organic and inorganic P sources 597 [Davies et al. 2016b], provided the ecosystem in question's limiting nutrient leans towards N or P 598 limitation. Furthermore, where access to organic P forms is likely to be lower, as in the limestone 599 grassland, model performance may improve. This could be further explored by allowing N fixation 600 limits in the model to adapt to P nutrient conditions or by attenuating the suppression of N 601 deposition on N fixation, to represent acclimatisation of N-fixers to greater N availability [Zheng et 602 al. 2018].
Ultimately, differences in modelled accessibility to organic forms of P enabled N14CP to distinguish 604 between the two empirical grasslands, and simulate the magnitude and pattern of data with 605 reasonable accuracy, albeit with the previously mentioned caveats. 606 607 4.2. Consequences of differential P access on ecosystem C, N and P 608 While the model's estimation of P CleaveMax for the acidic grassland is likely overestimated, the model 609 experiment has highlighted that differences in organic versus inorganic P availability are a key 610 determinant of an ecosystem's nutrient limitation, and consequently, how they respond to changes 611 in anthropogenic N and P availability. For instance, while being exposed to the same background 612 level of N deposition and the same magnitude of experimental treatment, the modelled acidic 613 grassland was able to stimulate growth in response to LN and HN treatment whereas the modelled 614 limestone grassland was negatively affected by it. and HN treatment, worsening P limitation in the limestone grassland, and depleting the SOP pool in 619 the acidic. As P cleaved from organic pools is the least bioavailable within the model hierarchy 620 (methods 2.2.3), this is indicative of increasing P stress in both grasslands. While SOP declined in 621 both grasslands, the responses of available and biomass P to nutrient treatments differed markedly 622 between the grasslands. Due to the higher rate of P CleaveMax in the acidic grassland, more P 623 accumulated in thewas in plant-available forms pool and hence P does not become the limiting 624 factor under N treatments (Table S9). Conversely, available and biomass P decline under LN and HN 625 addition in the limestone grassland (Table S13), highlighting how the grassland's P CleaveMax capability is 626 insufficient to meet increased P demand. 627 Such high access to organic P sources in the modelled acidic grassland likely led it to respond to 628 nutrient enrichment in an N-limited manner, increasing productivity in response to N deposition and 629 LN and HN treatments as the model's limiting nutrient stimulated plant growth. Detrital C inputs 630 from plant biomass are the primary source of SOC accumulation within N14CP [Davies et al., 2016b] 631 and as such, changes in SOC integrate long term trends in net primary productivity in systems where 632 external nutrients are supplied. The provision of additional N in the modelled LN and HN treatments 633 therefore led to large increases in biomass accumulation and consequently, almost linearly increased 634 SOC (Fig 4c). Despite its P-limited condition under the HN treatment (Fig 3c), the acidic grassland continued to 642 accumulate biomass with N addition as the grassland's greater access to topsoil SOP (Table S9) 643 allowed it to acquire sufficient P to stimulate additional growth but not necessarily to alleviate P 644 limitation. This is consistent with the acidic grassland at Wardlow, where N treatment stimulated 645 root surface phosphatases, likely supplying more SOP to plants [Johnson et al., 1999]. Our simulated 646 acidic grassland therefore supports the hypothesis that prolonged N deposition may increase SOP 647 access to such an extent that P limitation is alleviated and growth can be stimulated [Chen et al. 648 2020]. Organic P release from SOM and its potential immobilisation, is poorly represented in models 649 and we encourage further study aimed at quantifying these processes [Chen et al. 2020 Conversely, biomass C and SOC in the modelled limestone grassland responded positively to P 654 addition, via similar mechanisms to the N-response in the modelled acidic grassland. However, in 655 contrast to the acidic grassland, N addition caused declines in limestone biomass and SOC, the 656 former of which has been observed at the limestone grassland at Wardlow [Carroll et al., 2003]. 657 Reductions in limestone biomass C (and consequently SOC) in the model are a combined result of 658 reductions in bioavailable P (Table S13), occurring via N-driven increases in stoichiometric P demand, 659 in addition to an inability to access sufficient P from the SOP pool (Table S14). Plants therefore 660 cannot meet P demand and new biomass is insufficient to replace senesced plant material, 661 decreasing net biomass C input to the SOC pool. This suggests that in P-limited limestone grasslands 662 such as at Wardlow, where access to organic P forms may be comparatively limited, N deposition 663 may worsen pre-existing P limitation and reduce ecosystem C stocks [Goll et  While N14CP is a fairly simple ecosystem model by design, it is one of few models to integrate the C, 671 N and P cycles for semi-natural ecosystems and has been extensively tested against empirical NPP 672 and soil C, N and P data it is one of the first process-based biogeochemical models to explicitly 673 incorporate P with the C and N cycles for semi-natural ecosystems, and to simulate NPP and soil C, N 674 and P dynamics, for which it has been extensively tested [ which we aimed to do in this study by using long-term experimental data from contrasting P-limited 678 grasslands. 679 N14CP's simplified representation of plant nutrient pools and plant control over nutrient uptake, is 680 largely controlled by stoichiometric demand [Davies et al. 2016a], and does not incorporate many 681 plant strategies for P acquisition [Vance et al. 2003]. Indeed, by allowing P CleaveMax to vary to account 682 for empirical data, we attempt to somewhat increase plant control over organic P uptake. We 683 acknowledge earlier that such an approach likely underestimates the ability of soil surfaces and 684 microbes to protect newly-cleaved P from plant uptake. As such, where we may expect access to 685 organic P to be high, such as the acidic grassland at Wardlow, such modelled representation of 686 plant-mediated P access may lead to unrealistic depletions in soil P and increases in biomass and soil 687 C, and we would encourage further work aimed at improving model-representation of plant controls 688 on organic P cycling [Fleischer et al. 2019]. 689 While we feel incorporating a suite of plant strategies for acquiring P would represent over-690 parameterisation, we acknowledge that a modelled equivalent to P CleaveMax for accessing inorganic P 691 forms is lacking, such as carbon-based acid exudation to increase mineral P weathering [Achat et al.  Currently, N14CP assumes C to be in unlimited supply, with its uptake by plants and consequent 700 input into soil pools controlled by C:N:P stoichiometry, hence C availability has little effect on N and 701 P dynamics within the model. Increasing atmospheric CO 2 may increase nutrient availability, as 702 plants may reallocate additional carbon resources toward nutrient acquisition [Keane et al. 2020] or 703 elevated CO 2 (eCO 2 ) may increase limitation of other nutrients such as N [Luo et al. 2004]. The 704 inclusion of eCO 2 into N14CP poses a particularly enticing research opportunity, and we aim to use 705 this study as a foundation for future work to include this process. 706

Conclusions 708
We have shown that by varying two P-acquisition parameters within N14CP, we can account for 709 contrasting responses of two P-limited grasslands of differing soil P chemistry, and with reasonable 710 accuracy. However, such coarse representation of organic P cycling in the model likely overestimates 711 the ability of plants to use newly-cleaved P and limits our ability to simulate grasslands where N and 712 P interact to control plant productivity, including the potential for N inputs to alleviate P limitation. 713 Differences in organic P access was a key factor distinguishing the contrasting responses of the 714 modelled grasslands to nutrient manipulation, with high plant access allowing the acidic grassland to 715 acquire sufficient P to match available N from chronic deposition and prevent 'anthropogenic P 716 limitation'. In the acidic grassland, N treatment stimulated plant access of organic P, promoting 717 growth and C sequestration. However, the model suggests that this is an unsustainable strategy, as 718 the SOP pool rapidly degrades, and if N additions are sustained, P limitation may return. Conversely 719 in the limestone grassland, which was less able to access organic P, additional N provision 720 exacerbated pre-existing P limitation by simultaneously increasing plant P demand and reducing P 721 bioavailability. This reduced productivity and consequently C input to soil pools declined, resulting in 722 SOC degradation exceeding its replacement. 723 We further show that anthropogenic N deposition since the onset of the industrial revolution has 724 had a substantial impact on the C, N and P pools of both the modelled acidic and limestone 725 grasslands, to the extent where almost half of contemporary soil C and N in the model could be 726 from, or caused by, N deposition. 727 Our work therefore suggests that with sufficient access to organic P, long-term N addition may 728 alleviate P limitation. Where organic P access is limited, N deposition could shift more ecosystems 729 toward a state of P limitation or strengthen it where it already occurs [Goll et al., 2012], reducing 730 productivity to the point where declines in grassland SOC stocks -one of our largest and most labile 731 carbon pools -may occur.
Data availability: Data archiving is underway with the NERC's Environmental Information Data