Towards an integrated phosphorus, carbon and nitrogen cycling model for topographically diverse grasslands

Contemporary science on how livestock influence nutrient cycling in grazing systems is limited, particularly in topographically complex (i.e., slopes and aspects) hill country landscapes. Prominent slope and aspect variation affects primary production, animal behaviour and nutrient return. Here, we embed recent scientific advancements in nutrient dynamics across complex landscapes to (1) set up a soil organic carbon (SOC) saturation function to an existing SOC and total soil phosphorus (TSP) model (Bilotto et al. J N Z Grassl 81:171–178, 2019. https://doi.org/10.33584/jnzg.2019.81.397), (2) include total soil nitrogen (TSN) dynamics, and (3) establish if the model (herein the Grass-NEXT model) can simulate the spatial and temporal changes of TSP, SOC and TSN in hill country. A long-term P fertiliser experiment with contrasting different P fertilisation levels and associated sheep stocking regimes (herein, ‘farmlets’) was used for model testing. The Grass-NEXT model predicted TSP and SOC stocks with strong accuracy and precision (model performance), and TSN with a moderate performance across farmlets [Concordance Correlation Coefficient (CCC), 0.75, 0.72 and 0.49, respectively]. Grass-NEXT model simulated TSP, SOC and TSN distribution with moderate/strong performance across slopes (CCC, 0.94, 0.80 and 0.70) and aspects (CCC, 0.83, 0.67 and 0.51). Consistent with observed data, modelled changes in TSP and TSN were greater on low slopes and eastern aspects, but no clear pattern was observed for SOC stocks. The Grass-NEXT model provides an intuitive research tool for exploring management options for increasing SOC and TSN, as well as an instrument for monitoring and reporting on nutrient dynamics in complex landscapes.


Introduction
Globally, grasslands comprise a key land use vegetation cover, accounting for 26% of the ice-free land area (with about a third of this area in managed grasslands), store up to 343 Pg soil organic carbon (SOC) to 1 m depth (about 20% of the world's total soil carbon stock) and a rate of sequestration of 0.5 Pg C year −1 (Lorenz and Lal 2018; Abstract Contemporary science on how livestock influence nutrient cycling in grazing systems is limited, particularly in topographically complex (i.e., slopes and aspects) hill country landscapes. Prominent slope and aspect variation affects primary production, animal behaviour and nutrient return. Here, we embed recent scientific advancements in nutrient dynamics across complex landscapes to (1) set up a soil organic carbon (SOC) saturation function to an existing SOC and total soil phosphorus (TSP) model (Bilotto et al. J N Z Grassl 81:171-178, 2019. https:// doi. org/ 10. 33584/ jnzg. 2019. 81. 397), (2) include total soil nitrogen (TSN) dynamics, and (3) establish if the model (herein the Grass-NEXT model) can simulate the spatial and temporal changes of TSP, SOC and TSN in hill country. A long-term P fertiliser experiment with contrasting different P fertilisation levels and associated sheep stocking regimes (herein, 'farmlets') was used for model testing. The Grass-NEXT 2021). The availability of macro-nutrients, such as phosphorus (P) and nitrogen (N), are considered major drivers of the C cycle in grasslands soils (van Groenigen et al. 2017;Poeplau et al. 2018). While some studies have revealed that sustained P fertiliser inputs to temperate pastures increase soil C sequestration (Chan et al. 2010;Coonan et al. 2019), others have shown negligible effects (Cayley et al. 2002;Schipper et al. 2011;Condron et al. 2012;Young et al. 2016;Whitehead et al. 2018;Mackay et al. 2021a) or even a reduction of SOC stocks (Sochorová et al. 2016), especially under conditions of limited soil N (van Paassen et al. 2020). Including the effect of P fertiliser on soil N availability (from biological N fixation and other sources) might provide for a better understanding of the effects of P fertiliser on the cycling of C and SOC sequestration in grazed systems (Christie et al. 2018(Christie et al. , 2020Bilotto et al. 2021).
Several grassland models each in pursuit of specific research questions have been developed to simulate P, SOC and N dynamics in grazed systems [e.g., CENTURY (Parton 1996), DayCent (Parton et al. 1998;Del Grosso et al. 2008), DairyMod (Christie et al. 2018(Christie et al. , 2020, APSIM (Holzworth et al. 2014)]. Some models focus on improving the representation of biochemical, biophysical and ecosystem processes in natural grasslands [e.g., Grassland Ecosystem Model GEM (Fitz et al. 1996)], some on describing C, N and water cycles in grazed systems [e.g., Hurley Pasture Model (Thornley 1998), PaSim (Riedo et al. 1998), Soil Plant Atmosphere System Model SPASMO (Cichota and Snow 2009)] or on simulating multiple grass species which compete for water and nutrients [e.g., DNDC (Li et al. 2012)]. Soil models such as RothC (Coleman et al. 1997;Coleman 2008) require relative stock change factors based on land use, land management, soil type and climatic region (weather data), and monthly input of plant residues and livestock excreta. As far as we are aware no model captures the transfer and spatial redistribution of P, C and N from animal excreta in grazed agroecosystems. In addition to the factors above, animal grazing and excreta return is also influenced by topography, botanical composition, and proximity to shelter, shade and water (Langworthy et al. 2018;Chen et al. 2021). All these are conducive to greater spatial variability in excreta return as the land surface slope increases, which can also increase the variability in micro-climatic and agro-ecological conditions (Cottle et al. 2016).
New Zealand hill country landscapes (slopes > 15º) cover 4.1 million ha, accounting for more than 50% of the country's permanent grasslands and are the base for the sheep and beef pastoral industry (Kemp and López 2016). Modelling efforts to capture long-term soil nutrient dynamics under hill country grazing in New Zealand have been reported (Saggar et al. 1990;Hoogendoorn et al. 2011), but none of these reports capture the complex relationships between soil P, SOC and N stocks in soil and their interaction with topography. To account for some of these factors, Bilotto et al. (2019) built a spatial, process-based P and C cycling model following the P and sulphur (S) model of Saggar et al. (1990) and the C and N models by Hoogendoorn et al. (2011). In addition to slope, the model developed by Bilotto et al. (2019) accounted for aspect [the main compass direction (East, Northwest, Southwest) that each slope faces], increased soil depth (from 150 to 300 mm) to better account for P movement down the profile (Mackay and Costall 2016), and to explore temporal changes in P and SOC stocks and flows in these systems over time. The addition of aspect, extra soil depth and time components allowed for an improved prediction of P dynamics compared with the original model of Saggar et al. (1990). The model, however, overestimated SOC stocks to a 300 mm depth (Bilotto et al. 2019) compared with observed data (Mackay et al. 2021a). The objectives of this study are to i) include a C saturation function in the SOC model developed by Bilotto et al. (2019), ii) expand the model development to include the N dynamics of the grazed system (i.e., N pools and flows), and iii) establish if the new Grassland Nutrient EXplorer and Transfer (Grass-NEXT model) model can simulate the spatial distribution and temporal changes of P, SOC and N simultaneously in hill country as influenced by contrasting P fertiliser and associated sheep stocking and topography (slopes and aspects) over time.
Three slope classes and three aspect locations were used to describe the topography of the landscapes within the farmlets: Low (LS; 0°-12°), Medium (MS; 13°-25°) and High slope (HS; > 25°) on different aspects grouped relative to the true north as East (E; 35-155°), Northwest (NW; 155°-275°) and Southwest aspect (SW; 275°-35°) (Fig. 1). Land area in each slope and aspect combination (expressed as a proportion of the total area in each farmlet) is described in Table 1.
Soil samples from each farmlet were collected from permanently marked sites to a depth of 300 mm in 2003 and in 2020. A brief description of soil sampling procedures and soil sample analyses is presented in the supplement (Online Resource 1). The effect of farmlet, topographical feature (slope class and aspect location) and sampling year (and their interactions) on P, C and N concentrations, and P, C, and N stocks were subject to analysis of variance (ANOVA) according to a split-split-plot design (Online Resource 1). For more information on soil sampling procedures and sample analysis see Mackay et al. (2021a, b).

Model development
The new Grassland Nutrient EXplorer and Transfer (Grass-NEXT) model was built using the webbased, general-purpose simulation and modelling tool Insight Maker (Fortmann-Roe 2014). The model recognises the inherent differences in a wide range of soil attributes and the influence of topography on primary production and animal behaviour, as well as how these affect the dynamics of the biogeochemical cycles of P, C and N soil stocks over time (Saggar et al. 1990;Hoogendoorn et al. 2011;Bilotto et al. 2019). The effect of aspect as it influences primary production (Lambert and Roberts 1978) and animal behaviour including grazing, camping, and nutrient return (López et al. 2003a), were added to this framework. Mackay and Costall (2016) reported that under high P fertiliser inputs, P was found in significant amounts in the 75-150 mm soil depth. This prompted the need to extend soil depth of modelling of P to 300 mm, which is also in agreement and harmonizes with the soil depth used in global reporting on changes in SOC stocks (FAO 2018).
A similar approach to that of Saggar et al. (2015) was used to account for the effect of spatial variability on nutrient excreta return from livestock. Equations that relate to slope classes were generated and adapted to aspects based on Saggar et al. (2015). In addition to those published in Bilotto et al. (2019), model parameters and equations used to calculate changes in soil P stocks (herein total soil P; TSP) and SOC stocks to a depth of 300 mm are in the supplement (Online Resource 2 and 3, respectively). Model parameters and equations used to calculate soil N stocks (herein total soil N; TSN) to a depth of 300 mm are also in the supplement (Online Resource 4).

Accounting for soil carbon saturation
The physical capacity of a soil to retain soil organic matter (SOM) is determined by the amount of C associated with the clay and silt fractions in the soil. This in turn determines the amount of newly added C in plant root and plant and animal litter that is retained after decomposition (Hassink and Whitmore 1997). As such, a C saturation function was added to the C model parameters of Bilotto et al. (2019) leading to an equilibrium of SOC stocks in response to limitations of soil properties and climate, as opposed to the continuous accumulation of SOC stocks over time. The addition of a C saturation function had been signaled, as a prior iteration of the model predicted C accumulation over time across all three farmlets (Bilotto et al. 2019), when field measurements (herein observed data) of SOC stocks have shown minimal change over the last four decades (Lambert et al. 2000;Mackay et al. 2021b).
Plant and soil respiration, and SOC stock decomposition depend on the current level of SOC stocks and C saturation level (Kemanian and Stöckle 2010;Pravia et al. 2019), according to the following: where SOC sa i−1 is annual SOC change and the components C litter , C exc and C root and rhizo represent C inputs from plant litter, C returned as dung and C from root and rhizodeposition weighted by land proportion (Table 1) for slope and aspect combination (sa) and expressed in Mg C, respectively. Model parameters and equations used to calculate soil C stocks (herein total SOC) to a depth of 300 mm are also in the supplement (Online Resource 3).
The term ks sa is the use as saturation ratio to regulate SOM turnover rate, which can be expressed as: where kx sa is the maximum undisturbed C decomposition rate, ft sa is a tillage factor, fs sa is a factor that accelerates the turnover as the saturation ratio increases, and fe sa is an environmental factor.
ks sa = kx sa × ft sa × fs sa × fe sa kx = 0.0165 d r = soil disturbance factor. Long-term pastures have a null impact, a tool or operation rated lower for d r gently disturbs the soil, while a tool rated higher aggressively breaks up and mixes soil aggregates (Kemanian and Stöckle 2010).
f clay = 0.2 , which represents the clay concentration to a 300 mm depth.
where SOC s and SOC x are current SOC concentration and maximum SOC concentration (or soil C saturation), respectively.
Assuming a SOC saturation deficit of 15 mg g −1 (top 150 mm) for soils in New Zealand pastures (McNally et al. 2017) and considering initial conditions for SOC stocks and bulk density from Mackay et al. (2021b) extended to a 300 mm depth, the SOC saturation deficit (SOC x def ) for the soils of the long-term P fertiliser and sheep grazing experiment at Ballantrae was calculated as: The SOC saturation (SOC x , in Mg ha −1 ) or maximum SOC concentration was: The term hc refers to the organic humification rate, which represents the SOC gain due to plant residue decomposition and depends on the saturation ratio. When the soil is near saturation, the fractional stabilisation of microbial detritus into SOC or hc is 0 approaching the following equation: where hx is the maximum hc.
It is worth noting that SOC deficits and potential SOC stocks can be altered by certain factors such as soil type (defining factors), climate (limiting factors) and farm management (reducing factors) (Ingram and Fernandes 2001).

Modelling nitrogen dynamics
To increase our understanding of the influence of the interplay between P, C and N processes in hill country landscapes, N dynamics were simulated using data from the most relevant national and international studies published in peer-reviewed journals (Fig. 2, Online Resources 2-4). Bowatte et al. (2006) have revealed that despite large P fertiliser applications and the inputs of legume N, New Zealand hill country pastures exhibited N deficiencies. The biological N fixation (N fix ) is a product of clover DM production, total N concentration and the proportion of legume N (Ledgard and Steele 1992;Rawnsley et al. 2019). The levels of N fix through clover and other legumes vary among pastures with different botanical composition, grazing management, topography, soil quality and climate across the country (Hoglund and Brock 1978;Lambert et al. 1986;Lambert 1987). The amount of N fix observed in unimproved hill country pastures is around 13 kg N ha −1 year −1 (Grant and Lambert 1979) to 185 kg N ha −1 year −1 for improved lowland pastures (Hoglund et al. 1979). Ledgard et al. (1987) suggested different proportionality constant values for N fix at low slopes (N fix-low slope , 0.03 kg N fixed × kg DM clover production) and sloping sites (N fix-med or high slope , 0.04 kg N fixed × kg DM clover production which were assumed for the modelling purposes (Online Resource 4). Moreover, the sum of N fixation by non-symbiotic organisms and atmospheric deposition can be as important as legume N fixation. According to Lambert et al. (1982), non-symbiotic N fixation (N non-symb ) can be assumed as 13 kg N ha −1 year −1 and the atmospheric N deposition by rain (N rain ) contribute with an extra 3 kg N ha −1 year −1 . Considered as the main soil N output, the amount of plant N uptake (N uptake ) by pasture across slope and aspect combination was calculated from annual NHA and herbage N concentration. The values for herbage N concentration per farmlet, and slope class were taken from published data (Bowatte et al. 2006;Hoogendoorn et al. 2011). The calculation process was carried out by subtracting the contribution of N fix from the herbage N pool (see Online Resource 4).
The proportion of pasture utilisation for each slope and aspect combination was used to estimate the amount of pasture N intake by the animals (N intake ) and the N returned to the soil through N plant litter (N litter ). The proportion of dietary N removed as animal product (N ap ) was assumed to decrease as the N concentration of the diet increases, while the remainder of the N intake was excreted (Ledgard et al. 2003). According to Pacheco et al. (2018), N partitioning in dung (N dung ) and urine (N urine ) will directly depend on the N concentration of the diet and total N intake. Higher N intake and herbage N content would increase N partitioned to N urine .
Ammonia (NH 3 ) volatilisation (N volatilisation ) occurs from fertilisers, soils, fresh and stored animal excreta in grazing systems . Bowatte (2003) identified N volatilisation as a major pathway of N loss from urine patches in hill country pastures; these values ranged from 21 to 51% of the N in the urine patch. This model assumed that 33% of N excreted as urine is lost by N volatilisation from hill pasture. Based on the work of Hoogendoorn et al. (2011), N 2 O emissions (N Nitrous Oxide ) were assumed proportional to inputs of N in urine and dung with a proportionality constant that varied according to slope. To reflect the lower emission rates observed on MS and HS, the emissions of N Nitrous Oxide were assumed to be 0.5% of urine deposited on LS and 0.08%, on MS and HS (van der Weerden et al. 2020). In addition, based on the low nitrification rates observed from hill country (Bowatte 2003), the model assumes that 30% of the N excreted as urine is leached from LS (< 12°), and no leaching losses occur from MS and HS (Parfitt et al. 2009;Hoogendoorn et al. 2011). Although the mineral N pool may fluctuate from day to day, net mineralisation (N Min ) in the model was calculated on an annual basis as the difference between the inputs to the mineral N pool (urine, atmospheric deposition) and the losses from that pool [N uptake , N Nitrous Oxide , N volatilisation and leaching (N leach )].

Net herbage accumulation, herbage intake and utilisation
Annual NHA (kg DM ha −1 ) in response to changes in Olsen P values of the three farmlets (NF, LF, HF), slope classes (LS, MS, HS) and aspect locations (E, SW, NW) were simulated for the period 2003-2020 using data collected on site (Lambert et al. 2014;Mackay and Costall 2016;Mackay et al. 2021a). To represent this relationship, equations were developed with a triangular distribution fitted between the upper and lower limits estimated for the 95% confidence interval. According to Lambert et al. (1983), E was the most productive aspect over the nine-year period from 1972/73 to 1980/81 (9325, 9003 and 8388 kg DM ha −1 year −1 for E, SW and NW aspects, respectively). According to López et al. (2003b), low and medium slopes have a higher probability of defoliation than steeper slopes. To represent animal behaviour, the probability of pasture defoliation for low (PD low ), medium (PD med ) and high slope (PD high ) were assumed to be 0.48, 0.30 and 0.22, respectively. To keep selective defoliation constant, the allocation of Past int was linked to the composition of slope and aspect area per farmlet. Pasture utilisation (PU, %) for each slope and aspect combination was calculated as the ratio between Past int and NHA weighed by area (Online Resource 3 and 4).

Model evaluation
The model was evaluated by comparing observed vs predicted values obtained at the same time for TSP, SOC and TSN (in Mg ha −1 year −1 ) to a depth of 300 mm, across farmlets, and within slope and aspect classes. Observed field data on soil C and N concentration, and SOC and TSN stocks for each of three slope classes and aspect locations within farmlets have been published elsewhere (Mackay et al. 2021b), and a summary of this information is presented in Table 2. The association between observed and predicted values was assessed using the statistics concordance correlation coefficient (CCC; Lin 1989), Pearson correlation coefficient (ρ), and bias correction factor (C b ), as described in Tedeschi (2006). Briefly, the CCC evaluates the degree to which observed and predicted data pairs fall on a 45° line, and provides a measure of both accuracy and precision, and it is calculated as CCC = ρ × C b , where ρ is a measure of precision and C b is a measure of accuracy. The correlation coefficient ρ (also expressed as r) is a measure of how far each observation deviates from the bestfit line, whereas C b is a measure of how far the bestfit line deviates from the 45° line through the origin (Tedeschi 2006). Values of CCC closer to 1 indicate a better performance of the model (Lin 1989).
The significance test of ρ provided by the statistical package in R was performed to decide if the correlation coefficients ρ were "close to zero" or "significantly different from zero" (R Core Team 2022). If the p-value is less than the significance level (α = 0.05), there is evidence to conclude that a significant linear relationship between observed and predicted data exists and ρ is significantly different from zero. The analysis of CCC does not consider significance test but the interpretation of the absolute magnitude of the observed CCC can be described as for Pearson correlation coefficient (ρ) with < 0.2 as poor, > 0.2 and < 0.6 as moderate, > 0.6 and < 0.8 as strong and > 0.8 as very strong correlation between variables (Akoglu 2018).

Results
Changes in observed soil phosphorus, organic carbon, and nitrogen stocks Annual application of SSP rates increased TSP stocks (Table 2 and Online Resource 1). TSP stocks showed no change (853 to 841 kg ha −1 ), to a small (875 to 946 kg ha −1 ) and large increase (1331-1620 kg ha −1 ) on the NF, LF and HF farmlets, respectively, from 2003 to 2020 (Table 2). Both slope class and aspect location had significant (P < 0.001) effects on TSP stocks. These were greatest on LS (1333 kg ha −1 ) and least on HS areas (852 kg TSP ha −1 ) and were greatest on E (1187 kg ha −1 ) compared with the other two aspects (1070 and 977 kg ha −1 on SW and NW, respectively) (Online Resource 1). Farmlet (as a combined effect of P fertilisation and livestock regime) and time (year of sampling) had minimal impacts (P > 0.5) on SOC sequestration, but topographic features such as slope and aspect had a sizeable impact (P < 0.001) on SOC accumulation (Table 2 and Online Resource 1). Similarly, TSN stocks were similar across farmlets and time (P > 0.5; 8.8 Mg ha −1 ), but again varied significantly (P < 0.001) with slope and aspect. These were greater on LS and MS (9.5 Mg ha −1 ) than on HS areas (7.4 Mg ha −1 ), and greater on E and SW (9.3 Mg ha −1 ) than on NW aspects (7.8 Mg ha −1 ) ( Table 2 and Online Resource 1).

Model performance across farmlets, slopes and aspects
Overall, the precision and accuracy of the model for exploring changes in TSP, SOC and TSN stocks across the farmlets was moderate to strong, with CCC values of 0.75, 0.72 and 0.49, respectively (Fig. 3a-c). While the model predicts TSP and SOC stocks with high accuracy (C b ≥ 0.8) and precision (ρ > 0.7) (Fig. 3a, b), the biases in predicting TSN affected CCC values (Fig. 3c). Despite the strong to very strong precision and accuracy of the model for predicting TSP stocks with slopes and aspects across the farmlets, the model overestimated TSP accumulation in the HF farmlet (blue symbols in Fig. 3a). Topographically, the model provided a range of predictions, from a strong prediction of TSP stock changes by slope [CCC, p and C b values of 0.94, 0.99 and 0.95, respectively (Fig. 3d)], to a moderate prediction of TSN stock changes by aspect (corresponding CCC, p and C b values of 0.51, 0.75 and 0.68, respectively) (Fig. 3i). The latter can be explained by an overestimation of TSN stocks in the SW and E aspects and an underestimation of TSN stocks in the NW aspect. Moderate to strong predictions were observed for SOC stocks by slope (Fig. 3e) and aspect (Fig. 3h).

Modelled phosphorus, organic carbon and nitrogen dynamics across farmlets
The contrasting annual amounts of P fertiliser applied, and associated sheep stocking regimes had a minimal effect on SOC and TSN stocks (< 0.4% change per year), with decreasing trends across the farmlets (Fig. 4). The highest P fertilisation level, which has seen a 35% increase in TSP stocks over the last 17 years (HF farmlet), limited the decline in the SOC balance and TSN losses compared to the LF farmlet, with no significant changes in the C:N ratio (Fig. 4, Online Resource 5a, 5b, 5c). Interestingly, no P fertilisation (NF farmlet) has slightly decreased SOC stocks while maintaining TSN levels; this combined effect can be seen in the slight downward trend of the C:N ratio.
Phosphorus fertiliser applied and sheep stocking regime affected NHA (8.2, 8.8 and 12.6 Mg DM ha −1 for NF, LF and HF, respectively) and nutrient dynamics of N and C (Fig. 4). Fertiliser represented 43% and 48% of P inputs in LF and HF, respectively. The higher P fertilisation rates has almost tripled P losses in the HF farmlet, and P removed by animal products (P rap ). However, the HF farmlet has almost doubled the P recovered from plant litter (P litter ) and dung (P exc ) compared with the LF and NF farmlets, respectively.
Although the highest amount of soil C allocated from roots and rhizodeposition (C root and rhizo ) was observed in the HF farmlet, the proportion of C allocated to soils slightly decreased as the P fertility and livestock carrying capacity was increased across the farmlets (52-56%) (Fig. 4). The flux of C root and rhizo into the soil has shown to be the main source of soil C stocks (62-68%), followed by litter decomposition (C litter , 20-29%) and C from dung (C exc , 9-13%). But with higher NHA, the levels of C intake , C methane , and C respired by the animals (C res-a ) increased, leaving similar amounts of C litter and C fluxes through plant and soil respiration (C res-p ), which showed an exponential behaviour (the latter accounts for 80-85% of total C losses).
The increasing levels of TSP stocks were correlated with legume N fixation (N fix ) and N mineralisation per hectare (N Min ) (Fig. 4). Moreover, higher stocking rates and NHA in the HF farmlet doubled the amount of N intake, N retained as animal product (N ap ) and N excreted (N urine and N dung ) compared with the NF farmlet. Partitioning of animal excreta into urine varied between 58-60% based on the level of N intake and N concentration of the herbage. In addition, N losses in the HF farmlet (N volatilisation , N Nitrous Oxide and N leach ) were two-fold greater than those simulated for the NF farmlet.

Modelled soil phosphorus, organic carbon and nitrogen stock dynamics across slopes and aspects
Mean modelled annual changes in TSP stocks were greater on the LS and E aspect than the other two slopes and aspects (Fig. 5a-c), but no clear pattern was observed for annual SOC changes (Fig. 5d-f). Annual changes in TSN stocks (Fig. 5g-i) showed similar trends to TSP stocks. The net loss of TSP from the MS and HS classes and NW aspect (Fig. 5) is most likely a consequence of the transfer of P in animal dung from MS and HS areas to LS areas, and the amount of dung returned to the different aspect locations, respectively. The spatial distribution of P, C and N excreted showed that 61%, 29% and 10% were deposited on LS, MS and HS areas, respectively, with a similar pattern across the three farmlets. In terms of aspect positions, about 50%, 33% and 17% of P, C and N excreted were returned on the E, SW and NW aspects, respectively (data not shown).
The simulated NHA on the different slope class areas (14.0 ± 2.1, 10.0 ± 2.6 and 6.7 ± 2.2 Mg DM ha −1 for LS, MS and HS, respectively) resulted in higher proportion of C partitioning and translocation to C root and rhizo at the HS areas (55% at HS vs 51% at LS, Online Resources 6a, 6b, 6c). Moreover, in LS and HS areas, C res-p was up to 41% higher and up to 34% lower, respectively, than the average obtained across the farmlets (Fig. 4 and Online Resources 6a, 6b, 6c). In addition, the E aspect had a large proportion of C that remained in the above ground component (45-48% as C shoot ) and the NW aspect had the lowest C res-p (24% below the average obtained across the farmlets) ( Fig. 4 and Online Resources 6a, 6b, 7c).
Despite the increase of TSP, an inverse relationship was observed in TSN levels across the farmlets (Fig. 4). Within the farmlets there were similarities between the sites of increased annual TSP deposition and TSN changes (Fig. 5). No clear trend was observed in the TSN changes between aspects. There was an excess of N returned to LS areas. Most of the annual N losses occurred in MS, with N inputs barely covering N outputs in HS areas (Fig. 5g-i). The LS areas annual N input levels were from two to three, and four times greater than those observed on MS and HS areas, respectively (Online Resources 6a, 6b, 6c). While N litter was the main source of N inputs for MS and HS (32-63%), N excreta became the main source of N in LS areas (46-76% of N inputs). The MS areas were high N fix suppliers (32-69% of N fix across the farmlets) with large values per ha (11.5-61.7 kg N ha −1 ). Interestingly, this amount of N only represented 6-19% of N inputs across farmlets for this slope class. The N losses modelled for LS were two to four times higher than those observed at MS and HS areas, respectively (Online Resources 6a, 6b, 6c). Pasture N uptake was the main source of N outputs, being more important in relative terms at MS and HS (> 90% of N outputs) than LS areas (70-80%). Moreover, more than 50% of N emissions (N volatilisation and N Nitrous Oxide ) occurred on LS areas.

Modelling integrated phosphorus, carbon and nitrogen cycles
The Grass-NEXT model was able to simultaneously simulate with high (TSP and SOC) to moderate (TSN) accuracy and precision stock change and distribution in a topographically complex (hill country) landscape, as influenced by contrasting P fertiliser and sheep stocking regimes from 2003 to 2020 (Fig. 3). Modelled changes in TSP distributed across hill country landscapes were a reasonable approximation of actual changes in TSP stocks (Table 2 and Fig. 3) and dynamics across the farmlets in terms of (i) pasture species, P uptake and P intake ; (2) the amount of pasture P litter ; and (3) the amounts of P exc spatially deposited on slope and aspect combinations. Extension of the modelling of soil P to a 300 mm depth and the simulation of a longer time-period increased the precision of the model compared with the original version of the nutrient transfer model (Saggar et al. 1990). While the nutrient transfer model developed by Saggar et al. (1990) explained 89% and 79% (p, Pearson's correlation coefficient) of the variation between slopes and slope-aspect combinations across farmlets, corresponding TSP predictions in this study explained 99% and 87% of their variation (Fig. 3).
In comparison with the previous version of the C model (Bilotto et al. 2019), the addition of a C saturation component increased accuracy and precision of simulated SOC dynamics across farmlets, and topographically within slope classes (CCC's values from 0.67 to 0.80) and aspect positions (CCC's values from 0.47 to 0.67), respectively (Fig. 3). The N component of the model consistently calculated annual TSN stock changes in the LF and HF farmlets, capturing the observed trend towards TSN depletion over time (Fig. 4, Table 2). Consistent annual improvements in soil P fertility can simultaneously increase soil P availability, NHA (Mackay and Lambert 2011;Lambert et al. 2014) and legume content in the sward (Sakadevan et al. 1993;Gillingham et al. 1998;Parfitt et al. 2005;Rawnsley et al. 2019). On the latter, legume contents sampled (expressed as a proportion of edible sward DM; mean ± SE) were 0.039 ± 0.005, 0.036 ± 0.013, and 0.125 ± 0.027 on the NF, LF and HF farmlets, respectively, during the 2016-2021 period (data not shown).
While annual TSP stock change rates and animal production were dependent on P loading rate, this effect did not translate into significant temporal increases in SOC and TSN stocks (Fig. 4). In agreement with our findings, Schipper et al. (2011) did not find evidence to support that P fertiliser loading in grazing system altered SOC or TSN stock changes. Condron et al. (2012) also reported similar SOC stocks to a depth of 1 m across varying P fertiliser rates. These findings were attributed to greater stocking rates (Schipper et al. 2011) and a faster rate of SOM decomposition associated with greater herbage quality and earthworm activity (Condron et al. 2012) from an enhanced aboveground herbage biomass attained from P addition. Often, the cycling of macronutrients such as P, C, and N in grazing system is coupled through biologically mediated soil and plant processes involved in SOM turnover (Dungait et al. 2012), such that the extent of nutrient cycling can be closely related to CO 2 fluxes (Rawnsley et al. 2016). However, as TSP is progressively increased, SOC and TSN concentrations ruled by isometric stoichiometry appear to be progressively decoupled from TSP concentrations (Cleveland and Liptzin 2007) (Online Resource 5).
The highest predicted C and N inputs were provided via litter turnover and N fix , along with large amounts of dung and urine excreted in the LF and HF farmlets (Fig. 4). Based on international studies conducted on long-term temperate grazing systems, a linear increase of SOC and TSN stocks would be expected in soils with low C levels or after crop rotations (Chan et al. 2010;Coonan et al. 2019). However, under circumstances of near C saturation, large amounts of C and N inputs in long-term pastures could be counterbalanced by an exacerbation of their losses (White et al. 2018;Stoner et al. 2021) or a significant reduction of inputs (Schipper et al. 2011), maintaining or slightly declining SOC (Di et al. 2018) and TSN stocks (Bowatte et al. 2006). On the one hand, soil C saturation suggests that even when increasing C inputs, SOC stocks have reached a maximum level, and the SOC accumulation rate will become asymptotic. In fact, the maximum observed C stocks (Table 2) fluctuated around the plateau predefined in the model for SOC stocks (~ 150 Mg ha −1 ). During this process, the overall rate of plant respiration and OM decomposition in the soil equals the annual input rate of organic material (Stewart et al. 2007;Stoner et al. 2021). This is consistent with the faster dung disappearance rate found on HF sites by Vibart et al. (2021). On the other hand, any increase in legume content most likely had a small effect on the TSN stocks (Bowatte et al., 2006). Largely driven by increases in pasture intake and total nutrients excreted on pasture, higher livestock stocking rates to maintain similar grazing pressure and residual biomass tend to a shift towards a new equilibrium (Parsons et al. 2013), often with larger N losses to water and the atmosphere (Parfitt et al. 2009;Luo et al. 2019), C released as methane and loss of ground cover (Hoogendoorn et al. 2011;Harrison et al. 2012a, b) (Fig. 4). Therefore, the excess of P fertilisation co-limited by N application and high SOC stocks can lead to a decoupling of C:N:P ratios (Online Resource 5).
Impacts of topography, animal behaviour, and vegetation on soil phosphorus, organic carbon, and nitrogen stocks While disproportionate TSP and TSN accumulation was observed in LS and E and to a lesser extent SW aspects with nutrient depletion in HS and NW aspects, no clear pattern was observed for annual SOC changes (Fig. 5). The topographical features in hill country landscapes, slope class and aspect location, strong and divergent influences on TSP, SOC and TSN balances, can be attributed to multiple factors. Ruminant grazing and excretal behaviour show that a disproportionate amount of livestock excreta is deposited on LS areas (Saggar et al. 2015). Often, the nutrients taken up by livestock grazing HS pastures are transferred through urine and dung to LS and camping sites (Rowarth et al. 1992;Hoogendoorn et al. 2011). Soil P is an immobile anion generally transferred in animal dung and it is not subject to large losses by leaching or through atmospheric emissions like N or C. In agreement with Rowarth et al. (1992), higher P fertiliser levels and livestock carrying capacity have exacerbated animal P transfer from HS to LS in our study. Inversely related to the prevailing wind direction, sheep tend to camp on E and S aspects with larger amounts of excreta inputs (urine and dung), contributing to higher TSP and TSN stocks which allude to lower soil C:N ratios and a higher potential for mineralisation (Lambert and Roberts 1978). Temperature and moisture regimes by aspect also have a major influence on soil nutrient status (Lambert and Roberts 1976). As it was implicitly assumed in the model, moderate net radiation, soil temperature and potential evapotranspiration in E aspects create the best environment for NHA (Lambert and Roberts 1976;Lambert et al. 1983). With the addition of a SOC saturation function across slopes and aspects, the nutrient distribution by grazing animals is more influential in the redistribution of TSP and TSN than SOC (Fig. 5, Online Resources 6a, 6b, 6c). In line with our spatial analysis, Saggar et al. (1999) reported higher C respiration rates from soils in low slope areas. As well, lower SOC turnover rates have been suggested from observed slower dung disappearance rates on NW-facing slopes ). The greater capacity for C turnover of plant litter, root biomass and animal dung (C res-p ) with increasing P fertilisation levels may be associated with considerable changes in soil microbial biomass and fungal population size (Parfitt et al. 2010), and earthworm density, abundance and activity (Schon et al. 2019). Therefore, larger amounts of C inputs on low slopes created by higher P nutrient transfer and livestock carrying capacity were quickly vanished by an accelerated nutrient cycling and higher C losses (Mackay et al. 2021a). According to Bowatte et al. (2006), hill country pastures suffer a chronic N deficiency caused in part by the disproportionate animal N transfer to LS, conducive to high nitrification rates and subsequent N losses. The actual losses would vary according to slope, soil porosity and stocking density. The TSN depletion in high fertility farmlets was boosted through higher herbage N content leading to more N intake converted into urine (Ledgard et al. 2003). The higher partitioning of N excreta into urine was responsible for the higher N losses through N emissions to air and water in LS areas (Online Resources 6a, 6b, 6c), in line with reports from previous studies . To increase TSN stocks by stimulating clover growth and N fix , P fertilisers have been frequently applied at higher levels than recommended. However, this study has shown an inverse relationship between TSN and TSP stocks across the farmlets (Fig. 4). Recently, large-scale scale research in New Zealand has focused on urine patches on lower slopes as the main source of N losses to the atmosphere from sheep and beef hill country farms (ammonia volatilisation and N 2 O emissions) (Saggar et al. 2015;van der Weerden et al. 2020), overlooking the multiple impacts generated by animal behaviour and nutrient dynamics on the remaining hill country landscape. Interestingly, even though the large N fix levels occurred in MS, these slope areas were the main source of N losses across the farmlets. In combination, the large proportion of N transfer to low slope areas (61% of N excreta return in LS), the significant amount of plant N uptake to sustain the level of NHA in MS and the relative high N losses compared to HS, exerted a strong mining process in this topographical site (Bowatte et al. 2006;Hoogendoorn et al. 2011). Therefore, soil P fertilisation far in excess of plant growth needs could be considered as a debatable tool to increase SOC stocks in locations with high C concentrations and may potentially deplete TSN levels (Poeplau et al. 2016) by altering its spatial distribution.

Limitations of this study
Although the results obtained from Grass-NEXT model are a reasonable approximation of the TSP, SOC and TSN dynamics across the farmlets, several assumptions had to be made. These included (a) published probabilities of pasture utilisation in different slope classes based on selective grazing (López et al. 2003a, b), (b) 60% of the organic matter ingested was digestible and 40% of C from pasture intake ends up in soil, (c) the incorporation of a C saturation component in the model adapted to hill country conditions has increased accuracy and precision to simulate soil C dynamics but new refinements would include site specific information about soil respiration and SOC decomposition, (d) to better describe SOC in the profile, the main constraint is the absence of multiple SOC pools and depths as described in models like RothC (Coleman et al. 1997;Coleman 2008), (e) C root and rhizodeposition were solved by equations using limited datasets in number and scope (Saggar et al. 1997(Saggar et al. , 1999, and f) herbage N content was predefined by slope and farmlet according to relevant domestic research.
Stronger links between C saturation and N mineralisation, and the inclusion of a microbial pool, could improve TSN estimations in the model (White et al. 2014). Our study has shown lower net mineralisation with slightly lower C:N ratio in the NF farmlet compared to the LF and HF farmlets ( Fig. 4 and Online Resource 5). The C:N ratio is often used for estimating net N mineralisation rates but limitations in its predictive capability have been reported (Bonanomi et al. 2019). Lower N mineralisation rates can occur even at low C:N ratios due to low C accessibility by the soil microbiome, limiting decomposition. Further work is required to obtain additional quantitative data for, and examine the sensitivity and behaviour of, each of these variables. Emerging technologies including GPS monitoring and remote monitoring sensing techniques could improve the quantification of grazing activities and excreta deposition (e. g., Plaza et al. 2022).

Conclusions
The process-based Grass-NEXT model provides a basis for accounting, monitoring and reporting on changes in total soil P, organic C and N stocks in response to current management practices in complex grazed systems. The inclusion of aspect, in addition to slope, and the ability to simultaneously simulate soil P, C and N change in a topographically complex landscape has expanded our understanding of the dynamics of these three nutrients in hill country pastures. The addition of a C saturation function improved the model's ability to simulate SOC stocks change in a topographically challenged landscape (i.e., different slope classes and aspect locations). In agreement with observed data, the contrasting annual amounts of P fertiliser applied and associated sheep stocking regimes were reflected in increasing TSP predicted, but had a minimal effect on SOC and TSN stocks. Further developments of the model are proposed to explore the interactions between P and N fertiliser use efficiency in hill country, the contribution of legumes through biological N fixation, the spatial variability in greenhouse gas emissions and the impacts of predicted future climates on nutrients dynamics in topographically complex landscapes.
Acknowledgements Model development and analysis, and preparation of the manuscript was partly funded by the Sustainable Land Management and Climate Change fund of the Ministry for Primary Industries and the AgResearch Strategic Science Investment Fund. FB was funded by the New Zealand Government through the LEARN Awards programme, and hosted by AgResearch Grasslands, Palmerston North.
Funding Open Access funding enabled and organized by CAUL and its Member Institutions.

Conflict of interest
The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.