Modelling the response of yields and tissue C : N to changes in atmospheric CO2 and N management in the main wheat regions of western Europe

. Nitrogen (N) is a key element in terrestrial ecosystems as it inﬂuences both plant growth and plant interactions with the atmosphere. Accounting for carbon–nitrogen interactions has been found to alter future projections of the terrestrial carbon (C) cycle substantially. Dynamic vegetation models (DVMs) aim to accurately represent both natural vegetation and managed land, not only from a carbon cycle perspective but increasingly so also for a wider range of processes including crop yields. We present here the extended version of the DVM LPJ-GUESS that accounts for N limitation in crops to account for the effects of N fertilisation on yields and biogeochemical cycling.


Introduction
Nitrogen (N) plays an important role in plant productivity and physiology (Evans, 1989a) and is one of the main limiting nutrients for the functioning of ecosystems in many parts of the world (Grindlay, 1997;Gruber and Galloway, 2008;Vitousek et al., 2002), both in natural and agricultural ecosystems. Historically in agriculture, N limitation for crops has been overcome by the use of manure and N fixing legumes (Vitousek et al., 1997). Since the discovery of the Haber-Bosch process in the 1910s, humans have been able to effectively overcome N limitation by large-scale production and application of reactive N in the form of mineral fertilisers (Vitousek et al., 1997).

S. Olin et al.: Yield and C : N responses to N management and CO 2
The enhanced input of reactive N into agricultural ecosystems by fertiliser use, and deposition to the Earth's surface of nitrous oxides which are by-products from combustion, has together with other technical developments more than doubled global food production during the 20th century (Tilman et al., 2001(Tilman et al., , 2002. However, enhanced N input can also have detrimental effects on biodiversity and water quality, and lead to substantial emissions of N trace gases that affect air quality and climate (Galloway et al., 2004;Rockstrom et al., 2009;Tilman et al., 2002;Vitousek et al., 1997). Better understanding of N effects on yields, conjointly with other ecosystem processes, especially in a changing climate and CO 2 environment is therefore needed for a sustainable management of agricultural ecosystems, weighing enhanced productivity against detrimental side-effects.
N cycling by ecosystems is strongly interlinked with the carbon (C) cycle, which in turn has also undergone drastic changes during the 20th century (Ayres et al., 1994;Rockstrom et al., 2009;Vitousek et al., 1997), as the increased transport of C from the geo-and biosphere to the atmosphere through various human activities leads to an increase in carbon dioxide concentration ([CO 2 ]). Higher [CO 2 ] can have a positive effect on plant productivity -the reason for this is that CO 2 is the main substrate in photosynthesis. Elevated concentrations relative to O 2 in the intercellular spaces of leaves are known to reduce photorespiration resulting from fixation of O 2 by the enzyme Rubisco that catalyses the carboxylation step of photosynthesis (Long, 1991). In addition, enhanced levels of CO 2 result in increased water use efficiency in those plant species that lower stomatal conductance under elevated [CO 2 ], which limits transpirational water loss (Ainsworth and Rogers, 2007;Drake et al., 1997;Sun et al., 2014). However, CO 2 is also a greenhouse gas that leads to higher air temperatures which in turn can either increase or decrease plant productivity depending on the magnitude of the temperature increase. Several studies have assessed the effect that the already experienced environmental change has had on food production (e.g. Lobell et al., 2011;Olesen et al., 2011;Schlenker and Roberts, 2009;Tubiello and Ewert, 2002), and on the projected future changes, using crop models (e.g. Fischer et al., 2001;Rosenzweig and Tubiello, 2007;Rosenzweig et al., 2014). However, the magnitude of the CO 2 fertilisation of crop ecosystems is still under debate (Ainsworth et al., 2008;Rosenzweig et al., 2014;Sun et al., 2014).
A recent model intercomparison highlighted large uncertainties arising from treatment of warming effects vs. effects of CO 2 and N fertilisation on projections of global crop yields . In particular, differences between models, in the representation of a CO 2 fertilisation effect on productivity was highlighted as a key determinant of between-model differences, including globally applicable versions of traditional crop models (e.g. DAYCENT, Stehfest et al., 2007 andGEPIC, Liu et al., 2007), a crop management impact model (PEGASUS, Deryng et al., 2011) and also crop-enabled dynamic vegetation models (DVMs) LPJmL  and LPJ-GUESS (Lindeskog et al., 2013).
For the simulation of crop productivity, traditional crop models typically rely on empirical scaling factors to modify the radiation-use efficiency based on measurements in CO 2 fertilisation experiments (Boote et al., 2013). A mechanistic representation of the CO 2 response (as well as other processes) has been argued to be critical when modelling crop responses to climate change (Yin, 2013), as recently shown for state-of-the-art crop models (Boote et al., 2013). In contrast to crop models, which are optimised to simulate yields, DVMs are tools for exploring and predicting the coupled dynamics of ecosystem functioning, climate-carbon cycle interactions and biome distributions (Friedlingstein et al., 2013;McGuire et al., 2001;Prentice et al., 2007). In DVMs, photosynthesis and stomatal conductance are coupled and respond conjointly to changes in [CO 2 ] (Haxeltine and Prentice, 1996a). New developments in some DVMs in recent years are the inclusion of (1) land-use change and land management functionalities and (2) N cycling (see Arneth et al., 2010;Prentice et al., 2007;Smith et al., 2014;Thornton et al., 2002;Xu-Ri and Prentice, 2008).
The inclusion of N dynamics in DVMs has been found to alter future projections of climate and CO 2 interactions with the C cycle (Arneth et al., 2010;Friedlingstein et al., 2013;Thornton et al., 2009;Wårlind et al., 2014), while the land-use change functionality facilitates assessment of largescale patterns of changes in yields within a consistent model framework that can also address questions such as how management affects the land C sink Lindeskog et al., 2013;Shevliakova et al., 2013). Important management options in this context include decisions on when to sow and harvest, irrigation, residue removal, presence of cover crops, tillage, or fertilisation. For the DVM LPJ-GUESS, including cropland and managed grasslands notably improved phenology when compared with satellite data (Lindeskog et al., 2013). In a study on land-use change in Africa for the 20th century, Lindeskog et al. (2013) found that the impact of implementing land management decisions was of similar importance for the continental C budget as the effect of applying static vs. dynamic land use input data. Levis et al. (2012) showed that including an explicit representation of croplands in the Community Land Model changed both the patterns and amplitudes in the modelled climate compared to treating croplands as unmanaged grasslands. Still, only a few of today's DVMs account for crop processes and C-N coupling in crops (e.g. Arora, 2003;Drewniak et al., 2013), which is a prerequisite to accounting for fertiliser input, the associated effects it has on yields and the C cycle. These improvements will also facilitate global-scale modelling of soil processes such as nitrification and denitrification, because accounting for N uptake through plants will help to constrain ammonium and nitrate amounts, and will hence allow for modelling of soil N 2 O fluxes. While not the focus of this paper, the ultimate goal will be to assess how ecosystem fluxes affecting atmospheric composition and climate vary with changing environmental and socioeconomic conditions.
The standard version of LPJ-GUESS, which simulates potential natural vegetation, has recently been extended to include N dynamics, which has improved the model's ability to represent the biome distributions and productivity patterns globally Wårlind et al., 2014). Here we complement these developments to also encompass cropland dynamics, building on the approach of Lindeskog et al. (2013), including an enhanced temporal resolution of allocation of C and N between different plant compartments. We describe these model developments, and evaluate the model regarding the impact and uncertainty arising from differences in the timing and amount on N fertilisation. We analyse the model's ability to reproduce observed yields on different scales using data from detailed site experiments and regional wheat yield statistics in Europe as a case study. The overall aim of these developments was to find a reasonable level of complexity in processes governing physiology and management for global applications of the model.

LPJ-GUESS
LPJ-GUESS (Smith et al., 2001) is a DVM optimised for regional applications but also applicable globally based on a detailed individual-and patch-level representation of vegetation structure and dynamics. For global applications, vegetation is represented as a mixture of plant functional types (PFTs) that represent the globally most abundant growth strategies of woody and herbaceous vegetation. PFTs are distinguished in terms of growth form, phenology, life history strategy, allometry, photosynthetic pathway (C 3 or C 4 ), climate-dependent scaling of physiological processes and a limited set of bioclimatic limits (Haxeltine and Prentice, 1996b;Hickler et al., 2004;Sitch et al., 2003;Smith et al., 2014). The model uses climate, [CO 2 ], soil information and N deposition as input, and plant communities evolve dynamically through competition in response to these drivers.
Recent model development includes the incorporation of land-use change dynamics together with a crop module (Lindeskog et al., 2013), further developing approaches described in Bondeau et al. (2007) and Waha et al. (2011). In the crop module for global-scale applications, the dominant crop types, such as wheat, maize and rice are represented as crop PFTs, which differ amongst others in management-related parameters such as baseline sowing and harvest dates. Sowing and harvest decisions are modelled based on climate variability (Waha et al., 2011;Lindeskog et al., 2013) and climatic thresholds . Irrigation, residue removal after harvest and cover crops between the main growing-seasons, are further management options available in the model.
The present study builds further on LPJ-GUESS version 3.0 which includes N-cycle dynamics for the simulation of potential natural vegetation Wårlind et al., 2014). Soil C and N dynamics are based on the CEN-TURY model (Parton et al., 1993) which represents 11 soil organic matter (SOM) and litter pools that differ in their C to N ratios (C : N) and decay rates (K d ). Both C : N and K d are dynamic within certain limits -see Smith et al. (2014) for details. SOM decomposition depends on the C : N, K d as well as soil temperature and water content, and may result in either mobilisation or immobilisation of mineral N. Plant N uptake varies between PFTs which differ in their N demand and their competitive strength for N uptake. See Sect. 2.1.2 and Smith et al. (2014) for more details.
Allocation of the net primary productivity (NPP) to different plant organs is done on a yearly basis, based on a set of C allocation rules . If a plant experiences water or N stress during the year, the C allocation scheme is flexibly adjusted so that a larger proportion of the assimilates are distributed to the roots to alleviate these stresses during the following year.
However, for crops, growing periods are less than 1 year and an annual adjustment of the allocation and growth of different plant organs is not sufficient. Lindeskog et al. (2013) partly address this issue, incorporating a C allocation that operates on a daily time step. To allow for dynamic adaptation of the allocation as a response to stress, a more detailed representation of the allocation is developed in this study (see Sect. 2.1.1).
Below we describe and evaluate an updated version of LPJ-GUESS incorporating C-N interaction also for crops. The model allocates daily NPP based on the crop's development phase and allows for an adjustment of the allocation scheme based on the current nutrient and water status of the crop.

Crop development
Upon sowing, the development of a crop plant in LPJ-GUESS starts with a seedling that has an initial carbon mass in leaves and roots. The N content in the seedling is initiated with the highest N concentration ([N]) (the minimum C : N, C : N leaf,min ) allowed in the model assuming a seed with a high N density.

Development stage
In most ecosystem and crop models, plant phenological development is modelled based on weather conditions (Sinclair and Amir, 1992;Sitch et al., 2003;Smith et al., 2001;Wang and Engel, 1998;Zaehle and Friend, 2010), often accumulated over a certain time period such as heat units (HU) (Lindeskog et al., 2013;Bondeau et al., 2007). Here we define development stage (DS, Wang and Engel, 1998) as a number between 0 and 2 where: 0 < DS < 1 is the main vegetative phase, at DS = 1 anthesis occurs and DS > 1 represents the grain filling phase. Compared to the HU implementation currently in the model, the use of DS facilitates a more detailed division of the growing period into the different crop phenological stages (Wang and Engel, 1998). Periods when the plant is more susceptible to heat and nitrogen stress can thus be represented in a more precise manner.
DS at a given point in time (t) is a cumulative function of the maximal development rate d r (day −1 ) which differs between the vegetative phase (DS < 1, d r = d r,veg ) and the reproductive phase (DS > 1, d r = d r,rep ). Following Wang and Engel (1998), DS is also modified using dimensionless scaling factors dependent on temperature (f T ), vernalisation days (f vern ) and photo-period (f phot ): (1)

Daily carbon allocation
For the allocation of the plant's daily assimilates, and their partitioning to the plant organs during the growing-season, we use the established allocation scheme from Penning de Vries et al. (1989). This scheme differs from the one implemented in Bondeau et al. (2007) and Lindeskog et al. (2013) in that the allocation of C to the different organs is related to the daily NPP and to DS, as opposed to a function that meets a predefined target at the end of the growing-season. During the first part of the vegetative phase (DS 0.7 for winter wheat) most of the assimilates are used for root (R) and leaf (L) growth to maximise the uptake of water and nutrients and the absorption of radiation for photosynthesis, followed by a period when more of the assimilated C is allocated to the stem (St). After anthesis (DS > 1 for winter wheat), the grain-filling period starts, during which most assimilates are allocated to the storage organs. During this period, cereal crops reallocate some of their nutrients from the vegetative organs to the grains (Bertheloot et al., 2008).
When a plant experiences water or nutrient deficit during the vegetative phase, it starts to invest a relatively larger fraction of the assimilates into roots to overcome the stress (van Keulen and Seligman, 1989). It is thus important to be able to model the allocation to the roots separately from the other organs. The ratio between the allocation to leaves and stem (L : St), can be treated as constant during stress (Penning de Vries et al., 1989) and thus a relationship between the allocation to R and that to the vegetative parts (V = St + L + R) that is also valid under stress can be established. This approach also gives an opportunity for future implementation of dynamic adjustments in the allocation during the vegetative phase, which is lacking in the original allocation model (Penning de Vries et al., 1989).
Relationships between allocation to L, St, R and grains (Y) from the original allocation model of Penning de Vries et al. (1989) were established and fitted to a logistic growth function, a Richards curve (Richards, 1959), (Eq. 2): where f i is the daily allocation of assimilates to a plant organ relative to e.g. the shoot, a is the asymptote when DS → 0, b is the upper asymptote when DS → ∞, c the growth rate, and d is the DS of maximum growth.

Roots
The allocation to R (g R ) relative to the vegetative organs (g V ) (Eq. 3) is shown in Fig. 1a: (3)

Leaves and stems
Reflecting the shift from L (g L ) to St (g St ) allocation during the initial part of the vegetative phase as outlined above, a relationship between the two organs was derived (Eq. 4) which is illustrated in Fig. 1a: (4)

Harvestable organs, grains
Finally a relationship of the allocation to grains (g Y ) as the fraction of the whole plant (g Y + g V ) allocation (Eq. 5) was derived:

Dynamic allocation
These relationships between the allocation to the different organs of the plant can be applied to favour allocation to one organ over others. Combining Eqs. (3)-(5) yields which is illustrated for winter wheat in Fig. 1b and for spring wheat in Fig. A1b.

Carbohydrate retranslocation
Crops store an easily mobilised reserve of carbohydrates in L, St and R (for some crops also tubers) (van Ittersum et al., 2003;Penning de Vries et al., 1989). To represent this in the model, a labile C pool is filled with a fraction of the daily assimilates directed to the stem (g St ), set here to 0.4 for wheat (Penning de Vries et al., 1989). The labile C pool (M C,labile ) is constrained between 0 and 0.4M C,St . During days when the daily assimilated C is lower than respiration costs (negative NPP), these sugars are used to compensate the loss (Seligman et al., 1975). Additionally, during the grain-filling period the labile C pool is used to add to the grains and is reduced with a rate of 0.1 day −1 (Penning de Vries et al., 1989).

Daily nitrogen allocation
During the vegetative phase in which the leaves and roots are expanding, the plant seeks to maximise photosynthetic gain by having a leaf N content that optimises the carboxylation capacity (V max ) (Hirose and Werger, 1987;Kull, 2002). Following Haxeltine and Prentice (1996a) and Smith et al. (2014) this is done by calculating the V max that maximises canopy-level net C assimilation given the current temperature, water status and biomass C : N.

Leaf N content
Nitrogen associated with Rubisco, the key enzyme in photosynthesis, makes up more than 20 % of the total N in the leaves of wheat (Evans, 1989a), but N is also important for plant structural tissues (Anten et al., 1995;Hirose, 2005;Kull and Jarvis, 1995). However, the vertical distribution of N in the canopy is not even. Higher [N] is usually found in the upper part of the canopy, where leaves experience the highest levels of irradiance (Hirose and Werger, 1987;Hollinger, 1996;Evans, 1989b), compared to the more shaded leaves below. The decline in leaf [N] with the increase in cumulative leaf area index (LAI) from top to bottom typically follows an exponential decrease with a N extinction coefficient k N that is related to the light extinction coefficient (k L ) as follows: where β 0 and β 1 are regression coefficients taken from Yin et al. (2003). From theory on optimal N distribution in a crop canopy, Yin et al. (2000) derived a relationship between the LAI that can be supported given the amount of N that is currently in the leaves (LAI N ) and k N : where M N,L is the leaf N mass and N b is the minimum N requirement for the leaf to function: where C : N L,max reflects the minimum N required for photosynthesis and SLA is the specific leaf area (m 2 kgC −1 ). LAI N is then compared to LAI to determine the N status of the canopy, see Sect. 2.1.3.

Root N content
The N requirement of the root follows that of the leaves through the functional balance concept (van Ittersum et al., 2003;Smith et al., 2014;Zaehle and Friend, 2010): where M N,L denotes leaf N mass, M C,L leaf C mass, M N,R root N mass and M C,R root C mass. The theory behind the concept is that the activity of the roots (uptake and transport S. Olin et al.: Yield and C : N responses to N management and CO 2 of water and nutrients) is proportional to that of the leaves (photosynthesis). A high photosynthesis rate in the leaves (high [N] L ) implies a corresponding relative [N] in the roots to supply the demand of the leaves Zaehle and Friend, 2010).

Plant N uptake
Following Smith et al. (2014), plants take up N from the mineral N pool in the soil on a daily time step as the lesser of the plant demand versus the amount of mineral N in the soil accessible for the plant. N demand from leaves and roots depend on their current C : N status, as the plant seeks to reach optimal C : N in leaves and roots. The mineral N accessible for the plant depends on soil temperature and fine root biomass -see Smith et al. (2014) and Zaehle and Friend (2010) for details. For crops, we have expanded the soil N module so that the N available for uptake by the plant (M N,avail ) is related to the water content of the soil (Eq. 11), as proposed by Xu-Ri and Prentice (2008): where ϕ is the fraction of projected leaf coverage by the plant (proportional to the fine root area), M N,soil is the mineral N mass of the soil and θ is the mean water content of the soil profile.

Senescence
Senescence, the killing of cells, can be either genetically programmed and age dependent, or induced by stresses or environmental factors (Thomas and Stoddart, 1980). In the Conly original cropland version of LPJ-GUESS (Lindeskog et al., 2013), leaf senescence is a function of HU (see Sect. 2.1.1). We develop this further here with a dynamic response of plant senescence to its N status (Yin et al., 2000;Masclaux-Daubresse et al., 2010) and age (DS).

Leaf senescence
If the N status of the leaves is suboptimal, the plant tries to maximise the leaf N in the canopy by redirecting some of it from the shaded leaves towards those that are more sunlit (Hirose and Werger, 1987;Yin et al., 2000). This will eventually turn off the photosynthetic apparatus in the leaves from which all the non-structural N has been retranslocated ( Thomas and Stoddart, 1980). Senescence of part of the canopy in the model is induced when the N-determined leaf area index LAI N (Eq. 8) (Yin et al., 2000) is lower than the actual LAI (Eq. 12): where M C,L is the total C mass of the green leaves. For plants, senescence of leaves is not an instantaneous process. The time for the acclimation of the N content in crop leaves was estimated to be around 10 days (Seligman et al., 1975) which is within the range of what is observed for natural vegetation, 5-30 days (Kull, 2002). Implemented here is the proposed reduction of the leaf C mass as in Yin et al. (2000) m C,sen but with an inertia of 0.1 day −1 : The leaf C mass is then updated M C,L = M C,L − m C,sen and N accordingly using the minimum N content of the leaves, C : N L,max , M N,L = M N,L − m C,sen C : N −1 L,max . The senesced C and N is then transferred to a pool of dead leaves with a high C : N, currently set to 100 (van Ittersum et al., 2003) and the residual N is translocated to the labile N pool. In contrast to the labile C pool, N allocated to the labile pool is not determined as a fraction of the total allocation. The amount is constrained by the N translocated from senesced leaves (Eq. 14) and roots accordingly through the functional balance concept (Eq. 10). The N that is translocated to the labile N pool due to senescence of the leaf is the leftover after maximising the C : N L status: where C : N L,opt is the C : N below which a decrease has a small or no effect on photosynthesis which is estimated here as 3/4 of the range between C : N −1 L,max and C : N −1 L,min . In ageing leaves, observed enzyme efficiency is reduced. After anthesis, degradation of the enzyme Rubisco is higher than the de novo synthesis (Bertheloot et al., 2008). To reflect this in the model, a reduction of the leaf N content at rate of 0.1 day −1 (Bertheloot et al., 2008) starts at anthesis (DS > 1).
In order to avoid excessive allocation of C to the leaves while the plant experiences leaf N deficit (m C,sen > 0) during the vegetative phase, a rescaling of the factor that controls the flow of assimilates to the leaves (f 2 , Eq. 4) was implemented:

Root senescence
Root senescence is still an unexplored area (Kunkle et al., 2009). In the absence of a full mechanistic understanding, the dynamics of the root in the model are assumed to be coupled to those of the leaves through the functional balance concept (Eq. 10).
S. Olin et al.: Yield and C : N responses to N management and CO 2 the model as a transport of N from the leaves, roots and the labile N pool. In the model the plant tries to meet the demand from the grain: primarily by reducing the labile N pool, M N,labile .

Nitrogen retranslocation
If m N,Y,dem is larger than the labile N pool, the crop plant attempts to meet the unsatisfied N demand from the grains (m N,Y,dem = m N,Y,dem − M N,labile ) by a N transport from the donor organs (leaves and roots). These donor organs have a resistance to let go of their N, r (Seligman et al., 1975), to account for the fact that N is needed for maintaining organ processes (e.g. photosynthesis and maintenance respiration): where j denotes the organ, L or R. The actual transport of N (m N,retr ) is calculated by summing the individual organs' relative portion of the total N demand from the grains after the labile pool has been emptied (Eq. 18). If the demand on the organ is larger than the available N, it is reduced to its minimum N content (C : N j,max ): During the initial part of the grain filling period, only leaves contribute to fulfilling the grain N demand. Once more than half of the assimilates goes to the grain (DS > 1.15, see Eq. 5), the model can utilise part of the plant root N as well to fulfil the N requirements of the grains.

Updated soil water parameters
Soils are characterised by their ability to store and provide water to the plants; a parameterisation of these soil water characteristics based on fractions of grain sizes, available for the soils in the study area, was needed for this study. Soil water characteristics as used in LPJ-GUESS were derived from data on sand, silt and clay for the top soil layer taken from a map of soil mineral fractions. These fractions were then used as input to empirical relationships (Cosby et al., 1984, Table 3) for the following soil water characteristics: soil water pressure at saturation ( s ), volumetric water content at saturation (θ s ) and a shape parameter describing the response of the water retention curve to changes in water content (b). These parameters were then used in Eq. (19) to derive the volumetric water content under specific conditions: where i is the actual pressure head (m) and θ i is the actual volumetric water content (m 3 m −3 ). The percolation coefficient K (Haxeltine and Prentice, 1996b;Gerten et al., 2004), an empirical parameter used in the model to derive the daily percolated water, was fitted against b values for four of the soil classes from Haxeltine and Prentice (1996b) (coarse, medium-coarse, medium, fine) and resulted in 3 Experimental setups The model's ability to simulate yields was evaluated using data from fertiliser trials from the Netherlands, a Free Air CO 2 Enrichment (FACE) experiment from Germany and regional yield statistics from European countries. All simulations were performed using a 500-year spinup using [CO 2 ] and repeatedly cycled detrended climate input for the first years of the historic climate data to build up pools of C and N. In the simulations with N dynamics turned on, monthly N deposition input from Lamarque et al. (2010) that varies decadally was used. The values were interpolated using bilinear interpolation from the original resolution (1.9 • × 2.5 • ) to match the resolution of the climate data (0.5 • × 0.5 • ) Wårlind et al., 2014).

Fertiliser trials
To evaluate the model's ability to simulate phenology and yields, and sensitivity to N fertiliser additions, data from nitrogen fertiliser response trials with detailed measurements of dry mass and N mass allocation from the Netherlands (Groot and Verberne, 1991) were used. In the trials, winter wheat was grown with different fertiliser input in the years 1983-1985. The trials were conducted on three sites (The Eest, The Bouwing and PAGV), located in the central part of the Netherlands, see Fig. A2. At these locations, three different N treatments were carried out for two seasons (Table 1). Based on an initial calibration of the model using leaf phenology data from Trial I, the parameters a and b in the allocation function f 2 (Eq. 4) were changed from 0.88 and 0.09 to 0.8 and 0.2 respectively, for this application. Daily climate data for the years 1979 to 1984 were downloaded from the Haarweg weather station, Wageningen University 1 , located within 70 km from the sites. To initialise (spin up) N and C pools in the model, climate data for the year 1979 were repeated for 500 years. In Groot and Verberne (1991), there is no information on management practices in previous years, so we decided to implement a moderate level of 100 kg N ha −1 y −1 as a single application, 150 days after sowing for the spin-up. The year before the trials started Table 1. Site-and treatment-specific data after Groot and Verberne (1991). For all trials (I-VI), three experiments with different applications of N fertiliser were performed (1, 2 and 3). Their timing is expressed here by the development stage (DS).

FACE experiment
The ability of the model to simulate the observed response to elevated [CO 2 ] treatment on yields and C : N of cropland ecosystems was tested using a FACE experiment (Weigel and Manderscheid, 2012) Lindeskog et al., 2013) represented by wheat (spring and winter) in the model, therefore growth of cereals was simulated for all years. Four different trials from the experiment were simulated, high (100 % N) and low (50 % N) N input with ambient and elevated [CO 2 ] (378/548 ppm), see Table 2 or Table 1 in Weigel and Manderscheid (2012) for more details. Due to the lack of information on the timing and amount of the individual fertiliser applications, these parameters were set using the results from the regional comparison (Sect. 3.3); total amount of N added in the experiments are listed in Table 2. As climate input we used the WFDEI climate data set (Weedon et al., 2012) which is a bias-corrected reanalysis data set based on WATCH (Weedon et al., 2011) and Era Interim (Dee et al., 2011). During the spin-up period (500 years), 30 years of detrended data  were used together with the ambient [CO 2 ] from 1979.

Regional yields
To evaluate the model's ability to simulate wheat yields on a larger scale, 65 regions at the NUTS2 (Nomenclature of Territorial Units for Statistics in the EU; statistical administrative areas) level from northwestern mainland Europe and NUTS1 for southern England were selected based on their cereal fractions and yields taken from the EU statistics website, EUROSTAT 2 . The regions and their administrative names are shown in Fig. A2 and the mean reported yields together with the number of years for which there were data for each region are listed in Table A3.
As climate input, WFDEI (Weedon et al., 2012) with a spatial resolution of 0.5 • × 0.5 • was used. Spin-up with climate and [CO 2 ] was performed as in the FACE experiment (see above). Soil characteristics (mineral size distributions) for the top layer (0.3 m) were derived from Batjes (2005) and used for the whole column (2 m). Simulated yields of spring and winter wheat were assigned to the regions by their relative proportion of the grid-cell area (a grid-cell can belong to more than one region). Fractions of each grid-cell covered by spring and winter wheat as well as the area equipped for irrigation were derived from the MIRCA data set (Portmann et al., 2010).
Timing of the fertiliser applications were selected based on the mean development stages for the three N applications listed in Table 1 (0.18, 0.49, 0.89). A common practice is to apply some or all of the fertiliser at the time of sowing (Mahler et al., 1994). The timing of the three applications was therefore changed to DS = 0, 0.5 and 0.9 respectively.
To test the effect of timing and amount of fertiliser applied, an experiment with 50 model permutations (F T,A , timing (T) and application (A) varied), was conducted with five different fertiliser application rates (between 50 and 250 kg N ha −1 y −1 ). Fractions (0, 1/3, 2/3, 1) of the applied N were distributed at the development stages 0, 0.5 and 0.9, yielding 10 possible combinations.
The N managements (application rate and timing) that gave the best fit (lowest RMSE) were then selected for each region (F opt(T,A) ). To test whether a mean N management can be representative for the whole region, the mean timing and amount from F opt(T,A) were derived and simulated (F t,a , timing (t) and application (a) fixed). Additionally, to test the relative importance of timing and amount, an experiment with simulations where timing was fixed using the mean development stages as in F t,a together with varying input as in F T,A (F t,A ) and an experiment where timing was varied as in F T,A but with a fixed N application, were performed (F T,a ). Furthermore a simulation using a gridded data set of annual N input for wheat  together with the mean timing from F t,a was performed (F t,I ).
To test whether adding C-N dynamics in the model increased the overall model performance, simulations using the C-only version of LPJ-GUESS were performed (F C ). In these simulations, C allocation was as described in Lindeskog et al. (2013).
A short description of the model setups used, together with their abbreviations, are listed in Table 3.

Statistical methods
In order to quantify the degree of agreement between simulations and the associated observations, two indices were calculated, the Willmott index of agreement (W i , Eq. 21) (Willmott et al., 2012) and the root mean square error (RMSE) (Eq. 22). W i is calculated as where M denotes the sum of absolute differences between the modelled and the observed mean, O is the sum of absolute differences between the observations and the observed mean and c is scaling constant here set to 2 (Willmott et al., 2012). W i is without unit and ranges from 1 to −1, where 1 is a perfect agreement between the modelled and observed variances and −1 means that there is a low or no agreement between the modelled and observed variances. RMSE is calculated as where n is the number of observations, o is the observed value and m is the modelled value.

Conversion factors
To convert plant C to total dry matter, a conversion factor of 0.446 was used (Osaki et al., 1992). Dry weight was converted to wet weight (used in the regional statistics) by assuming a wet fraction of 0.15 in the grains (Fader et al., 2010).

N fertiliser response
In the N fertiliser experiments from Groot and Verberne (1991), 18 trials with N input ranging from 0 to 240 kg N ha −1 y −1 applied at various crop development stages were performed, resulting in grain C production from 1 to 3.5 ton C ha −1 . During the growing-season, leaf C in the field trials increased until peaking around June, after which senescence commenced and leaf C decreased again (Fig. 2ac). Simulations with LPJ-GUESS at these sites, and reproducing the applied fertiliser scheme broadly captured these seasonal dynamics, and the response to the different levels of N applications (Fig. 2a-c). Modelled grain and above-ground C mass per kg N applied (19 and 46 kg C kg N −1 ) were in line with the observed response of 22 and 42 kg C kg N −1 , which indicates appropriate sensitivity of yield and growth to nitrogen addition. Differences between simulated and observed leaf C values were largest towards the end of the growingseason (Fig. 2d), especially at the highest-fertilised trial sites and in the second growing-season. As seen from the example time series in Fig. 2a-c, rates of senescence in the simulations were too slow, compared to measurements, which resulted also in underestimated dead-leaf C (see Fig. A3a-b). The model generally simulates the observed above-ground production of biomass (C and N) well at the sites, with more accuracy in the medium-and high-input trials (2 and 3) (Fig. 3c-d). C content in grains and above-ground biomass is captured reasonably well, with some underestimations for the lowest N trial (Fig. 3c). This picture was mostly similar also for simulated N content, as a consequence average C : N were not biased towards too high or too low values (Fig. 4a). The C : N of the grains (Fig. 4b) in response to dif- Table 3. Description of the setups used in the comparison with regional statistics, with abbreviations used throughout the paper.

Setup
Description Timing N app.
F T,A 50 permutations, with timing and application rate varied varied varied F opt(T,A) optimised timing and rate based on model fit (RMSE) in comparison to regional yield statistics opt.
opt. F t,a mean timing and application rate over all regions from F opt(T,A) fixed fixed F T,a 10 permutations with timing varied as in F T,A , application rate from F opt(T,A) varied fixed F t,A 5 permutations with timing from F opt(T,A) , application rate from F T,A fixed varied F t,I Timing from F opt(T,A) , application rate from  fixed input F C C only The difference between observed and simulated leaf C for three different levels of fertiliser application for the Bouwing, the Eest and PAGV (Netherlands) for seasons 1982-1983and 1983-1984(Groot and Verberne, 1991. Blue symbols indicate lowest levels of fertilisation; red represent medium and black symbols a high N fertiliser input. Open symbols are for the season 1982-1983, and closed symbols are for the season 1983-1984. ferences in N treatment was better captured than the total C mass (Fig. 3c). Over the growing-season there is an underestimation of the N content in the grains, especially so for the low input treatments (Fig. 3a-b).

Response to elevated [CO 2 ]
At the Braunschweig FACE experimental site the mean observed yields under ambient [CO 2 ] (∼ 378 ppm) were 8 and 6 ton ha −1 for the sufficiently fertilised (100 % N) and the treatment receiving 50 % N respectively, whereas the simulated yields were 9 and 7 ton ha −1 , with the same pattern of higher simulated than observed yields also for the elevated [CO 2 ] (∼ 548 ppm) treatments (Fig. 5). Grain yields were simulated to rise by 19 % as a response to elevated [CO 2 ] for both the 100 and 50 % N treatments (Fig. 5). The observations show a similar response with a rise of 14 % (9-19 % for 100 % N and 5-24 % for 50 % N), and neither simulations nor measurements indicated a clear impact of N treatment on the CO 2 effect.
Under elevated [CO 2 ], increased C sequestration and yields were not balanced by grain N rising at the same rate as grain C, leading to enhanced grain C : N at elevated CO 2 . In the observations, this increase was on average 16 % for both N treatments, whereas in the simulations the increase was 24 (100 % N) and 20 % (50 % N) (Table 4).

Regional wheat yields
In order to test model performance on spatial scales beyond field trials, regional modelled wheat yields were compared with yield statistics provided by EUROSTAT. The simulations were also set up to test the effects of different N management regimes (Table 3) in order to derive an implementation of fertiliser application that can be adopted for large-  N (b), and total biomass C (c) and N (d) for the Eest, the Bouwing and PAGV, the Netherlands over the growing-seasons 1982-1983 and 1983-1984. Symbols are the same as in Fig. 2 with blue symbols for low input of N fertiliser, red for a medium input of N and black for high N input.

Figure 4.
A comparison between modelled and observed C : N in harvested above-ground biomass (a) and grains (b), for the Eest, the Bouwing and PAGV, the Netherlands for the seasons 1982-1983 and 1983-1984. Symbols are the same as in Fig. 2 with blue symbols for low input of N fertiliser, red for a medium input of N and black for high N input.

Figure 5.
Effect of CO 2 fertilisation on observed and simulated grain yield, comparing wheat grain yields grown at elevated CO 2 (∼ 548 ppm) with those grown at ambient CO 2 (∼ 378 ppm). Simulated yields are depicted by solid lines and filled circles, observations are depicted by dashed lines and markers, shown for treatments with sufficient N fertiliser input (100 % N, blue), and treatments that received half of that amount (50 % N, red). Observations are from (Table 4 Weigel and Manderscheid, 2012).
scale models even when exact information on fertiliser timing is not available.
To do so, from the 65 regions chosen from the EUROSTAT database, a set of 50 permutations of timing and amount of N application (F T,A ) was used to identify the management strategy resulting in the best agreement (lowest RMSE) with reported time series of yields (F opt(T,A) ). The optimised timing and associated amount of N input for each of the regions are listed in Table A3. This "optimised" simulation deviated only marginally (Fig. 6b) from the observed yields (Fig. 6a). The interannual variability in yields for individual regions was captured best for regions with a low productivity (model performance based on the Willmott index (W i ) and RMSE, Table A3), whereas the interannual variability for the more productive regions was captured less well.
In the optimised simulation (F opt(T,A) ), the mean N application rates for spring-sown wheat across the entire region were a total of 129 kg N ha −1 y −1 , applied in fractions of 0.11, 0.50 and 0.39 at the three development stages described in Sect. 3.3, with the main application in mid-spring. For winter wheat, on average 172 kg N ha −1 y −1 were applied with fractions 0.08, 0.19 and 0.73 for the three development stages, with the main application in late spring or early summer.
For simulation F t,I , the mean timing obtained from the optimisation was combined with a gridded data set of N application rates , resulting in a reasonable agreement with the observed yield but with larger spread compared to F opt(T,A) (Fig. 6c). In particular, overestima-tions were found in parts of the Netherlands, Belgium and southwestern France, and a considerable underestimation in northern France. Some of the deviations between modelled and reported yields were likely due to a lack of spatial variability in the fertiliser data input, e.g. a constant value (110 kg N ha −1 y −1 ) was applied for all regions in France (Table A3).
Despite the spread between the model and the observations for individual years, the temporal average of the optimised set showed a good agreement (Fig. 7). Even the simulation that applied the mean timing with reported fertiliser rates (F t,I ) showed a reasonable agreement, but generally overestimating the yields in low-productive areas, primarily because of a higher fertiliser application rate (regional mean N input of 188 kg N ha −1 y −1 for F t,I , compared with 169 kg N ha −1 y −1 for F opt(T,A) ). Applying the same timing as in F t,I with the lower constant rate of 169 kg N ha −1 y −1 (F t,a ) resulted in a much better agreement for the mean response. However, both simulations (F t,I and F t,a ) have a smaller range of simulated yields from low-productive to high-productive areas than reported in the statistics data.
The spatial variation of the observations was captured well in all simulations except for those with maximised yields (max(F T,A )) and the C-only version (F C ). As expected, optimised N management (F opt(T,A) and F t,A ) improved the fit of the model results to spatial variation in the data, but all the C-N enabled simulations except for max(F T,A ) increased the agreement between the modelled and observed variance (Table 5), when compared to F C . In order to address the spatial variability in timing and rates of fertiliser application, the optimised simulation (F opt(T,A) , red line in Fig. 7) was compared with two additional optimisations. In the first of these (F opt(t,A) ), timing was prescribed using the same as for F t,a while application rates were varying. In the second optimisation (F opt(T,a) ), application rates were prescribed as for F t,a while timing was varying. The grid-cell average yield over the region and all permutations in the F opt(T,A) simulations was 5.2 t ha −1 y −1 , ranging from 2.4 to 10.3 t ha −1 y −1 between the different application rates and timing. For F opt(T,a) the same measures were 5.5 t ha −1 y −1 (3.1-8.7 t ha −1 y −1 ) and for F opt(t,A) , 5.2 t ha −1 y −1 (3.2-8.6 t ha −1 y −1 ). The average yields for all simulations were of the same order of magnitude. For F opt(t,A) and F opt(T,a) the ranges in yield were also of similar size whereas the range for the F opt(T,A) was larger although smaller than the sum of the ranges of F opt(t,A) and F opt(T,a) . Most importantly, both the optimisations with either fixed timing or application rate, resulted in a better agreement with the reported yields than when only using a mean uniform N management over the region (F t,a , Table 5), but optimising the application rates gave a considerably better fit than optimising the timing. While timing had a large effect, these results imply that highest priority is to obtain data on application rates. Table 4. Comparison of modelled and observed grain C : N from a FACE experiment where wheat was grown in ambient CO 2 (∼ 378 ppm) and elevated CO 2 (∼ 548 ppm). The observed C : N were compiled using Tables 4 and 5 in Weigel and Manderscheid (2012), observed C values were derived from dry matter using the conversion described in Sect. 3.5. The range of the observed CO 2 effect is estimated from the standard errors listed in Tables 4 and 5 in Weigel and Manderscheid (2012

Discussion
Accounting for C-N dynamics in the crop version of LPJ-GUESS (Lindeskog et al., 2013) together with the new flexible allocation scheme resulted in good overall agreement when compared against site-scale observations and regional yields statistics. The simulated response to N management was also in line with the observed dynamic responses in a fertiliser trial and wheat grown under elevated [CO 2 ].

Model performance, fertiliser trials
Modelling the seasonality of growth, phenology, and the response to fertiliser is a prerequisite not only for model projections of crop yield responses to management in a changing environment, but also to aid assessments of management-  7. Reported vs. modelled yields for the 65 regions used in this study. Grey triangles represent all available years for each region, with the fertiliser management that gave the best agreement with data for each region (F opt(T,A) , all); the red markers are their means. Dark green markers show the mean for each region with the C-only version of the model (F C ), pale green represents the management that gave the maximum yields for each region (max(F T,A )). Orange markers show the result from using a mean N management over the region, blue represents simulations using the same timing for each region but with a spatially explicit data set of N application  (F t,I ). Lines are fitted linear regressions, see Table 5 for details. The number of years included for each region is listed in Table A3. related detrimental effects (Hawkesford, 2014;Lokupitiya et al., 2009). Models that operate on regional to global scales are not designed to be suitable tools for detailed local-scale decision making on fertiliser use, but the improved representations of crop C-N coupling and phenology are necessary for simulating regional to global land-use-related surfaceatmosphere exchange fluxes, to evaluate models against observations, and to contribute to analyses of the effects of land-use change in the climate system, including assessment of how multiple ecosystem services are affected following land conversions (e.g. Kelley et al., 2013;Rounsevell et al., 2014). For these types of questions, a chief challenge remains in representing phenology and growth responses to fertiliser application and climate in a way that is suitable for large-scale models, but still reproduces realistic results (Boote et al., 2013;Lokupitiya et al., 2009). In LPJ-GUESS, crop phenological events, and especially the growth of the leaves throughout the growing-season were captured (Fig. 2) and at the end of the growing-season, senescence was induced as a result of N retranslocation to grains from the leaves, albeit with a response that was a little too weak compared to the measurements (Fig. 2). The modelled mean harvest index (HI = Y / (Y + V)) of 0.57 was in line with the value obtained from observations (0.52; Fig. 3a,c). Modelled grain yields and above-ground biomass were on average only slightly lower than the observations, and the overall tissue as well as grain C : N also agreed well with the corresponding ratios derived from the measurements. But there was a discrepancy in the modelled C allocated to the grains during the early parts of the grain filling period that disappeared towards the end of the growing period (Fig. 3a).
When making these comparisons, the need to convert dry matter to mass of C, and vice versa, added a level of uncertainty that is not associated with modelled processes, since site data often are reported as biomass (dry or wet). For instance, we used the published observations of [C] in the biomass of 44.6 % (Osaki et al., 1992) when converting the experimental site data. By contrast, in LPJ-GUESS a C content of 50 % is assumed , which leads to slightly higher C density in modelled output compared to the observation-derived values.
Underlying the good agreement of the tissue C : N, was an underestimation of the absolute levels of both C and N in the grains at the end of the growing-season, which were underestimated by about 30 % in many cases (Fig. 3a). The same pattern could also be seen for the above-ground C mass (Fig. 3c), suggesting that the productivity generally is too low in the model. One potential explanation for the underestimation of [C] in grains at the end of the growing period is C-retranslocation from leaves and roots to the grains during the growing-season, which has been observed in forms of sugars (Osaki et al., 1991;Penning de Vries et al., 1989). Since in the current version of the model C supply originates solely from NPP, one potential source of carbon to grain filling therefore is missing. In addition, the model does not presently include other forms of molecular transport whereas the retranslocation of N from leaves to the grains is mostly in the form of amino-acids (van Keulen and Seligman, 1989) with a relatively low C : N (e.g. 4.0 in wheat; Thornley, 1990;Osaki et al., 1991). Accounting for such an amino-acid transport would suggest that for every unit of N that is transported to the grain, four C would have to be supplied as well. This, too, cannot be captured by the current implementation that is based on allocation of the daily NPP.
Since the C : N of the various organs plays an important role for determining N uptake and its N (re)allocation between organs, variable C : N have been suggested to be crucial when modelling vegetation C-N dynamics. This aspect might become especially important under changing climate and [CO 2 ] environments (Yin and Struik, 2010) as these can affect the chemical composition of the plant (e.g. Ainsworth and Long, 2005;Yin, 2002), discussed in Sect. 5.2. For present-day conditions, the simulated biomass C : N agreed best with observations from the two high-N-input treatments (2 and 3) (Fig. 4), although the overall agreement across treatments was acceptable, keeping in mind the assumptions that had to be made when deriving the observation-based values. Likewise, growing-season green leaf [N], which together with LAI is essential for predicting photosynthesis (Boote et al., 2013), was also well reproduced (results not shown). By contrast, [C] as well as [N] in the dead leaf pool were either over-or underestimated (see Fig. A3), most likely because the C : N in dead leaves is set to a constant value.
Tissue C : N in dead leaves will affect C and N input to the soil through litter decomposition, and hence is arguably an important feature for representing the N cycling in ecosystems. We concentrate here on above-ground processes such as live-tissue growth, element concentrations and yields and it remains to be tested whether an implementation of C : N limits that vary over the course of the growing-season (Asseng et al., 2004) will be a necessary improvement for simulating agricultural soil processes. Moreover, a dynamic adjustment of the C and N allocation to root growth under stress (see Sect. 2.1.1) might be an important future development, since as a result, less C would be partitioned towards aboveground growth, while at the same time, more N could be extracted from the soil. For the observations used here for model evaluation, no data are available on the N or C content of the roots. A flexible adjustment of the root : shoot growth is already implemented in the model as a response to N and water stress for natural vegetation, but it operates so far only on a yearly basis . Studies elsewhere have demonstrated dynamically adjusted shoot : root growth and C : N (Poorter et al., 2012), and some published crop models like GECROS (Yin and Van Laar, 2005) do so at higher temporal resolution. Next development phases in LPJ-GUESS will explore how to include these into the model.

Model performance, FACE comparison
In a recent intercomparison of crop models that can be applied globally, Rosenzweig et al. (2014) identified the CO 2 fertilisation response to be one of the major sources of uncertainty of how yields might change in a future environment. When disregarding advances in breeding, the CO 2 fertilisation response of crops (in particular: of the C 3 photosynthesis-type) is fundamental, as this response can counteract or at least dampen effects of climate change. The chief principles are related to the carboxylation reaction of Rubisco being stimulated by enhanced levels of CO 2 , and by plants being able to operate at lower levels of stomatal conductance, thereby increasing the efficiency of gaining carbon per unit of water lost (Ainsworth and Rogers, 2007;Drake et al., 1997;Sun et al., 2014). The magnitude of the positive CO 2 response on yields will depend on the degree of leaf-level acclimation response, and how such an acclimation would translate to the whole-plant level (Ainsworth and Long, 2005;Drake et al., 1997;Smith and Dukes, 2013).
In existing crop models, the implemented CO 2 response is typically based on empirical relationships between an increase in [CO 2 ] and plant productivity, derived from e.g. FACE experiments (Boote et al., 2013;Rosenzweig et al., 2014). Adapting a more mechanistic representation of the CO 2 response by replacing a radiation-use efficiency formulation with a coupled photosynthesis and stomatal conductance description (as in LPJ-GUESS) increased the overall model performance of the DSSAT model and its ability to capture responses to the combined effects of different factors affecting plant productivity (Boote et al., 2013).
In the Rosenzweig et al. (2014) study, LPJ-GUESS showed a very strong yield response to enhanced levels of atmospheric CO 2 . The model version used in that study did not include C-N interactions, and hence this strong response was expected, since the simulated underlying physiology was not constrained by N availability. With this constraint in place, LPJ-GUESS in the present study reproduced crop growth and productivity under elevated [CO 2 ] with different N fertiliser treatments, as well as responses in dry matter C : N observed at the different levels of CO 2 , with C : N at the 100 % N treat-ment being somewhat above the observed values. A positive CO 2 effect on yields and C : N, as well as the higher variability between years was captured in treatments that received less N input. The results shown in Sect. 4.2 were also in line with observations from previous FACE experiments (Ainsworth and Long, 2005) that showed a mean increase of 14 % (−2 to 33 %) for wheat yields, compared to the modelled response of 20 and 24 % for the two different N treatments.
The modelled yields (ambient and elevated [CO 2 ], 50 and 100 % N) exceeded observations only slightly. This might be explained by insufficient information on the fertiliser management at the experimental plots but also by, for instance, a too-high allocation to grains compared to other tissues. Moreover, under elevated [CO 2 ] a relative increase in root biomass (Pritchard and Rogers, 2000) often occurs. Currently, as discussed in the previous section, there are no explicit mechanisms implemented that would yield such a dynamic growth response. A step forward, that remains to be tested, would be to set a flexible root : shoot allocation via a modification of f 1 in response to water or N stress, see Sect. 2.1.1, which would result in such a response in situations when additional leaf C mass in response to elevated [CO 2 ] induced a N demand that cannot be met by soil uptake.
Moreover, as mentioned in Sect. 5.1, SLA in the model is treated as a constant, even though it has been observed to vary over the growing-season, and in response to elevated [CO 2 ] (Poorter et al., 2009). Ainsworth and Long (2005) found in a review of CO 2 effects on different plant traits that elevated [CO 2 ] was associated with a reduction in SLA, as the increase in leaf mass was not accompanied by a proportional increase in leaf area. These results are in line with the observed increases in leaf C : N under elevated CO 2 (Yin, 2002). Currently in the model the lower limit of leaf N (N b ) is a function of SLA (Eq. 9), and whether or not variable SLA would help to constrain the C : N response at variable N and CO 2 treatment remains to be tested. While some studies have found elevated CO 2 to also change the chemical composition of the plants and thus the C : N (Yin, 2002), others (e.g. Gifford et al., 2000) found no evidence of changes in the C : N of the senesced leaves grown under higher [CO 2 ] implying that the fixed C : N limits and in turn the N b are valid also under elevated [CO 2 ], but this remains to be explored.

Regional yields, model performance and implications for large-scale modelling
The implementation of C-N dynamics improved the ability of LPJ-GUESS to simulate yields not only at local scale but also across larger regions, especially when all permutations of N managements were combined (Fig. 7). As expected, the comparison of the simulated yields with reported ones was best captured when considering time-averaged values (accounting only for the spatial variation) compared to the full temporal variability (difference between "mean" and "all" in Table 5). The discrepancy between results from the C-N and the C-only versions of the model was striking and clearly demonstrates the need to consider C-N interactions when modelling crop processes . In addition to the improved phenology and C-N coupling, which was seen also at the FACE site, on the regional scale the representation of soil texture proved to be an additional important aspect. By incorporating the WISE map of soil mineral fractions, we were able to increase the model performance from a W i of 0.715 using the standard soil map in LPJ-GUESS simulations with fixed mineral fractions , to 0.807. With the more detailed soil information, the heterogeneity in the growth response from fertiliser applications due to differences in physical properties could be better captured. Still, when applied across large spatial domains, there can be multiple reasons for a disagreement between modelled and observed variability that go beyond process representation linked to the basic physiology of C-N interactions. For instance, extreme heat or freezing, pests, or water logging of soils are frequent events detrimental for crop production (Reichstein et al., 2013). Effects of extreme weather events are difficult to account for, partially because of the smoothing effect of a daily time step, and also because aggregation averaging in the production of the gridded climate input data tends to remove weather extremes. Likewise, local management decisions that are not based on weather variability, and their effects on crops and environment, are difficult to capture with the current setup.
Several approaches to modelling N limitations in agricultural ecosystems over large region scales are available in the literature. CLM (Drewniak et al., 2013) e.g. includes N limitation for crops, and this model also simulates the retranslocation of N during the grain-filling period based on prescribed C : N of the plant organs pre-and post-anthesis. The C allocation scheme implemented in CLM has the same origin (Penning de Vries et al., 1989) as the one implemented here for LPJ-GUESS. In LPJmL (Fader et al., 2010), an implicit nutrient limitation on yields is considered by varying production parameters (LAI max , HI, and α a ) to match country or region statistics from e.g. the FAO (United Nations Food and Agriculture Organisation). This approach has the advantage that it can be applied without knowledge as to the management (fertilisation) practices that are common in a region, but lacks the possibility to assess effects of future changes in e.g. N fertiliser availability and changes in management. The crop model GEPIC (Liu et al., 2007) -like LPJmL  -has LAI max and HI as input parameters but also includes N dynamics with empirical response functions as modifiers of productivity and yields. In contrast to LPJmL, GEPIC is based on a site-scale model approach which can be extended in space via a GIS interface that facilitates input of spatially explicit management, where available. WOFOST, a detailed crop-growth model that has recently been expanded to work at larger regions (Boogaard et al., 2013), has demonstrated skills in simulating local and regional yields with a mechanistic approach to modelling crop photosynthesis as well as physiology, but the model approach assumes a continuation of current management practices, and thus cannot be applied for future simulations.
LPJ-GUESS shares several aspects of the approaches of these models: the sowing algorithm (LPJmL Waha et al., 2011), temperature limits for different CFTs (GEPIC), and a large portion of the mechanistic formulations of crop physiology (WOFOST and GECROS are "School of de Wit" models, van Ittersum et al., 2003). Because of the dynamic PHU calculations (Lindeskog et al., 2013) and dynamic sowing and harvest calculations (Waha et al., 2011), LPJ-GUESS can be applied to all areas where suitable temperature and soil moisture conditions allow wheat growth. In addition, we have shown here also that it is possible to find a suitable timing of N fertilisation on the regional scale by relating the N application to those stages in the crop growing period when it is most needed. For globally applicable models this is an important result, since information on fertiliser application often includes total amounts per year, but typically lacks information about the seasonal distribution.

Conclusions
The approach chosen here to implement C-N dynamics in the crop module of LPJ-GUESS seeks to adopt mechanistic process implementations, which has been advocated in literature to be able to fully capture the effects of climate change on ecosystems. The modelling framework demonstrably responds realistically to different N fertiliser treatments and [CO 2 ] and produces results that are in line with observations from site scale to a larger region.
These findings support the aim of this study, to find a level of complexity in the implementation of the N management that can be applied on larger regions, and that is ultimately also applicable under climate and other environmental changes. As long as spatially estimated total N applications are available, adopting mean treatment of timing of the N in the model that is based on the development stage appears sufficient for representing the mean and variance of regional yields.
Representing dynamic C-N interactions within a consistent terrestrial modelling framework provides the capacity to predict changes in global C and N pools and fluxes in historic or future land-use change scenarios, as well as to quantify and explore the effect of different managements on the global C and N budgets, considering hindcasts and projections of land-use change (e.g. Hurtt et al., 2011;Klein Goldewijk et al., 2011), climate change and historic or future N fertiliser application rates (e.g. Potter et al., 2010;Bouwman et al., 2013).

S. Olin et al.: Yield and C : N responses to N management and CO 2
Appendix A  Wang and Engel (1998) N b Minimum leaf N content 0.0011 kg N m −2 β 0 and β 1 Regression coefficients k L ∼ k N 0.01, 0.52 Yin et al. (2003) 1982-1983 and 1983-1984. Blue symbols are with a low input of N fertiliser; red are with a medium input of N; and black are with a high input. Open symbols are for the season 1982-1983, and closed symbols are for the season 1983-1984.