Impact of two centuries of intensive agriculture on soil carbon, nitrogen and phosphorus cycling in the UK Science of the Total Environment

Roth-CNP model estimates C, N and P cycling within the UK agriculture for 1800 – 2010. • Simulated crop yields were comparable to the yields of UK's agricultural statis- tics. • Simulated SOC stock decreased under arable and increased under improved grassland. • Simulated N and P losses increased under both arable and grasslands. • Results shows the effect of local agriculture in a larger context of space and time. This paper describes an agricultural model (Roth-CNP) that estimates carbon (C), nitrogen (N) and phosphorus (P) pools, pool changes, their balance and the nutrient ﬂ uxes exported from arable and grassland systems in the UK during 1800 – 2010. The Roth-CNP model was developed as part of an Integrated Model (IM) to simulate C, N and P cycling for the whole of UK, by loosely coupling terrestrial, hydrological and hydro-chemical models. The model was calibrated and tested using long term experiment (LTE) data from Broadbalk (1843) and Park Grass (1856) at Rothamsted. We estimated C, N and P balance and their ﬂ uxes exported from arable and grass-landsystemsona5km×5kmgridacrossthewholeofUKbyusingtheareaofarableofcropsandlivestocknum- bers in each grid and their management. The model estimated crop and grass yields, soil organic carbon (SOC) stocks and nutrient ﬂ uxes in the form of NH 4 -N, NO 3 -N and PO 4 -P. The simulated crop yields were compared to that reported by national agricultural statistics for the historical to the current period. Overall, arable land in the UK have lost SOC by − 0.18, − 0.25 and − 0.08 Mg C ha − 1 y − 1 whereas land under improved grassland SOC stock has increased by 0.20, 0.47 and 0.24 Mg C ha − 1 y − 1 during 1800 – 1950, 1950 – 1970 and 1970 – 2010 simulated in this study. Simulated N loss (by leaching, runoff, soil erosion and denitri ﬁ cation) increased both underarable(


H I G H L I G H T S
• Roth-CNP model estimates C, N and P cycling within the UK agriculture for 1800-2010. • Simulated crop yields were comparable to the yields of UK's agricultural statistics. • Simulated SOC stock decreased under arable and increased under improved grassland. • Simulated N and P losses increased under both arable and grasslands. • Results shows the effect of local agriculture in a larger context of space and time.

G R A P H I C A L A B S T R A C T
a b s t r a c t a r t i c l e i n f o

Introduction
Agriculture in the United Kingdom (UK) has a long history of human settlement and development which dates back to 6000 years ago when humans began domesticating plants and animals in Neolithic times (Edwards and Hirons, 1984;Woodbridge et al., 2014). By 900-700 BCE, settled agriculture was established in the UK with crop rotations, pasture and coppiced woodlands. By 100-350 CE, natural forest was largely cleared with large estate-based farming systems with cattle, sheep and arable production. The UK's countryside then further changed dramatically with the majority of the population living in small farmsteads under subsistence farming. By 1300 CE, increasing demand for food brought the subsistence farming system under huge pressure because of increasing population as the land area available for agriculture was already in use. However, between 1300 and 1800 average crop yields increased in the UK due to improvements in crop management such as mixed husbandry (by combining crop and livestock), grass and arable rotation, crop rotation by including fallow and legumes leading to a British agricultural revolution during 1700-1850 (Allen, 1999;Allen, 2008;Overton, 1996). With the industrial revolution in 1850s, technological improvements also happened in the agricultural sector, for example, switching from draught animals to machines in early 1900s (Whetham, 1970). Much of the agricultural growth during this period came about as a result of increases in the area of crops and grass, which peaked in mid 1880s. After this, the agricultural area underwent a steady decline as farms became more intensive and the availability of labour diminished. During the second half of the 20th century (Musel, 2009), agricultural intensification driven by new high yielding varieties, mineral fertilizer application, chemical pest control and improved methods of cultivation (Marks and Britton, 1989) led to increase in agricultural production many-fold. Per-hectare yields of wheat almost tripled whilst barley, potato yields and milk yields per cow more than doubled (DEFRA, 2014;Marks and Britton, 1989). The total cattle population increased sharply after the middle of 20th century although there has been a decline since 1974. About 170 million ton of animal excreta (slurry) are produced annually in the UK. In terms of farm inputs, mineral nitrogen (N) fertilizer used in the UK increased five times between 1950 and 1978 (Cooke, 1980). Greater use of N and P fertilizers during this period has led to an increased loss of these nutrients into our rivers and ground water through leaching, runoff (Hood, 1982;Hooda et al., 2000), and increased atmospheric emissions of ammonia, nitrous oxide and other reactive N compounds. Agricultural land contributes 70% and 28% of the N and P load to the UK waters (Hunt et al., 2004;White and Hammond, 2007). Losses of these nutrients are associated with excessive or poorly timed applications of N or P or both (Dungait et al., 2012). Pretty et al. (2000) calculated the annual external cost of agriculture for the UK in 1996 as £2343 M (£208/ha), with the major costs associated with contamination of drinking water by pesticides, nitrate and phosphate and increased greenhouse gas (GHG) emissions, soil erosion and organic carbon losses.
Numerous spatially-variable, interacting factors such as land-use, vegetation type, weather, catchment topography and total nutrient inputs over time determine the nutrient stocks and fluxes at a farm, landscape or catchment scale. For example, nutrient concentrations in groundwater under agricultural land have been found to be several times higher than that under semi-natural vegetation (Nolan and Stoner, 2000). Growing vegetables and crops such as potatoes and oilseed rape intensively has led to high rates of nitrate leaching (Stuart et al., 2011). Nutrient concentrations in ground water have been found to be highly variable and related to changes in the weather (Rozemeijer et al., 2009) and increased as a result of land-use change (Whitmore et al., 1992). There is a strong influence of catchment slope on water quality due to slope-dependent seasonal waterlogging, which determines the fate of dissolved substances produced within and moving through the catchment (D' Arcy and Carignan, 1997). Temporal dynamics of these nutrients depend on the relative occurrence of the nutrients in different pools at different points in time. Nutrients are retained during the dry summer months as a result of bioaccumulation and adsorption in case of P, and during the wetter autumn to spring periods, these nutrients are released and transported from the floodplain into the river channel (Bowes et al., 2005).
Understanding the processes that have led to the build-up of C, N, and P in soil, ground water and surface water from the past to the present is essential to understand how to manage the supply and utilization of these nutrients into the future. This will contribute to the long-term goal of achieving a sustainable agricultural system by increasing or maintaining crop yields whilst minimising impacts on other ecosystem services (Powlson et al., 2011). It is also important to understand how these nutrient cycles (between atmosphere, terrestrial ecosystems including agriculture and hydrological systems) operate at large spatial scales across the whole UK in response to climate change and management options. A model that can summarise essential processes of soil and plant growth and their interactions and that can be applied over long timescales with readily-available driving data (climate, land-use, nutrient inputs) is essential to investigate the temporal and spatial responses in soil macronutrients at the national scale. Such a study should help both farmers and policy makers to see the effect of agriculture at the local scale in a larger context of space and time.
There are many agroecosystem models in the literature that can simulate the C, N, and P cycling under crop and grassland systems at the field scale EPIC (Jones et al., 1991); DNDC (Li et al., 1992), APSIM (McCown et al., 1996); DAYCENT (Parton et al., 1998); CropSyst (Stöckle et al., 2003); DSSAT (Jones et al., 2003). A recent review assessed the comprehensiveness of underlying processes in nine widely used C and N models found that one of the major weakness of these models is their scalability over time and space (Brilli et al., 2017). To describe soil C and N dynamics at higher spatial and temporal scales we need models of lower complexity and with longer timesteps (McGill, 1996;Manzoni and Porporato, 2009). In the literature, we can find many process-based models that were applied at national or global scale to estimate C stocks for terrestrial systems (Al-Adamat et al., 2007 (Century and RothC); Kamoni et al., 2007 (Century and RothC); Ogle et al., 2003(IPCC method);Smith et al., 2010 (ECOSSE); Tian et al., 2015 (MsTMIP); Wang et al., 2017 (RothC). Many studies estimate N and P balances at the national or global scale using spatially explicit data on crop area, livestock population, fertilizer input rates and empirical models to estimate nutrient stocks and fluxes (Smil, 1999;Lesschen et al., 2007;Liu et al., 2010;Pathak et al., 2010;MacDonald et al., 2011;Bouwman et al., 2013;Cui et al., 2013). However, only a few studies use process models to simulate coupled C, N and P models to estimate the stocks and fluxes of these macronutrients for the terrestrial systems at the national or global scale (Wang et al., 2010;Zaehle, 2013). However, in all these models, the concept of the cropping system is often simplified by omitting details of crop management, such as rotation (Leenhardt et al., 2010). In this study we use a simplified agricultural model that aggregates processes at a monthly timestep as compromise between detailed process modelling and computation time when run at 5 × 5 km grid for the UK. We chose RothC (Coleman et al., 1997) model extended for N and P (Coleman et al., 2017) as the model to describe the C, N and P dynamics because of its simplicity and its applicability in the agricultural system in the UK. Plant (both crop and grass) growth processes were based on the LINTUL model (Shibu et al., 2010;Wolf, 2012) and were simplified to suit the monthly timestep of RothC. Atmospheric and hydrological models are run at finer timesteps and atmospheric N deposition, evapotranspiration, runoff, drainage and erosion scaled up to the monthly timestep.
This paper estimates C, N and P pools, pool changes, their balance and the nutrient fluxes exported from arable and grassland systems in the UK during the historical to current period (1800-2010) using an agricultural model that was developed as part of an integrated model to analyse and simulate long-term and large-scale (LTLS) interactions of C, N and P in the UK land, freshwater and atmosphere (http://www. ltls.org.uk/). This integrated model is referred here as LTLS-IM (Bell et al., in prep), which loosely couples terrestrial (semi-natural and agricultural), hydrological and hydro-chemical models and driven by atmospheric deposition (Fig. 1). Our emphasis in this paper is to present this integrated modelling approach with a major focus on estimates of UKwide historical yield, SOC changes and nutrient fluxes.

Methodology
This study focuses on the CNP stock changes in soil in the agricultural system driven by the environmental variables, plant and livestock management. At a spatial scale of a 5 × 5 km grid, we did not consider field or farm scale dynamics or exchange of biomass or nutrients between the grid cells. We assume that lateral flow of nutrients across farms within a cell or between cells will be a 'continuous' process and will be balanced by the influx and the outflux within each cell. The agricultural model referred to as Roth-CNP model was developed by simplifying the Landscape Model (LM) (Coleman et al., 2017) to an appropriate level of detail. The LM which works on a daily timestep, simulates the biophysical processes of an agroecosystem at the field/farm scale taking into account the spatial interactions between the fields or farms across a landscape. The Roth-CNP model presented here aggregates the essential processes within the LM on a monthly timestep without any spatial interactions between the spatial units. We briefly describe here the main features of Roth-CNP model together with any major changes from the LM (Appendix I). In Roth-CNP, we use the same parameters as the LM but adapting for a monthly timestep. We tested the Roth-CNP model using the data from Broadbalk and Park Grass long-term experiments (LTEs) at Rothamsted (http://www.era.rothamsted.ac.uk/), South-East England before undertaking the historical simulation for a continuous period from 1800 to 2010 for the whole of the UK. For these simulations, the whole land area in the country was divided into 5 km × 5 km cells on a square grid with improved grass present in in 91% and arable land in 76% of the grids cells (see SI, Fig. S2.1).

Model description
The Roth-CNP model ( Fig. 1) has two major subunits: soil and landuse. In the soil module, the soil profile (can be of any depth, but in this study, it ranges from 30 cm to 150 cm) is divided into three layers. Depths of soil layers can be variable, but for this application the first and second layers were set to 15 cm each to enable a spatial comparison of CNP pools as most of the soil management activities affect the top 30 cm. The depth of the third layer is variable depending on the actual soil profile, which varies spatially across the UK. The soil unit consists of organic C, N and P, mineral N and P modules. Variables such as actual evapotranspiration (AET), soil drainage, runoff and soil moisture are treated as inputs that are calculated by a hydrological model, which is a simplified version of the G2G model (Bell et al., 2009). However, potential evapotranspiration (PET) for each landuse was estimated in a crop module (as it varies with crop type and developmental stage of the crop) based on the Penman's method (Penman, 1948). The PET estimated by the crop model was compared to the PET estimated by the hydrology model (using MORECS PET for grass assuming variable leaf area index (LAI) for summer and winter (Hough and Jones, 1997) for a few selected sites and were found to be comparable (not reported here) with a difference in PET of up 5% in winter to 12% in summer. The hydrology model within the LTLS-IM calculates components of the water balance (runoff, drainage, AET and soil moisture) for each 5 × 5 km grid-cell in the UK. Soil moisture for the entire profile (mm of water/profile depth) was used to estimate moisture content in each soil layer within the Roth-CNP model. Soil organic carbon dynamics inherited from the RothC model has been described elsewhere (Coleman and Jenkinson, 1999;Smith, 2000;Jenkinson and Coleman, 2008). The model was extended for organic N and P with similar pool structures as that for carbon determined by the C/N or C/P ratios of the incoming organic materials for Decomposable Plant Material (DPM), Resistant Plant Material (RPM), and a fixed C/N or C/P ratios for Microbial Biomass (BIO) and Humified Organic Matter (HUM) (Coleman et al., 2017). Additional temporary pools of dissolved organic carbon (DOC), nitrogen (DON) and phosphorus (DOP) were created in the model in order to estimate the loss of dissolved organic C, N and P that enters soil solution. These pools are not linked to the main C, N, P pools and become active only when manures are added to the soil. In agricultural soils, added organic amendments such as farm yard manure (FYM), slurry and other animal manures are the major sources that contribute to DOC (Bhogal et al., 2010). Since we could find little information on the export of DOC from soils under agriculture, we assume that soil organic carbon (SOC) itself contributes only a negligibly small amount to DOC and therefore, its loss from agricultural lands was ignored in this study. In the model, we assume that when organic substrates are added, a fraction (FYM-4.6%; slurry-51%, and poultry manure-6.6%) of these goes directly to the DOC, DON and DOP pools (Bhogal et al., 2010) and is lost by leaching and/or runoff immediately before the reminder enters the SOC, SON and SOP pools.
Mineral N and P species exist in single (vertically integrated) stores without partitioning them between different soil layers to co-exist with the dynamics of soil water which estimates the water balance for the whole profile. All the N and P mineralised from soil organic matter (SOM 1 ) in the three soil layers is transferred to these mineral nutrient stores. Mineral N consists of NH 4 − N and NO 3 − N pools and mineral P includes available and fixed pools. Mineral N dynamics comprises N inputs (through atmospheric deposition, biological N fixation, fertilization), transformations (nitrification and denitrification) and losses (through plant uptake, denitrification, runoff, leaching and erosion) (Coleman et al., 2017). Similarly, P dynamics comprises P inputs from fertilizers, chemical P fixation and release, crop uptake, runoff, leaching and erosion. P contribution from weathering is not simulated in the model separately, but it is considered as a part of the fixed pool. Dynamic processes leading to an equilibrium between P in fixed and available pools are described elsewhere (Coleman et al., 2017).
The rate of nitrification and denitrification depends on the relative nitrification (0.99 month −1 ) and denitrification rates (0.20 month −1 ), soil temperature, moisture and pH (Coleman et al., 2017). In the model we assume that biological N fixation (BNF) occurs only in grassland systems and on an average about 30% of grassland is a leguminous clover mix and can fix N biologically (Lüscher et al., 2014;Sanderson et al., 2013). In the model, BNF rate is calculated as a function of potential maximum N fixation rate and the rate modifying factors for temperature (f T ),soil moisture (f m ) and inorganic N (f N ) (Liu et al., 2013).
where Nfix rate and Nfix max are the actual and maximum rates of BNF (g N m −2 month −1 ). Potential maximum BNF rate depends on the live shoot biomass (g DM m −2 ), fixation rate per unit standing biomass (g N g −1 DM month −1 ) and root growth rate (g DM month −1 ). See Liu et al. (2013).
Increases in mineral N (NH 4 N and NO 3 N) concentration reduce the BNF rate in the model and we assume N that is fixed is directly transferred to the NH 4 N pool.
Mineral N and P losses occur either through runoff (in water phase) or through soil erosion (particulate) and leaching. Loss of these nutrients through runoff depends on both the nutrient (NO 3 N, available P) concentration (kg mm −1 ) at the surface (calculated as a function of depth) and the runoff (mm of water month −1 ). Since nutrient distribution in the profile may depend on the soil water and other soil profile characteristics (e.g. soil organic matter distribution, P weathering) it will be difficult predict their distribution in the profile over time. We assume nitrogen largely follows the pattern of SOC distribution with more NO 3 -N on the surface compared to the subsurface. Therefore, to estimate the NO 3 -N content for the 15 cm layer, we used an exponentially decreasing function, which will decrease with increase in soil depth. Leaching depends on the nutrient concentration (kg mm −1 ) in the soil solution and the drainage rate (mm of water month −1 ). Estimated rates of runoff and drainage were input from the hydrology model (see Section 2.2.5).
A generic plant growth model, which uses the light use efficiency (LUE, g dry matter MJ −1 ) based approach (Monteith, 1990;Monteith and Moss, 1977) is used to simulate crop and grass growth within the landuse module. The rate of biomass production depends on the incoming solar radiation in terms of photosynthetically active radiation (PAR, i.e. 50% of the global radiation), crop/grass specific LUE and growth affecting factors such as moisture and nutrient stresses (Coleman et al., 2017). The biomass formed is partitioned between roots, stem, leaves and storage organs based on the development stage (DVS) as described by Wolf (2012). In principle, crop phenology is expressed in terms of crop development stage (DVS), which is a function of temperature sum or growing degree days and includes the effect of vernalisation and/or photosensitivity of the crop (deVries, 1989), which are variety specific and may vary across the country. As the model works on a monthly timestep, and the flowering and maturity of the crop falls within a given month for a given crop across the whole country, we used a simple growth function to represent the DVS for each crop. We calculated DVS for each crop by applying the Landscape model for Rothamsted site for several years  and generated a general growth curve for each crop (Fig. 2). For grass, we assume the plant remains in vegetative phase (DVS b 1) throughout its growing period because it is continuously grazed or cut with sufficient frequency.
Crop specific parameters for different crops were taken from the model database (http://models.pps.wur.nl/glossary/l). For older varieties of wheat (developed before 1970), a few parameters were changed by calibration (see the Supplementary information (SI 5) for more details). Insufficient water and nutrients (N & P) lead to stress that affects crop growth and reduces the biomass production and yield as described by Coleman et al. (2017). The crop is harvested at maturity (when DVS = 2) and the crop yield is the weight of storage organ (g m −2 , which could be grains, seeds or tubers). Straw or trash yield depends on the stem biomass and the method of harvesting. During the historical period (prior to 1950), we assume that crops were harvested manually and were left with less residue (15% of stem biomass) compared to the current period when machines were used (30% of stem biomass). The grass model differs from the crop only in the assumption that grass is perennial in growth and is managed differently by allowing livestock grazing or frequent cutting.
A reasonably good agreement in the model results to the measurements from Broadbalk and Park Grass (see SI 5 for model calibration and testing) for plant yield, SOC, total N and NO 3 N leaching over the last 160 years with a relative RMSE of 4-71% (see SI, Tables S5.2 and  S5.3) provides confidence in Roth-CNP's ability to estimate crop and grass yields, SOC and SON for the historical to current period of 1800-2010.

Model inputs
Input data driving the model come from multiple sources and are described in the following sections. Some of these driving data come from readily available datasets (weather and soil). Other inputs include time-dependent datasets that have been created specifically for the LTLS-IM (landcover, landuse, livestock population, atmospheric deposition, fertilizer and manure input rates) and outputs from the loosely coupled models used either to initialise Roth-CNP (N14CP (Davies et al., 2016;Tipping et al., 2012;Tipping et al., 2017): provides initial states for soil C, N and P) or to provide dynamically-changing variables (the hydrological model provides actual evapotranspiration, soil moisture, erosion and runoff).

Atmospheric deposition
Atmospheric N (NH 4 -N and NO 3 -N) input to arable and grassland systems were estimated for different time slices: 1800,1900,1950,1970,1990 and 2010 (see SI 1) at a 5 km × 5 km grid resolution across the UK, using land-cover dependent deposition velocities. Nitrogen deposition values for each land cover type in each grid square were interpolated for the whole period within these time slices to calculate the annual deposition rates and were averaged to calculate the monthly input rates. Nitrogen deposited in the forms of NH 4 -N and NO 3 -N was added directly to the mineral NH 4 -N and NO 3 -N pools, respectively. Atmospheric P deposition is ignored in this study as P deposition is insignificant in the UK as most of the atmospherically transported P is natural in origin (ocean, tropical forests, peatland, dust from deserts) and are redeposited to its origin (Tipping et al., 2014).

Weather
For weather, data from several sources were combined to derive an observation-based dataset from 1800 to present. For rainfall (mm month −1 ), daily observations, which are available back to the 19th century, were used, although network coverage ranges from only 2 rain gauges in 1853 to thousands in the late 20th century. National daily rainfall estimates for each a 5 km × 5 km UK grid square were derived from any daily observations available for the period 1853 to 1910. From 1910 onwards, gridded 1 km resolution rainfall observations from CEH-GEAR (doi:https://doi.org/10.5285/5dc179dc-f692-49ba-9326-a6893a503f6e) were used. No observed daily rainfall values were available prior to 1853, so daily rainfall for a median year (1904) was assumed to be representative for the period from 1800 to 1853.
For all other weather variables (temperature (°C), shortwave radiation (W m −2 ), wind speed (m s −1 ) surface pressure (Pa) and specific humidity (kg kg −1 ), monthly mean values were used. These were obtained from the WATCH forcing dataset (http://www.eu-watch.org/) for the period 1800 to 1970 (data for the "median" year 1904 was assumed to be representative for the pre-1901 period), and similar data from the UK Met Office (http://www.metoffice.gov.uk) were obtained for the later period 1970 to 2010. The Met Office dataset provided observation-based estimates of minimum and maximum temperature (°C), sunshine hours (h), wind speed (m s −1 ) and vapour pressure (hPa). Prior to 1901, these data were not available for the whole country, and again, data from a nominal year (1904) is used. Note that although the WATCH forcing dataset provides weather estimates from 1901 to 2001, prior to 1958 these consist of statistically re-ordered data from 1958 to 2001, and as such do not provide a historically accurate weather record for the earlier period. In the UK, only rainfall data are available at a daily time-step with national coverage before 1958 (CEH-GEAR) and these have been used in preference to WATCH rainfall estimates. Any inconsistency incurred through the use of CEH-GEAR rainfall is assumed to be small compared to the use of the statistically correct but historically-inaccurate WATCH rainfall series. Short wave radiation in W m −2 was converted to MJ m −2 . Surface pressure and , which is a function of temperature sum or growing degree days and includes the effect of vernalisation and/or photosensitivity of the crop, estimated for different crops and grass as a function of their growing months (for winter wheat, 1-11 growing months = October-August; for potato, 1-5 growing months = April-August; for spring barley, 1-6 growing months = March-August; for Oilseed rape, 1-11 growing months = September-July; for fodder maize, 1-4 growing months = May-August, for grass, growing months are indefinite. Growing months are based on MAFF (1998)). specific humidity were used to calculate vapour pressure (KPa) (Nievinski, 2009).

Land cover and landuse
A land cover history for the UK was constructed using contemporary landuse datasets and the few historical maps available (see SI 2). Livestock populations and agricultural landuse data were estimated for four time slices: 1900, 1950, 1970 and 1990 and assumed constant between these dates. Five major crops (winter wheat, spring barley, oil seed rape, potato and fodder maize) were selected, which represented five major groups of crops (winter cereals, spring cereals, Oil seed crops, tuber crops and fodder crops) in the UK. The area under each of these crops represented the sum of the total area of all the crops within each of these groups. For example, the area under winter wheat represented the total area under winter wheat, winter barley and winter oats. Similarly, spring barley represented the area under both spring barley and spring wheat. Area under potato represented the area under potato and sugar beet and all the fodder crops under fodder maize. Livestock classes considered were beef, dairy, sheep, pig and poultry. Estimates were based on historic agricultural census data and were distributed using the AENEID model (Dragosits et al., 1998;Hellsten et al., 2008) explained in SI3 with the land cover data summarised in SI 2.
The arable area in each grid cell was assumed to grow up to a maximum of five representative major crops: winter wheat, spring barley, potato, oil seed rape, and fodder maize depending on their presence or absence in that grid cell.
Under improved grass (i.e. anthropogenically managed grassland), four types of grass land management: dairy, beef, sheep and silage (ungrazed) were simulated according to the livestock population at that location.
To estimate the area under each of these livestock management systems (A i ), we used the livestock numbers in each grid (N i ) and the standard stocking rate (D i , animals/ha) for different species of livestock where i represents livestock species such as dairy, beef and sheep.
Stocking rates may have been different in the past especially when the livestock population was much lower than today. Due to lack of any such information for the past, we use the current standard stocking rates which are 2 (dairy), 3.3 (beef) and 20 (sheep) (Nix, 2003) for the entire study period. Any grass areas left after allocating to different livestock management were assumed to be ungrazed (hay or silage). In locations where the grass area was smaller than that estimated based on the livestock population, the model stocking rate was increased to achieve the observed population.
After 1950, further expansion of agriculture occurred with more of the semi-natural land converted to improved grass and improved grass converted to arable whilst a modest area of arable land became improved grass.

Soil
Soil texture and soil profile depth maps for 5 km × 5 km grid cells required by Roth-CNP were created from the Harmonised World Soil Database (HWSD) (http://webarchive.iiasa.ac.at/Research/LUC/External-World-soil-database/HTML/). Soil organic C, N, P and mineral P pools were initialised with the outputs from semi natural model, N14CP (Davies et al., 2016;Tipping et al., 2012) at the point of their transition to agriculture on 1800 and 1950 (See SI, Fig. S2.1). The N14CP model assumes three SOM pools (fast, slow and passive) to describe SOC, N and P dynamics compared to four active pools within the Roth-CNP.
Total SOC from N14CP was distributed between Roth-CNP's carbon pools for both surface and subsurface layers according to the RothC initialisation as follows: Here TSOC N14CP , SOC fast , SOC slow , SOC passive refer to the total, fast, slow and passive N14CP SOC pools and DPM C , RPM C , BIO C , HUM C represent the carbon redistributed to DPM, RPM, BIO and HUM Roth-C pools (Coleman et al., 1997).
Total organic N and P were redistributed to Roth-CNP pools based on the C/N or C/P ratios of fast, slow and passive pools of N14CP model as follows: The BIO pool of Roth-CNP is largely microbial in nature and is assumed to have a fixed C/N (8.5) and C/P (50) ratios. In a similar way, SOP was also allocated to different Roth-C pools. In this way, fast and slow pools C, N and P from N14CP were allocated to the corresponding pools within the Roth-CNP model, without creating or losing C, N and P.

Hydrology
Hydrological inputs such as AET (mm month −1 ), soil moisture (mm), drainage (mm month −1 ) and runoff (mm month −1 ) on a 5 × 5 km square grid covering the UK were estimated by the hydrology component of the LTLS-IM (Bell et al., in prep). The hydrology model is summarised in Supplementary information (SI 4).

Fertilizer
Manure and fertilizer application rates to arable and grass land were calculated based on the information available from various sources. During the period 1800 to 1840s, sewage in the form of "night soil" was applied to crops and grass (Naden et al., 2016). After 1840, imported N fertilizers (seabird guano, Chilean nitrate) and superphosphate were applied in small amounts. Average N fertilizer input to agricultural land during this period was calculated based on the total fertilizer use (Archer, 1985) and the total area under agriculture (see Section 2.2.3). The average per hectare fertilizer use increased from 7.2 to 13.1 kg ha −1 for N and 4.4 to 16.2 kg ha −1 for P during 1840 to 1940 (Archer, 1985), with 75% of these nutrients were assumed to be applied to arable and 25% to the grass. Chemical fertilizers were applied from 1940s and their rates increased over the years (Archer, 1985;DEFRA, 2011b). For example, N fertilizer application in winter wheat increased from 19 to 195 kg N ha −1 and 4 to 100 kg N ha −1 for grass during 1943 to 2010 (Fig. 3). Mineral N fertilizer application before 1940 was small.

Manure
Manure contribution by deposition of grazing animals (beef, dairy and sheep), slurry, and poultry is calculated based on livestock population and their daily manure (dung and urine) excretion rate.
Carbon and nutrient contributions from deposition of grazing animals depend on the frequency of manure deposition, dry matter (DM) content, carbon, organic-N, NH4-N, and P content of the urine and dung for different livestock species (Table 1). Carbon and nutrient concentrations in dry matter are estimated by the equations where C dep, j, k and eC dep, j, k , are the carbon content (g) and carbon concentration in the dry matter (g event −1 ). F dep, j, k is the frequency of occurrence of dung or urine event (month −1 ) for different animal species. N dep, i, j, k , is the nutrient, i (NH 4 -N, NO 3 -N, organic N, inorganic P and organic P) deposited (g animal −1 month −1 ) in the form of dung or urine (j) of different livestock species (k) and eN dep, i, j, k is the nutrient deposited by urine or dung event (g event −1 ). Slurry is collected when cattle are housed during winter (for dairy and beef). Slurry production depends on the slurry volume, density, DM content, and the nutrient content (Table 1) of livestock species (beef, dairy and pig) as follows: where C sl, i, k is the carbon (g C animal −1 month −1 ) in the slurry of livestock species k. The variables V sl, k , fDM sl, i, k and D sl, k represent the volume (m 3 month −1 ), volume fraction of DM (m3 m −3 ) and density (g m −3 ) of the slurry collected from each animal for a given livestock species k, and cN sl, i, k is the nutrient concentration (g nutrient kg −1 of DM) of the slurry for a given livestock species k. Poultry manure is collected during the whole year and of rate of C (C man, k ) and nutrients (N man, i, k ) produced (g animal −1 month −1 ) is given by depends on the manure DM production (DM man, k , g DM month −1 ) and C (cC man, i, k , g C g −1 DM) and nutrient concentration ( cN man, i, k , g nutrient g −1 DM) ( Table 1). Nutrient input from animal excreta is not directly linked to the nutrient concentration of the grass the animals eat because it is likely that animals also eat food supplements that could affect the nutrient concentration in their excreta. A part of NH 4 -N is lost through volatilization (from manure management practice (housing, manure storage & application to land) and is found to be 0.09 and 0.60 for NH 4 -N in the urine and dung deposition by cattle and sheep, respectively (McGechan and Topp, 2004;Whitehead, 1995). For slurry, volatilization fractions for dairy, beef and sheep are 0.6, 0.31 and 0.6, respectively. For poultry manure the volatilization loss fraction is 0.3.

Historical to current simulation
Prior to 1800, landuse in the UK was largely semi-natural, dominated by natural forests and low input agriculture. Before 1800, a semi-natural model (N14CP, Davies et al., 2016) was used to simulate all landcover including agriculture and to estimate the carbon stocks and nutrient fluxes from 10,000 BP onwards. Following the agricultural revolution (assumed to take place in 1800), Roth-CNP was used to provide a dedicated simulation of agricultural practices in arable and improved grasslands while N14CP continued to simulate nutrient pools and fluxes in semi-natural areas. The soil macronutrient initial condition required by Roth-CNP to allow simulations to start at 1800 was provided by N14CP. A further expansion of agriculture into previously semi-natural areas in the mid-20th Century, assumed to take place in 1950, necessitated a second exchange of soil nutrient pools between N14CP and Roth-CNP in these areas.
Here we aim to estimate C, N and P pools, pool changes, their balance and the nutrient fluxes exported from arable and grassland systems in the UK on a 5 × 5 km grid across the whole of UK during the historical to current period (1800 to 2010). Crop and grass landuse models were run separately on the arable and grassland area within each 5 km grid cell. Based on the five crops in landuse (Section 2.2.3), we identified a maximum of five possible crop rotations with the actual number of rotations dependent on the number of crops in each grid cell (Fig. 4). The number of simulations in each grid cell depends on the number of these rotations and the model variables are re-initialised each time during these simulations. In this way, all the crops that are present in each grid cell are simulated in each year. To calculate the mean yield of a crop we took the weighted average of the yield for each crop for each year by multiplying the area of the crop in each rotation at the end of all the simulations for all the rotations.
Similarly to crop rotations, the model runs for different grazing management systems after re-initialising the model variables for soil and plant growth at each time and the nutrient fluxes are calculated as weighted averages of the area under each grazing management. Under different livestock management systems, animals graze from April to September and the rate of manure (urine and dung) input and the grass removed depends on the stocking rate and animal species (Coleman et al., 2017). During winter when animals are housed, the manure is collected, stored and applied in March in the form of slurry. Nitrogen and P fertilizers are also applied and their rate increases over time, peaking in the late 20th century before starting to decline in more recent years (Fig. 3). All of the P fertilizer is applied in spring whereas N fertilizer is applied in splits (up to 6 in 1990 compared to one single application in 1950 (DEFRA, 2010).
In 1950, widespread landcover change in the UK resulted in different landcover histories depending on the location (grid cell) (Fig. 5). For computational simplicity, Roth-CNP soil variables were reinitialised in 1950 with the outputs from the semi-natural model N14CP (Davies et al., 2016) to incorporate new landcover histories applied from 1950 onwards.
Simulated model results were analysed in three different periods: historical (1800-1950), transition (1950-1970) and current (1970-2010), which are distinct in terms of landuse and agronomic practices. During the historical period, agriculture was more traditional with local varieties and manure and/or slurry based fertilizer inputs. During the transition (post war) period, widespread land cover changes occurred alongside increased use of chemical nitrogen fertilizers in agriculture. The current period is characterised by the so called 'green revolution' effect where improved crop varieties, mechanisation, increased livestock population with higher inputs of chemical fertilizers and pesticides were used, supplementing but also disturbing the natural cycle of C, N and P (Galloway et al., 2004). To calculate the average of a nutrient variable (e.g.: C input, NO 3 -N leached) in each grid cell, we calculated the weighted average for each variable for each year by multiplying with the area under arable or grass land in the cell. The overall C, N and P balance for the whole of UK was calculated by averaging the mean values for these different variables for different time periods across all the grid cells.  (Sakadevan et al., 1993;Whitehead, 1995) NH 4 -N deposited, ( g N event −1 ) 6.07 11.07 0.05 − − (Lantinga et al., 1987;Sakadevan et al., 1993;Whitehead, 1995) Total-P deposited, (  In summary, the changes of SOC (ΔSOC, kg N ha −1 y −1 ), mineral N (ΔN, kg N ha −1 y −1 ) and P (ΔP, kg P ha −1 y −1 ) averaged for the whole of UK are then calculated as where C plant and C animal are the overall mean average annual carbon input through plant and animal sources (Mg C ha −1 y −1 ), CO 2, loss is the loss (Mg C ha −1 y −1 ) of SOC in the form of CO 2 through microbial respiration, DOC loss is the loss of SOC (Mg C ha −1 y −1 ) in the dissolved form through leaching and runoff, and POC loss is the loss through soil erosion in the particulate form. N dep , N min , N BNF , N fert are the overall N inputs through atmospheric deposition, SOM mineralisation/immobilisation, biological N fixation and fertilizer N application (all in kg N ha −1 y −1 ). N loss , N denit and N uptk are loss of nitrogen through leaching, runoff and soil erosion and N removed from soil by plant uptake (all in kg N ha −1 y −1 ). P min and P fert are the overall P inputs through SOM mineralisation and fertilizer P application and P loss and P uptk are loss the of P through leaching, runoff and soil erosion and P removed from soil by plant uptake (all in kg P ha −1 y −1 ). In the model, nutrient (N and P) inputs through litter and those removed by plant uptake are calculated separately. Plant nutrient (N and P) uptake is the cumulative uptake during the whole growing season and a part of these nutrients goes back to soil when plant residues are returned after the harvest (by keeping track of the nutrient concentration in different organs). Litter enters soil SOM pools and undergoes decomposition and forms NH 4 -N by mineralisation and forms NO 3 -N by nitrification.

Historical to current simulation
3.1.1. Historical period  Simulated wheat yields ranged from 0.3 to 1.9 Mg DM ha −1 with an overall mean average yield of 1.0 Mg ha −1 (Fig. 6; Table 2). Simulated potato yields were similar to those of wheat and ranged from 0.1 to 2.0 Mg DM ha −1 with an overall mean average yield of 0.9 Mg ha −1 and simulated fodder maize yield had an overall mean average yield of 4.9 Mg DM ha −1 . For both wheat and potato, simulated yields ranged between 0.08 and 2.0 Mg ha −1 (Fig. 6) lower than that reported for this period in national statistics (Table 2). Simulated grass yields varied widely across the UK from 1.3 to 16 Mg ha −1 (Fig. 6), with the lowest yields occurring mostly in Northern Scotland and Northern Ireland, where SOC was lower than elsewhere. A lower SOC indicates lower SON and SOP and lesser availability of N and P for plant uptake through their mineralisation.
For arable land, simulated average annual SOC change during the historical period is small (−0.08 to 0.12%) (Fig. 7) across whole of the UK with an overall mean net carbon change of −0.18 Mg C ha −1 y −1 (Table 3). During the same period, there was a general build-up of simulated SOC with an overall mean net carbon change of 0.2 Mg C ha −1 y −1 under grass land with a change in carbon ranging from −0.2 to 0.17% annually (Fig. 8). Simulated plant and animal C input to the grass land was greater (2.9 Mg C ha −1 y −1 ) compared to that under arable (1.0 Mg C ha −1 y −1 ). About 93% of simulated total (plant plus animal) carbon input under grass was decomposed, resulting in the build-up of C by 0.2 Mg C ha −1 y −1 .
For arable land, the major part of estimated N input during the historical period was from soil organic matter mineralisation (39 kg N ha −1 y −1 ) Fig. 6. Simulated average wheat, potato and Grass (grazed and/or cut) yields (Mg DM ha −1 ) at different time periods (1800-1950, 1950-1970 and 1970-2010) across the whole UK.
( Table 3). Simulated N was removed from soil mainly by crop offtake (36 kg N ha −1 y −1 ) followed by losses through leaching, surface runoff and soil erosion. Simulated N loss varies across the country with an overall mean average of 15 kg N ha −1 y −1 (Table 3; Fig. 7). However, simulated N loss through denitrification was relatively smaller (0.3 kg N ha −1 y −1 ). For grassland, overall simulated total N input was about 164 kg N ha −1 y −1 with the major contribution from N mineralisation (67 kg N ha −1 y −1 ) and BNF (47 kg N ha −1 y −1 ). Simulated N loss ranged across the country with an overall mean loss of 18 kg N ha −1 y −1 (Fig. 8; Table 3). The net rate of change of N under grass was almost double of that under arable land.
Phosphorus balance takes account of similar components to the N balance except that for atmospheric deposition and BNF (Table 3). Simulated total annual P input includes P from weathering, SOM mineralisation and fertilizer application. Under both arable and grassland, simulated P offtake and P loss through leaching, runoff and soil erosion were less than the P input and resulted in a P build up in the soil at a rate of 2.6 and 2.8 kg P ha −1 y −1 (Figs. 7 and 8; Table 3).
To summarise, estimated crop yields in this period were underestimated when compared to the reported yields between 1880 and 1950. SOM under arable was lost but accumulated under grassland. Mineralisation from SOM was major source of mineral N and P (and BNF in case of N in grass) and most of these nutrients were removed through crop uptake with little remaining for loss through leaching, runoff and/ or soil erosion.  Simulated wheat (0.8 to 4.0 Mg ha −1 ) and potato (1.2 to 5.2 Mg ha −1 ) yields during the transition period were greater than that under historical period with an overall mean average yield of 2.1 Mg ha −1 and 3.4 Mg ha −1 , respectively (Fig. 6, Table 3). Nevertheless, these yields were less than the reported average yield for the whole UK for 1950-1970 (Table 3). Simulated overall mean average fodder maize yield (7.4 Mg ha −1 ) increased by half compared to that during the historical period. Simulated grass yield also increased during this period especially in the western parts of the country (Fig. 6) with overall mean average annual yield increasing by 16% compared to the historical period (Table 3).

Transition period
In arable land, there was a marginal increase in simulated SOC stock particularly in areas of England where grass was converted to arable land in 1950 and elsewhere, SOC changes were less apparent or even decreased (Fig. 7). In grassland, there was a decrease in simulated SOC stock in large parts of England and a marginal increase in the rest of the UK (Fig. 8). Plant derived C was the major source of C input under arable and grassland (Table 3). Under arable land, simulated SOC loss by decomposition exceeds the total C input resulting in decrease in SOC stock during this period. Simulated overall mean average annual changes during this period were − 0.25 and 0.47 Mg C ha −1 y −1 under arable and grasslands, respectively.
For arable land, the major contribution of mineral N in the model is from fertilizer application followed by N input through mineralisation and atmospheric deposition (Table 3). A large part of this nitrogen is taken up by the crop (about 70%) and 25% is lost through leaching, runoff and erosion. Simulated N loss varies across the country with greatest losses occurring in the western England and Northern Ireland (Fig. 7). For grassland, simulated overall average total N input was 206 kg N ha −1 y −1 with the major contribution of N from mineralisation and BNF followed by fertilizer application (Table 3). Simulated N removal by grass offtake is about 85% of this total N and about 11% was lost through leaching, runoff and erosion. Similarly to the arable land, N loss was greater in the western parts of the country (Fig. 8).
Simulated overall average P input (33-34 kg P ha −1 y −1 ) and P offtake (22 kg P ha −1 y −1 ) under both arable and grass were very similar resulting in an overall annual P build up at a rate of 11 kg ha −1 y −1 during this period (Table 3). This period is characterised by an increase in yields compared to the historical period, but yields are still underestimated by the model when compared to the reported values. Soil organic carbon in the model continues to be lost under arable but built up under grassland. Mineral fertilizer became the major source of N and P in arable land whereas under grassland soil mineralisation and BNF continued to be the major sources of N. Although nutrient loss by leaching, runoff and/or erosion has increased during this period, plant nutrient uptake was the major process of nutrient removal from the soil system.

Current period (1970-2010)
Simulated crop yields increased substantially during the current period compared to the transition period (Table 2). Winter wheat yields ranged from 1.4 to 8.7 Mg ha −1 with maximum yields occurring in the South and South east England (Fig. 6). Simulated overall mean average wheat yield more than doubled in the first half of the current period  and increased further during 1990-2010 by another 30% (Table 2). Potato yields increased in some parts of the country to about 8.8 Mg ha −1 with an overall mean average yield of 6.0 Mg ha −1 during this period. Similarly, overall mean average yields for spring barley and oilseed rape has also increased whereas fodder maize yields decreased slightly. However, simulated yields for winter wheat, potato and spring barley were lower by −6%, −25% and −12% than the reported average yields for these crops for whole of the UK during these periods. Simulated grass yields increased especially in the western parts of the UK (Fig. 6; Table 2).
For arable land, simulated SOC decline during the current period was relatively small suggesting that SOC was approaching an equilibrium. The carbon inputs from plant and animal sources were marginally increased and SOC decomposition was slightly less than in the transition period (Fig. 7, Table 3). Under grassland, SOC stock continued to build up during this period with increased carbon input through plant and animal sources with a reduced of loss C through SOC decomposition. (Fig. 8; Table 3).
In both arable and grassland systems, the major contribution of simulated N was from fertilizer (Table 3) during this period with nitrogen offtake by grass more than double of that of crops. Simulated N loss by leaching, runoff and erosion increased and were greatest during this Table 2 Overall mean average simulated crop/grass yields (Mg dry matter ha −1 ) compared to that reported by national statistics a for different time periods in the UK.
The overall mean average annual P fertilizer application increased under arable land (35 kg P ha −1 y −1 ) and decreased slightly under grassland (15 kg P ha −1 y −1 ) during the current period compared to Fig. 7. Simulated soil organic carbon change, average annual N and P losses (leaching + runoff) at different time periods (1800-1950, 1950-1970, and 1990-2010) under arable land for the whole UK. the transition period. As a result, simulated P continued to build up in the soil at a higher rate especially under arable land and there was an increase in overall mean average P loss by leaching, runoff and soil erosion under both arable and grass (Table 3).
Simulated crop yields in the present period were comparable to those in the historical period. This was especially true for winter wheat, the management of which is characterised by high increases in the rates of N fertilizer application (i.e. almost double under general arable and four times grassland compared to the transition period). Simulated N and P uptake and their losses increased during this period with the increases in fertilizer application. Although the rate of P fertilizer application decreased towards the end of this period under grass, the rate of uptake of P and its loss through leaching, runoff and erosion continued to increase because of the large amounts of P in soils that accumulated in the past.

Discussion
The Roth-CNP model developed from Landscape model by simplifying the processes for a monthly time-step reproduced the observed results with varying, but satisfactory degree of goodness of fit for different fertilizer treatments in Broadbalk and Park Grass LTE (See SI 5). In general, simulated wheat yields for Broadbalk and grass yields for Park Grass followed the measured trend although yields were slightly overestimated with an overall RMSE of 1.7 and 1.3 Mg ha −1 , respectively (Table S5.3).
For all the crops, simulated crop yields for the whole UK show a greater yield in the north west of England for different time slices (Fig. 6). This trend is similar to the potential yield of cereals estimated by Sylvester-Bradley and Wiseman (2005) for whole UK, in which this region is characterised by greater summer rainfall and day length compared to the rest of the UK. The terrain and shallow soils may however, limit the actual production in this region. Simulated wheat and potato yields for the whole UK when compared to the national yield statistics reported by DEFRA every year since 1890s show that the model underestimated these yields during the historical and transition periods but agreed well during the current period ( Fig. 9; Table 2). Availability of N for crop uptake is the major yield limiting factor in the model particularly during the historical period. This argument is based on the simulated responses of crop and grass yields to different levels of fertilizer application (please see the supplementary information, S6). Simulated winter wheat yields for Broadbalk averaged 3.5 t ha −1 during the historical period and 7.4 t ha −1 during the current period for a fertilizer application of 144 kg N ha −1 (Plot 8). Whereas simulated grass yields (from Park Grass) were about 7-7.5 t ha −1 during the historical to the current period with a fertilizer application of 96 kg N ha −1 . Assuming a NUE of 50%, 72 kg and 48 kg N ha −1 will be taken up by wheat in Broadbalk plot 8 and by the grass in plot 14 of Park Grass in addition to the N contribution from mineralisation (~50 kg N ha −1 ). This could result in a nutrient uptake of about 122 kg N ha −1 in wheat (plot 8) and 100 kg N in grass (plot 14). Comparing that to the simulated average N uptake at the national scale (Table 3), we can see that N rates (36 kg N ha −1 during historical period) were not enough to achieve the potential yield of 3.5 t ha −1 in wheat yielding only 1 t ha −1 instead, whereas in grass, we can see that simulated yields are closer to the potential of 6 t ha −1 .
Our model uses a simplified approach to development (DVS) as a function of time ( Fig. 2) with detailed processes of water and nutrient dynamics to simulate the crop yields. There is a risk that a fixed DVS based on monthly instead of daily temperature will fail to capture any effect of weather or climate change on the crop duration. On balance, this should not be a problem because the effect of an increase in temperature on crop yield (a decrease of 2.5 to 10% (Hatfield et al., 2011;Lobell et al., 2011) is most likely to be offset by an increase in yield of the same magnitude of up to 11% when CO 2 increased from 280 to 390 ppm during 1850 to 2010 (Long et al., 2005;Long et al., 2006). However, for the future simulations this may be an issue if the negative effect of increases in temperature outweigh the positive effect of an increase in yield due to increase in CO 2 itself and vice versa.
In grass, sufficient N is taken up and a large fraction (1/3rd) of this come through BNF during the historical period. Some BNF is undoubtedly occurring in arable land too, either through leguminous crops and/or free living bacteria (Bohlool et al., 1992;Herridge et al., 2008). However, in the model we did not include either leguminous crops in the rotation or any other form of BNF, which might potentially increase modelled yields more. After 1950, simulated yield increased with the increase in fertilizer application as reported by national yield statistics, but was still underestimated. Biological N fixation might be an additional source of N. Powlson and Jenkinson (1990) reported that BNF could be contributing as much as 25 kg N ha −1 under fertilized treatments in Broadbalk.
An accurate comparison of the simulated and actual grass yield is not possible as the actual grass yield is removed by the livestock species and depends on grazing and grazing intensity in different parts of the country. In general, the simulated grass yields increased over the years since 1800 at an average annual growth rate of 0.6% (Fig. 9). The average grass yields (grazed or cut) estimated by the model (9 Mg ha −1 y −1 ) for the current period show that they are greater than the national average yield (6 Mg ha −1 y −1 ) (Morris et al., 2005). A greater overall mean N uptake, which is more than double that under arable system (Table 3) results in higher yields in the model. For example, during the historical period overall mean N uptake by grass is 144 kg N ha −1 y −1 when N uptake in wheat is only 35 kg N ha −1 leading to a higher yield in the grass (6.4 t ha −1 ) compared to that under wheat (1 t ha −1 ). However, the model does not simulate any effect of changing nutrient, variety or management regimes on grass or crop quality. Table 3 Overall mean average annual soil carbon, nitrogen and phosphorus balance a (for the whole profile) for arable and grass lands estimated based on simulation results for different time periods for whole of the UK.
Components 1800-1950 1950-1970 1970-2010 Arable Simulated SOC for whole of the UK changed in different parts of the country at different rates ( Figs. 7 and 8). Overall, arable lands in the UK have lost SOC during the historical to current period because carbon losses were always higher than the carbon inputs from plant residues compared to improved grasslands which gained SOC during the same period. In grassland the carbon input from plant and animal sources was greater than the loss through decomposition (Table 3). The losses of SOC (and SON and P) via erosion may be underestimated, as it is difficult to capture the effect of extreme, high intensity short duration events in erosion models, particularly at this scale. However, the simulations here provide a broad indication of the scale of C, N and P loss via erosion in comparison to other pathways. Conversely, SOC losses Fig. 8. Simulated soil organic carbon change (%), average annual N and P losses at different time periods (1800-1950, 1950-1970 and 1970-2010) under grass land for the whole UK. may be overestimated in agricultural lands on floodplains, as additions of SOC via sediment deposition during times of flood have not been considered. These are difficult to predict at the spatial and temporal scales of this approach. During the current period , average annual loss of SOC was at a rate of 0.22% for top soil (0-30 cm) (Fig. 9). This is similar to the SOC loss (0.38% y −1 ) reported by Reynolds (Fig. 9), which was higher than that reported by Reynolds et al. (2013 for improved grasslands (0.03% y −1 ). However, Bellamy et al. (2005) reported a higher loss of SOC at a rate of 0.6% y −1 for 15 cm depth for most soils under most landuses in England and Wales. Hopkins et al. (2009) studied the SOC trends under two long-term grassland experiments (that included Park Grass LTE used in this study to test the model) and found that there were no significant trends in SOC as these plots were showing declines, no net changes or increases in SOC. Prior to 1800, a large fraction of the grassland area simulated was under semi natural systems with relatively smaller SOC contents (see SI, Fig. S2.1) led to a build up with change in landuse from semi-natural to improved grass. A peak in SOC under arable in 1950 is due to the assumption of historical-scenario that all the landuse change which happened between 1800 and 1950 occurred in the model in 1950. As a result, a large area of improved grass land has been ploughed up then converted to arable. Similarly, a large fraction of land area changed from semi natural and arable to improved grass in 1950 resulting in a dip initially and a build-up of SOC thereafter (Fig. 9). It is quite possible that when a soil with little SOC is planted to permanent grass, SOC builds up and takes about 100 years to reach an equilibrium (Johnston et al., 2009). However, the effect of the discontinuity in the historical trend around 1950 ( Fig. 9) are short-lived for components apart from SOC. A difference in the structure of the SOC model pools in N14CP and Roth-CNP may also contribute to apparent carbon build up or depletion as a result of placing some fraction of carbon from N14CP model in slow or fast decomposing pools respectively in the Roth-CNP model at the point of landuse transitions in 1800 and 1950. Further work will be needed to investigate the uncertainty in model estimates arising from the simplified historical scenarios we have assumed in this work and distribution of SOM pools between different models at their point of transfer during landuse change.
Under both arable and grass land, mineral N dynamics was dominated by different components at different time periods (Table 3). During the historical period, under arable land, the productivity was mainly determined by soil's natural fertility through N mineralisation and then fertilizers during the transition and current periods. Nitrogen mineralisation depends on SOM content and its rate of decomposition. As a result, a greater N input through mineralisation occurs under grass (65-81 kg N ha −1 ) compared to arable (39-41 kg N ha −1 ). Under grassland, BNF was always a major source of N input to soil during all periods. Legume-based N fixation can vary depending on the grass management and proportion of clover (assumed to be 30% all over in this study). The model estimated overall mean average annual N fixation to vary from 43 to 53 kg N ha −1 for different periods. The quantity of N fixed by high fertilizer, clover-rye grass mixture was 31-72 kg N ha −1 and was less than that of a low fertilizer system (120-160 kg N ha −1 ) (Høgh-Jensen and Schjoerring, 1994). There is always some uncertainty in the rates of natural biological nitrogen fixation (Galloway et al., 2004).
Simulated mineral nutrient pools and fluxes during different time periods are mainly determined by the rates of fertilizer application. Both uptake by plants and losses of nutrients increase with amounts applied. Simulated N loss by leaching, runoff and soil erosion increases through different time periods under both arable (15-52 kg N ha −1 ) and grass (18-36 kg N ha −1 ). These figures for the current period were comparable to those reported by Lord et al. (2002), who estimated N surpluses (i.e. the amount of N that could be potentially lost by leaching, runoff and denitrification) for arable and grassland were 51 kg N ha −1 and 23 kg N ha −1 (after discounting for N removal by grass) in 1995. Overall mean average N loss by denitrification in the model was negligibly small for both arable (0.3-1.5 kg N ha −1 y −1 ) and grassland (0.3-0.6 kg N ha −1 y −1 ) during different time periods. Annual denitrification is variable depending on the N-fertilizer application rate and grazing or slurry application (Whitehead, 1995). Global estimates of denitrification for different combinations soil drainage and N fertilizer application shows 10 and 14 kg N ha −1 y −1 for upland crops and grass for a fertilizer application in the range of 75-150 kg N ha −1 (Hofstra and Bouwman, 2005). A lower denitrification rate occurs in the model because soil is rarely saturated as the soil water is uniformly distributed in the profile as a result of averaging across the whole soil depth. This is a weakness of our approach where soil water is not integrated within the soil model and the total soil moisture storage (mm) estimated by the hydrology model is averaged for the profile depth in the soil model. Although total moisture is same, its distribution within a profile (in different soil layers) may vary depending the season: relatively more water stored at the surface layer during the autumn (rainy season) and at the lower layers during the summer (dry season). Nitrogen offtake estimated by our model for arable (128 kg N ha −1 ) and grass (284 kg N ha −1 ) were higher than that estimated for arable (100 kg N ha −1 ) and grass (116 kg N ha −1 ) land for the whole UK (Lord et al., 2002). Intensively managed grassland, which is harvested by cutting or grazing may yield between 8 and 15 Mg ha −1 y −1 of DM and contain 200-550 kg N ha −1 (Whitehead, 1995). In that case, for an average simulated yield of 9 Mg ha −1 y −1 (Table 2), the grass may well take up N250 kg N ha −1 y −1 . A small loss estimated for denitrification may also contribute to high N uptake in the model especially under grass.
Phosphorus loss varies across the UK with maximum losses found in the North-west England where soils are shallow (Figs. 7 and 8). Similarly to N, overall mean annual P loss through leaching, runoff and soil erosion increased over the years ( Fig. 9) with increase in the P fertilizer application (Fig. 3). Simulated P builds-up in soil during different time periods under both arable and grassland. Withers et al., 2001 estimated the P balance for the whole of UK both under arable and grassland systems for 1993 showing that there was a surplus of 19 and 12 kg P ha −1 . Simulated overall mean average P build up was comparable to that reported by Withers et al. (2001) for arable (18.05 kg P ha −1 y −1 ) but was underestimated for the grassland (3.61 kg P ha −1 y −1 ) ( Table 3). Other studies also found greater P surplus for grasslands in the UK ranging from 14 to 26 kg P ha −1 y −1 from farm to region (CAS, 1978;Brouwer et al. 1995;Smith et al., 1995). The difference in P balance between the simulated and that reported by Withers et al. (2001) is mainly due to high simulated P uptake by grass. However, simulated annual P uptake (30 kg ha −1 ) is similar to that reported elsewhere for grassland systems (Haygarth et al., 1998).

Conclusions
This paper describes an agricultural model (Roth-CNP) that was developed as part of an Integrated Model (LTLS-IM) to simulate the cycles of C, N and P for the whole of UK, comprising terrestrial, hydrological and hydro-chemical model over the long-term period from 1800 to the present. The Roth-CNP model summarises the CNP cycling in an agricultural ecosystem by aggregating soil and crop processes using a daily to monthly timestep. The model simulated crop and grass yields and estimated SOC stocks, DOC and POC losses, and nutrient fluxes (NH 4 -N, NO 3 -N and PO 4 -P) spatially across the whole UK taking into account the biophysical characteristics at each location. The simulated trends of crop yield are comparable to those reported by national agricultural statistics for the same period. Overall, arable land in the UK lost SOC between 1800 and the present day whereas under grassland, SOC stock increased over the same period. This is due to the fact that SOC builds up when a soil with a low initial SOC is planted to permanent grass and it may decrease under arable crops. Simulated N losses were comparable to losses/surpluses reported in the literature. Similarly, P dynamics including P loss and P surpluses were comparable to the literature reports although the P surplus was underestimated for the grass. The model results from the historical to current period show an increase in crop and grass yields especially due to increases in rates of mineral fertilizer application. This has resulted in large positive mineral N and P balances in the soil particularly during the post-war years and the period of the green revolution . Fertilizer inputs largely stabilised or decreased thereafter especially in improved grass and hereon there was increasing impact on nutrient losses rather than on yield. These results clearly show that the model can be used to explore the implications of different management options on crop or grass yields and the nutrient stocks and fluxes at the regional and national scales. The spatial variability of different variables (yield, SOC stock changes and nutrient fluxes) can be attributed to a combination of soil and weather factors. Although crop yields did not show any distinct spatial pattern between different periods of the simulation, in general, simulated crop yields were higher in the central and south west England where there was high rainfall and higher incidence of solar radiation during the growing season. Grass yields were higher in the western parts where rainfall is well distributed throughout the year. Soil organic carbon change did not follow any specific pattern with respect to SOC stock. But the N loss follows the spatial variation of rainfall with a higher loss in the western parts of the UK where rainfall is also higher c). However, no such pattern has been observed for P loss. In summary, the relatively simple agriculture model described in this paper was able to capture variability in the dynamics of CNP at the national scale once coupled to other large scale models of hydrology and soil erosion and driven by atmospheric deposition. The simulation results presented in this study can help farmers to see the effect of agriculture in their region within larger contexts of time and space to understand the consequence of their activities. For a policy maker, this study may help to design policies targeted at the regional level to curb negative impacts of agriculture on the sustainability of the environment. The model could be applied at subnational or catchment scale to attempt to optimise multiple stakeholder interests and for projecting forwards the plausible outcomes under different scenarios of climate and management.