The Joint UK Land Environment Simulator (JULES), Model Description-Part 2: Carbon Fluxes and Vegetation Dynamics

The Joint UK Land Environment Simulator (JULES) is a process-based model that simulates the fluxes of carbon, water, energy and momentum between the land surface and the atmosphere. Many studies have demonstrated the important role of the land surface in the functioning of the Earth System. Different versions of JULES have been employed to quantify the effects on the land carbon sink of climate change, increasing atmospheric carbon dioxide concentrations, changing atmospheric aerosols and tropospheric ozone, and the response of methane emissions from wetlands to climate change. This paper describes the consolidation of these advances in the modelling of carbon fluxes and stores, in both the vegetation and soil, in version 2.2 of JULES. Features include a multi-layer canopy scheme for light interception, including a sunfleck penetration scheme, a coupled scheme of leaf photosynthesis and stomatal conductance, representation of the effects of ozone on leaf physiology, and a description of methane emissions from wetlands. JULES represents the carbon allocation, growth and population dynamics of five plant functional types. The turnover of carbon from living plant tissues is fed into a 4-pool soil carbon model. The process-based descriptions of key ecological processes and trace gas fluxes in JULES mean that this community model is well-suited for use in carbon cycle, climate change and impacts studies, either in standalone mode or as the land component of a coupled Earth system model. Correspondence to: D. B. Clark (dbcl@ceh.ac.uk)


Introduction
Terrestrial ecosystems play an important role in land surface energy and trace gas exchange with the atmosphere. They currently absorb almost one third of the anthropogenic carbon dioxide (CO 2 ) emissions Le Quéré et al., 2009), although the locations and mechanisms for these terrestrial carbon sinks are debated and uncertain (Ciais et al., 1995;McGuire et al., 2001;Stephens et al., 2007;Phillips et al., 2009). Furthermore, land-atmosphere exchange of non-CO 2 greenhouse gases, such as Methane (CH 4 ), Ozone (O 3 ) and Nitrous Oxide (N 2 O), affect atmospheric chemistry and climate (Arneth et al., 2010). Vegetation and soils also exert a strong control on the surface energy balance and the physical state of the atmosphere. Anthropogenic climate change has been projected to radically alter the structure and function of terrestrial ecosystems (Cramer et al., 2001;Sitch et al., 2008). Future shifts in vegetation, such as a northward migration of the boreal forest into tundra, are likely to impact the climate via both biogeophysical and biogeochemical feedbacks. This spurred the development of Dynamic Global Vegetation Models (DGVMs; Cox, 2001;Sitch et al., 2003;Prentice et al., 2007) which describe the structure and function of the major global terrestrial ecosystems.
Advances in recent years have seen the inclusion in land surface models of first a carbon cycle  and a nitrogen cycle (Thornton et al., 2007;Sokolov et al., 2008). Using the TRIFFID DGVM coupled to a General Circulation Model (HadCM3LC), Cox et al. (2000) were the first to show the possibility of a positive climate-land carbon cycle Published by Copernicus Publications on behalf of the European Geosciences Union. feedback, with the counteracting effects of climate and atmospheric CO 2 on ecosystem function. A reduction in terrestrial carbon in response to climate change leads to higher atmospheric CO 2 levels, and thus accelerated climate change. This has major policy implications for climate change mitigation . In the C 4 MIP study Friedlingstein et al. (2006) extended this work using 11 coupled climatecarbon cycle models. All models simulated a positive land carbon cycle feedback but of widely varying strengths and there was little consensus among models on the underlying mechanisms.
The land surface scheme used by Cox et al. (2000) was the Met Office Surface Exchange Scheme (MOSES; Cox et al., 1999;Essery et al., 2003) combined with the TRIFFID DGVM (Cox, 2001). The representation of plant and soil processes in this model, and the implications for the modelled carbon cycle, have been the subject of several subsequent studies. In Cox et al. (2000) the positive feedback is associated with enhanced soil respiration, and drought-induced forest dieback in Amazonia Cox et al., 2004Cox et al., , 2008. Subsequent studies investigated the structural uncertainty in future projections associated with the soil carbon representation (Jones et al., 2005), the role of tropical ecosystems in the control of atmospheric CO 2 on the interannual timescales , and evaluated the coupled model against atmospheric data, proposing a prototype benchmarking methodology for coupled climate-carbon cycle models (Cadule et al., 2010). Jones et al. (2005) replaced the one-pool soil decomposition model with the more elaborate 4-pool model of RothC (Jenkinson, 1990;Coleman and Jenkinson, 1999) and concluded that the projection of a positive feedback between climate and carbon cycle is robust, however, the magnitude of the feedback is dependent on the structure of the soil carbon model. The multi-pool carbon dynamics of RothC cause it to exhibit a slower magnitude of transient response to both increased organic carbon inputs and changes in climate compared with the one-pool model. Gedney et al. (2004) developed an interactive wetlands scheme model that was calibrated using present-day atmospheric CH 4 variability. They predicted increases in global CH 4 flux between present day and 2100 of 75 % with an increase in emissions from northern wetlands (> 30 • N) of 100 %, despite an estimated 10 % reduction in wetland extent. This wetland response corresponds to an amplification of the total anthropogenic radiative forcing at 2100 by approximately 3.5-5 %. Sitch et al. (2007) showed how elevated future tropospheric O 3 concentrations would have detrimental effects on plant productivity and reduce the efficiency of the terrestrial biosphere to sequester carbon, constituting a large indirect radiative forcing of tropospheric O 3 on climate. Mercado et al. (2009) showed how changes in surface irradiance over the global dimming and subsequent brightening period, 1960-2000, associated with changes in anthropogenic scattering aerosols and cloud cover, led to enhanced global plant productivity and carbon storage. Scattering aerosols change both the quantity and quality (partitioning between direct and diffuse) of surface irradiance. Diffuse light is able to penetrate further into the canopy than direct light, stimulating production in light-limited understorey leaves. Mercado et al. (2009) found this diffuse radiation fertilisation effect was larger than the negative effect of reduced irradiance on global plant production. However Mercado et al. (2009) also showed local site optima in the relationship between photosynthesis and diffuse light conditions; under heavily polluted or dark cloudy skies, plant productivity will decline as the diffuse radiation effect is insufficient to offset decreased surface irradiance. Huntingford et al. (2011) took advantage of several of these model developments to study how higher atmospheric CO 2 , climate change, higher near-surface O 3 and lower sulphate aerosols affect net primary productivity and runoff. For equivalent amounts of radiative forcing, the different forcing mechanisms varied markedly in their impacts, suggesting a need for more informative metrics of the impact of changing atmospheric constituents that go beyond simple radiative forcing, and underlining the importance of considering a range of trace gases when modelling land processes.
A comprehensive understanding and description of key ecological processes and nutrient cycles is needed in Earth system models. These include the cycles of carbon, nitrogen and phosphorus; the ecophysiological response of vegetation to changes in atmospheric composition (e.g. plant response to elevated CO 2 and O 3 , N deposition, aerosol radiation effects); the response of vegetation and soils to drought and elevated temperatures; wetland processes and methane exchange; permafrost; and wildfire disturbance. Currently, no single land surface model adequately describes all these processes.
This paper describes modelling of carbon fluxes and stores, in the vegetation and soil, as represented in version 2.2 of the Joint UK Land Environment Simulator (JULES). JULES was initially based on MOSES and TRIFFID and consolidates the improved representations of key processes gained from the studies summarised in the preceding paragraphs. A companion paper , hereafter referred to as Part 1) describes how JULES models fluxes of heat and moisture. Although they are presented separately, the fluxes of moisture and carbon are intimately linked, in particular through the stomatal resistance of the vegetation. The performance of JULES is assessed in Blyth et al. (2010).
Section 2 provides a brief overview of JULES before Sect. 3 describes the photosynthesis model, which has been substantially augmented since Cox et al. (1999) with the addition of an explicit description of the interception of direct and diffuse light at different canopy levels, leading to a multi-layer approach to scaling photosynthesis from leaf to canopy scale. The parameterisations of plant respiration and the effect of ozone on leaf photosynthesis are also covered in that section. The phenology model described in Sect. 4 is essentially unchanged since Cox et al. (1999). Section 5 gives details of the dynamic vegetation model, TRIF-FID (Cox, 2001). Section 6 outlines the simulation of soil carbon, which has changed with the introduction of a 4-pool model and the possibility of choosing between alternative descriptions of the response of heterotrophic respiration to soil temperature. A parameterisation of methane emissions from wetlands is also presented. Finally Sect. 7 describes the approach typically used to bring the initial vegetation and soil carbon pools to an equilibrium state.

Model overview
JULES describes the vegetation in a gridbox using a small number of Plant Functional Types (PFTs). The default is to use five PFTs: broadleaf trees, needleleaf trees, C 3 (temperate) grasses, C 4 (tropical) grasses and shrubs. The surface fluxes of CO 2 associated with photosynthesis and plant respiration are calculated in the physiology component of JULES, as described in Sect. 3 on each JULES timestep (typically 30 to 60 min). The accumulated carbon fluxes are passed to the vegetation dynamics model (TRIFFID, described in Sect. 5) and the area covered by each PFT is updated on a longer timestep (typically 10 days) based on the net carbon available to it and on the competition with other vegetation types, which is modelled using a Lotka-Volterra approach (Cox, 2001). Leaf phenology (bud-burst and leaf drop) is updated on an intermediate timescale of 1 day, using accumulated temperature-dependent leaf turnover rates (Sect. 4). Litterfall from vegetation is input to a model of soil carbon (Sect. 6) which calculates the rate of microbial soil respiration and the consequent flux of CO 2 back to the atmosphere. This part of the model has changed since Cox et al. (1999) with the introduction of a 4-pool model and the possibility of choosing between alternative descriptions of the response of heterotrophic respiration to soil temperature. Methane emissions from wetlands are also calculated. After each call to TRIFFID the land surface parameters required by JULES (e.g. albedo, roughness length) are updated based on the new vegetation state, so that changes in the biophysical properties of the land surface, as well as changes in terrestrial carbon, may feed back onto the atmosphere. The land surface parameters are calculated as a function of the type, height and leaf area index of the vegetation, as described in Sect. 5.2.
The state (or prognostic) variables required to describe the vegetation and soil carbon in JULES are presented in Table 1. Further surface state variables which affect the terrestrial carbon cycle, such as soil moisture and soil temperature, are discussed in Part 1. Appendix A lists the variables used in this paper, along with their units.

Photosynthesis
The photosynthesis model used in JULES is based upon the observed processes at the leaf scale, which are then scaled up to represent the canopy. There are several options available in JULES for the treatment of radiation interception and scaling up to the canopy scale, from a simple big leaf approach to a multi-layer canopy.

Leaf biochemistry
JULES uses the biochemistry of C 3 and C 4 photosynthesis from Collatz et al. (1991) and Collatz et al. (1992), as described by Sellers et al. (1996) and Cox et al. (1999), to determine potential (unstressed by water availability and ozone effects) leaf-level photosynthesis. This is calculated in terms of three potentially-limiting rates: 1. Rubisco-limited rate (W c ) where V cmax (mol CO 2 m −2 s −1 ) is the maximum rate of carboxylation of Rubisco, c i (Pa) is the leaf internal carbon dioxide partial pressure, (Pa) is the CO 2 compensation point in the absence of mitochondrial respiration, O a (Pa) is the partial pressure of atmospheric oxygen, and K c and K o (Pa) are the Michaelis-Menten parameters for CO 2 and O 2 , respectively.
2. Light-limited rate (W l ) where α is the quantum efficiency of photosynthesis (mol CO 2 mol −1 PAR), ω is the leaf scattering coefficient for PAR and I par is the incident photosynthetically active radiation (PAR, mol m −2 s −1 ). where P * is the surface air pressure.
The default values of PFT-specific parameters for leaf biochemistry and photosynthesis are given in Table 2. The parameters V cmax , K o , K c , and are all temperature dependent. JULES uses the temperature dependencies from Collatz et al. (1991Collatz et al. ( , 1992. V cmax at any desired temperature is calculated from the maximum rate of carboxylation of the enzyme Rubisco at 25 • C (V cmax25 ) assuming an optimal temperature range as defined by PFT-specific values of parameters, T upp and T low , as: where T c is canopy (leaf) temperature ( • C) and f T is the standard Q 10 temperature dependence: The default value of Q 10 leaf is 2. V cmax25 is assumed to be linearly related to leaf nitrogen concentration, n l : where n e is a constant that has values of 0.0008 and 0.0004 mol CO 2 m −2 s −1 kg C (kg N) −1 for C 3 and C 4 plants, respectively. These values were derived from Schulze et al. (1994) assuming that leaf dry matter is 40 percent carbon by mass and that the maximum rate of photosynthetic uptake is 0.5 V cmax for C 3 plants and equals V cmax for C 4 plants (Cox, 2001). n l is set equal to n 0 , the leaf N concentration at the top of the canopy, unless variation within the canopy is specified (see Sect. 3). , the photorespiration compensation point, is found as: where τ is the Rubisco specificity for CO 2 relative to O 2 : with Q 10 rs = 0.57.
K c and K o are calculated as: with Q 10 Kc = 2.1 and Q 10 Ko = 1.2. The rate of gross photosynthesis (W ) is calculated as the smoothed minimum of the three potentially-limiting rates: where W p is the smoothed minimum of W c and W l , and β 1 = 0.83 and β 2 = 0.93 are "co-limitation" coefficients. The smaller root of each quadratic is selected. Leaf dark respiration (R d ) is calculated as: where f dr is the dark respiration coefficient. The net potential (i.e. unstressed) leaf photosynthetic carbon uptake (A p ) is then calculated as: Leaf photosynthesis is linked to stomatal conductance via the internal CO 2 concentration, which is calculated using the Jacobs (1994) formulation. The Jacobs formulation shares similarities with the stomatal conductance formulations of Ball et al. (1987) and Leuning (1995). A description of the coupled stomatal conductance-photosynthesis model is given in Part 1, with further details in Cox et al. (1998Cox et al. ( , 1999. To account for soil moisture stress, the potential (nonstressed) leaf photosynthesis A p is multiplied by a soil water factor (Cox et al., 1998): where A l is leaf-level photosynthesis. (Note that the effect of O 3 is also included as a factor on the right-hand side of this equation -see Sect. 3.3 -but is omitted here for clarity.) β is the moisture stress factor which is related to the mean soil moisture concentration in the root zone, θ, and the critical and wilting point concentrations, θ c and θ w , defined as the moisture levels at which photosynthesis first falls below the potential rate and is zero respectively, as follows:

Scaling photosynthesis from leaf to canopy
The description of within-canopy radiation interception and scaling from leaf-to canopy-level photosynthesis has been developed considerably since Cox et al. (1999). JULES2.2 includes a process-based scaling-up of leaf-level photosynthesis to the canopy level, with alternative methods to calculate canopy radiation interception and canopy-level photosynthesis. There are two options available in JULES for scaling up from the leaf-level to the canopy scale: (i) the canopy is considered as a big leaf and (ii) a multi-layer canopy.
Within the multi-layer option, JULES has four variations that depend on considerations of vertical gradients of canopy photosynthetic capacity, inclusion of light inhibition of leaf respiration, inclusion of sunfleck penetration and splitting canopy layers into sunlit and shaded leaves. All options are described below and summarised in Table 3.

Radiation interception and scaling up to canopy level
In the big leaf approach, incident radiation attenuates through the canopy following Beer's law (Monsi and Saeki, 1953): where I c is irradiance beneath the canopy, I o irradiance at the top of the canopy, k is a light extinction coefficient and L c is the canopy leaf area index. Leaf-level photosynthetic capacity is assumed to vary proportionally with the vertical distribution of irradiance (Sellers et al., 1992), therefore leaf photosynthesis can also be expressed as a function of the top of the canopy leaf photosynthesis (A o ), leaf area index (L) and the light extinction coefficient: Canopy photosynthesis is calculated as the integral of leaflevel photosynthesis over the entire canopy leaf area index: Canopy-level conductance and respiration are estimated using similar expressions.

Radiation interception
In all variants of the multi-layer approach (options 2-5 in Table 3), the canopy is divided into a number of layers (n, typically 10) of equal leaf area increments dL c = L c /n. JULES adopts the two-stream approximation of radiation interception from Sellers (1985) to calculate surface spectral albedos (Essery et al., 2001) and the absorbed incoming radiation for each canopy layer. The absorbed incident PAR at each layer varies with solar zenith angle, incident direct and diffuse radiation at the top of the canopy, canopy leaf angle distribution and leaf radiation properties in the visible and near-infrared wavebands. JULES explicitly describes absorption and scattering of both direct and diffuse radiation fluxes separately in the visible and near-infrared wavebands at each canopy layer, which leads to the calculation of upward and downward diffuse fluxes of scattered direct beam radiation (I dir ↑ i , I dir ↓ i ) and scattered diffuse radiation (I dif ↑ i , I dif ↓ i ) per canopy layer, normalised by the incident direct and diffuse fluxes respectively above the canopy. The normalised fluxes are used to calculate the direct and diffuse fractions of absorbed incident PAR, (FAPAR DIR i and FAPAR DIF i ), at each canopy layer i: Multi-layer radiation with Two stream Constant through canopy no two classes (sunlit and shaded) for photosynthesis 4 Multi-layer Two stream Decreases through canopy yes 5 Multi-layer including sunlit and Two stream with Decreases through canopy yes shaded leaves in each layer sunfleck penetration A comparison of the vertical profile of absorbed incident PAR calculated with the two-stream approach against the profile estimated with Beer's law showed that the results were similar only when the incident PAR was a direct beam coming from a high sun angle, otherwise the fraction of absorbed PAR at any canopy level is higher when calculated using Beer's law (Jogireddy et al., 2006). The two-stream approach provides a vertical profile of intercepted radiation within the canopy which allows estimation of photosynthesis and leaf respiration for each leaf area increment within the canopy. When option 3 is selected, rather than calculating photosynthesis for each canopy layer, the leaves in each layer are considered to be either lightlimited or not light-limited (sunlit or shaded) according to whether the light exceeds a threshold (Jogireddy et al., 2006). The photosynthesis calculations are performed separately for each class (sunlit or shaded) using the average light in the class.

Sunfleck penetration
A further improvement to the estimation of absorbed radiation fluxes within the canopy considers penetration of sunflecks through the canopy (option 5 in Table 3), which corresponds to the direct component of the direct beam radiation, i.e. it excludes the scattering component. Such a term is not included in Eq. (20). Attenuation of I b , the non-scattered incident beam radiation per unit leaf area at canopy depth L, normalised by the incident direct beam radiation above the canopy, is calculated as (Dai et al., 2004): where (1 − ω) is the non-scattered part of the incident beam (i.e. what is absorbed) and k b is the canopy beam radiation extinction coefficient. Following Dai et al. (2004) as implemented in Mercado et al. (2009), radiation fluxes are split into direct beam radiation, scattered direct beam and diffuse radiation, and it is assumed that sunlit leaves absorb all types of radiation while shaded leaves absorb only diffuse radiation. The fraction of sunlit leaves (f sun ), is defined as: For each canopy layer i with leaf area increment dL c , the fraction of sunlit leaves, fraction of absorbed direct beam radiation (I b i ), fraction of scattered direct beam (I bs i ) and fraction of absorbed diffuse radiation (I d i ) are: The fractions of the incident radiation above the canopy which are absorbed by sunlit leaves (I sun i ) and shaded leaves (I sh i ) in each leaf area increment within the canopy are calculated as: where f d is the fraction of PAR which is diffuse radiation. I sun i and I sh i are used to calculate the radiation absorbed in each canopy layer by sunlit and shaded leaves by multiplying by the incident radiation above the canopy, and thus to estimate photosynthesis from sunlit (A sun i ) and shaded leaves (A sh i ) for each canopy layer.

Scaling up to canopy-level
For all multi-layer options, canopy-scale fluxes are estimated as the sum of the leaf-level fluxes in each layer, scaled by leaf area. Hence canopy-level photosynthesis is estimated from layer leaf-level photosynthesis (A li ) as follows: with A i as photosynthesis from each canopy layer. When including sunflecks and accounting explicitly for photosynthesis by sunlit and shaded leaves (option 5), A i is calculated as: Canopy-level respiration and conductance are estimated from layer values in a similar manner.

Photosynthetic capacity at each canopy layer
The multi-layer scheme has been applied and tested against eddy correlation flux measurements (Mercado et al., 2007) using different assumptions for the vertical distribution of leaf nitrogen. Such a distribution is a proxy for the vertical distributions of photosynthetic capacity V cmax and leaf respiration through the canopy. The distributions tested were a uniform distribution with leaf N constant through the canopy, and a distribution with leaf N decreasing from top to bottom of the canopy. In the latter case, the vertical profiles of leaf N, photosynthetic capacity and leaf respiration within the canopy were estimated following de Pury and Farquhar (1997) using measured vertical profiles from a rainforest site in the Amazon Basin (Carswell et al., 2000) and prescribed in JULES to decrease exponentially with increasing canopy depth (Mercado et al., 2007) for options 4 and 5 in Table 3. Photosynthetic capacity at each canopy layer i is calculated assuming that the reference value varies as: with n 0 the leaf N concentration at the top of the canopy and k n a nitrogen profile coefficient estimated to be 0.78. (Note that this is the same form as Eq. (6) but with the vertical variation of leaf nitrogen specified.) Vertical profiles of V cmax remain to be tested further and evaluated for other vegetation types.
The multi-layer options 4 and 5 in JULES also account for inhibition of leaf respiration in light. Mercado et al. (2007) tested the inclusion of inhibition of leaf respiration by light from Brooks and Farquhar (1985) as implemented by Lloyd et al. (1995) for a rainforest site in the Amazon. Once JULES was correctly parameterised for canopy photosynthetic capacity at this site, the inclusion of this inhibition allowed much better predictions of observed rates of net photosynthetic uptake. Light inhibition follows Mercado et al. (2007) when option 4 is used, while option 5 uses a 30 % inhibition of leaf respiration at irradiance levels greater than 10 µmol quanta m −2 s −1 (Atkin et al., 1997(Atkin et al., , 2000(Atkin et al., , 2006.

Assessment of big-leaf and multi-layer approaches
JULES was evaluated using the multi-layer approach and eddy correlation data for a temperate coniferous forest site in the Netherlands (Jogireddy et al., 2006) and a tropical broadleaf rainforest site in the Brazilian Amazon (Mercado et al., 2007). Both studies demonstrated the superior performance of the multi-layer approach with the two-stream canopy radiation interception compared to the big-leaf approach in simulating canopy scale photosynthetic fluxes, specifically both the simulated light response and diurnal cycles of photosynthesis. A further advantage of a layered scheme is that it differentiates between direct and diffuse radiation, which is not possible using the Beer's law approach. Further evaluation of the multi-layer approach at eddy correlation sites and globally is presented in Blyth et al. (2010).
Using a 10 layer model, Jogireddy et al. (2006) examined the impact of calculating photosynthesis at a coniferous forest site for only two classes (light-limited or not) of leaves rather than for each of 10 canopy layers (options 3 and 2 in Table 3). The 2-class option gave a good fit to the 10layer simulation and, because this option is computationally efficient, it was recommended for large scale applications if computer resources are limited.
Allowing leaf nitrogen, canopy photosynthetic capacity and leaf respiration to vary through the canopy (options 4 and 5) provides a more realistic representation of canopy and total plant respiration in JULES; the description of stem and root respiration in JULES is a dependent function of canopy respiration and their respective nitrogen contents. This is especially apparent in tropical ecosystems, where simulations that assume a uniform vertical distribution of leaf N, and therefore photosynthetic capacity, produce very large respiration fluxes from leaves in the shaded understorey. Observations of decreased leaf N and photosynthetic capacity within canopies (Meir et al., 2002) and decreased leaf respiration in the light (Brooks and Farquhar, 1985;Atkin et al., 1998;Hoefnagel et al., 1998;Atkin et al., 2000) support their inclusion into JULES. Figure 1 shows evaluation of Gross Primary Productivity (GPP) simulated by JULES against eddy correlation data from a temperate broadleaf (Knohl et al., 2003) and a needle leaf (Rebmann et al., 2010) site. Observations are compared with runs of JULES that used the big leaf approach and the multi-layer option with vertical variation of photosynthetic capacity and inclusion of light inhibition of leaf respiration (options 1 and 4 in Table 3). Simulated photosynthesis using the big leaf approach shows light saturation at low levels of radiation (left panels) and a "flat" response around midday (right panels), unlike the observations. This is because with the big leaf approach the simulated photosynthesis is light-limited only with very low levels of radiation and thereafter is limited by carboxylation of Rubisco. With the multi-layer approach for light interception and photosynthesis, JULES simulates competition between light-limited and Rubisco-limited photosynthesis at each canopy layer, resulting in increased Rubisco limitation towards the top of the canopy and increased light limitation lower in the canopy. In addition to a better representation of light response and diurnal cycles of canopy photosynthesis (Fig. 1), the multi-layer approach also leads to improved simulation of stomatal and canopy conductance (not shown).
The inclusion of sunflecks by Mercado et al. (2009) provided a platform to study the differential effects of direct and diffuse radiation on carbon and water exchange. The 10layer model, including sunfleck penetration, the vertical decrease in photosynthetic capacity within the canopy and the inhibition of respiration in the light, is the recommended default setting (option 5 in Table 3) for applications at all scales from individual sites to global modelling, because it can provide the most realistic representation of plant physiological processes. However, the inclusion of a vertical profile of photosynthetic capacity through the canopy is likely to require specific parameterisations for each PFT. Field data describing how nutrients and related physiological processes vary through the canopy are available for relatively few vegetation types and locations. The existing description in JULES is largely based on analysis of data for broadleaf trees in Amazonia (Carswell et al., 2000, as implemented in JULES by Mercado et al., 2007). Based on data from a variety of Amazon broadleaf trees, a recent study (Lloyd et al., 2010) derived an equation describing the relationship between vertical gradients in photosynthetic capacity within tree canopies and the photosynthetic capacity of the upper leaves. We anticipate using this type of information in future development of JULES, although field studies of other vegetation types and in other regions are also required. In addition, the description of inhibition of leaf respiration in light is based on limited data for a small number of species. More observational data are needed to refine this inhibition and how it varies across plant functional types.

Ozone effects on photosynthesis
Ozone causes cellular damage inside leaves which adversely affects plant production, reduces photosynthetic rates and requires increased resource allocation to detoxify and repair leaves (Ashmore, 2005). JULES uses a flux-gradient approach to model ozone damage, following Sitch et al. (2007). It is assumed that ozone suppresses the potential net leaf photosynthesis in proportion to the ozone flux through stomata above a specified critical ozone deposition flux, so that the actual leaf-level net photosynthesis (A l ) is given by: where A * l is leaf-level net photosynthesis in the absence of O 3 effects and the reduction factor: is given by the instantaneous leaf uptake of O 3 (F O 3 , nmol m −2 s −1 ) above a PFT-specific threshold F O 3 crit , multiplied by a PFT-specific parameter (a), following Pleijel et al. (2004). The cumulative effect of leaf damage and early senescence is implicitly accounted for in our calibration of a (by giving reduced photosynthesis during the growing season). The flux F O 3 is calculated by analogy with Ohm's law as: where [O 3 ] is the molar concentration of O 3 at reference level (nmol m −3 ), r a is the aerodynamic and boundary layer resistance between leaf surface and reference level (s m −1 ), g l is the leaf conductance for H 2 O (m s −1 ), and κ O 3 = 1.67 is the ratio of leaf resistance for O 3 to leaf resistance for water vapour. The uptake flux is dependent on the stomatal conductance, which is dependent on the photosynthetic rate in JULES. Given g l is a linear function of photosynthetic rate (Eq. 13, Cox et al., 1999), from Eq. (34) it follows that: where g * l is the leaf conductance in the absence of O 3 effects. Through this mechanism the direct effect of O 3 deposition on photosynthesis also leads to a reduction in stomatal conductance. As the O 3 flux itself depends on the stomatal conductance, which in turn depends upon the net rate of photosynthesis (Cox et al., 1999), the model requires a consistent solution for the net photosynthesis, stomatal conductance and the ozone deposition flux. Equations (35)-(37) produce a quadratic in F which is solved analytically.
Data from field observation Pleijel et al., 2004) are used to calibrate plant-ozone effects for the five standard PFTs described by JULES (see Sitch et al., 2007 for details of the calibration procedure and Table 4 for parameter values). Sitch et al. (2007) presented "high" and "low" parameterisations for each PFT to represent species sensitive and less sensitive, respectively, to ozone effects.
The default parameter values in JULES are the "low" sensitivity values. Threshold values, F O 3 crit , are taken at 1.6 and 5 nmol m −2 s −1 for the woody and grass PFTs, respectively. Although a threshold of 5 implies a smaller O 3 dose for grasses, the gradient of the dose-response function (a), is larger, and therefore grasses may become more sensitive to ozone exposure than trees at high ozone concentrations. For shrubs we assume the same plant-ozone sensitivity as broadleaf trees.

Plant respiration
The gross primary productivity ( G ) is: where βR dc is the soil moisture-modified canopy dark respiration. The net primary productivity ( ) is: where R p is plant respiration. R p is split into maintenance and growth respiration (Cox et al., 1999): Growth respiration is assumed to be a fixed fraction of the net primary productivity: The growth respiration coefficient (r g ) is set to 0.25 for all PFTs. Leaf maintenance respiration is equivalent to the soil moisture-modified canopy dark respiration, βR dc , while root and stem respiration are assumed to be independent of soil moisture, but to have the same dependences on nitrogen content and temperature. Thus total maintenance respiration is given by: where N l , N s and N r are the nitrogen contents of leaf, stem and root, and the factor of 0.012 converts from (mol CO 2 m −2 s −1 ) to (kg C m −2 s −1 ). The nitrogen contents are given by: where n m is the mean leaf nitrogen concentration (kg N (kg C) −1 ), R and S are the carbon contents of root and respiring stem, L is the canopy leaf area index and σ l (kg C m −2 per unit of LAI) is the specific leaf density. The nitrogen concentrations of roots and stem are assumed to be fixed (functional type dependent) multiples, µ rl and µ sl , of the mean leaf nitrogen concentration: µ rl = 1.0 for all PFTs, µ sl = 0.1 for woody plants (trees and shrubs) and µ sl = 1.0 for grasses. The respiring stemwood is calculated using a "pipe model" approach (Shinozaki et al., 1964a,b) in which live stemwood is proportional to leaf area L and canopy height, h: The constant of proportionality η sl is approximated from Friend et al. (1993).

Leaf phenology
Leaf phenology is modelled as described in Cox (2001). Leaf mortality rates, γ lm , are assumed to be a function of temperature, increasing from a minimum value of γ 0 , as the leaf temperature drops below a threshold value, T off : d T is the rate of change of γ lm with respect to temperature, and equals 9 by default (Table 5) meaning that the leaf turnover rate increases by a factor of 10 when the temperature drops 1 • C below T off . Equation (47) describes how leaf mortality varies with temperature, but it is not sufficient to produce realistic phenology. A variable, p, is introduced which describes the phenological status of the vegetation: where L is the actual LAI of the canopy, and L b is the balanced (or seasonal maximum) LAI as updated by the dynamic vegetation model (TRIFFID) via the inverse of Eq. (58). The phenological status, p, is updated typically on a daily basis assuming: -leaves are dropped at a constant absolute rate (γ p L b ) when the mean value of leaf turnover, as given by Eq. (47), exceeds twice its minimum value -budburst occurs at the same rate when γ lm drops back below this threshold, and "full leaf" is approached asymptotically thereafter: where γ p = 20 yr −1 . The effective leaf turnover rate, γ l , as used later to calculate litterfall (see Eq. 59), must also be updated to ensure conservation of carbon when phenological changes are occurring: Taken together, Eqs. (47), (49) and (50) amount to a "chilling-days" parameterisation of leaf phenology. Colddeciduous behaviour can effectively be disabled for any of the PFTs by setting parameter d T to zero for that PFT. A similar approach exists in the JULES code for drought-deciduous phenology involving equations similar to Eq. (47) but for leaf turnover as a function of soil moisture, with parameters M off and d M . However, this is considered to be still under development and default parameters of d M = 0 for all PFTs mean that this is effectively switched off for general use. Calculation of leaf phenology is independent of the calculation of the evolution of vegetation coverage and can be included even when the dynamic vegetation component, TRIF-FID, is turned off.

Vegetation dynamics
The dynamic vegetation model used in JULES is TRIFFID (Cox, 2001). Each land grid box is assumed to be either completely covered by permanent ice (in which case TRIFFID is not used) or to have no ice cover. At these non-ice points, the urban and lake surface types (if present) are assumed to have time-constant fractions while the 5 PFTs compete for the remaining coverage as simulated by TRIFFID. 1 The final surface type, bare soil, is the remaining space after simulating the coverage of the vegetation types.
The 5 PFTs of TRIFFID were chosen as a minimal set to represent the variation in vegetation structure (e.g. canopy height, root depth) and function (e.g. C 3 versus C 4 photosynthesis) for inclusion of both biophysical and biogeochemical vegetation feedbacks in Earth System Models. The number of PFTs defined in DGVMs typically ranges from 2 (Brovkin et al., 1997) to 10  and depends on the ecological processes explicitly represented in the model and the availability of field data for defining parameter values, which is often incomplete. The latter is especially relevant to land surface models within Earth System Models as these require detailed information on both ecophysiological and physical parameters. The original 5 PFTs of TRIFFID represented a pragmatic choice balancing comprehensiveness against computational expense and limited data availability for the paramaterization of each PFT. The choice and definition of PFTs is an area of active research within the land surface modelling community.

Vegetation growth and competition
The vegetation carbon density, C v , and fractional coverage, ν, of a given PFT are updated based on the carbon balance of that PFT and on competition with other PFTs: 1 Note that although the number of PFTs in JULES is generally flexible, a run that uses TRIFFID to simulate vegetation competition can only use the 5 standard PFTs (broadleaf trees, needleleaf trees, C 3 grasses, C 4 grasses and shrubs) as the competition coefficients are hardwired.
where is the net primary productivity per unit vegetated area of the PFT in question, l is the local litterfall rate, and γ v is the large-scale disturbance rate. ν * = max{ν, }, where = 0.01 is the "seed fraction". Under most circumstances ν * is identical to the a real fraction, ν, but each PFT is "seeded" by ensuring that ν * never drops below the seed fraction. A fraction λ of this NPP is utilised in increasing the fractional coverage (Eq. 52), and the remainder increases the carbon content of the existing vegetated area (Eq. 51). Default values of PFT-specific parameters for TRIFFID are given in Table 6.
The competition coefficients, c ij , represent the impact of vegetation type j on the vegetation type of interest. These coefficients all lie between zero and unity, so that competition for space acts to reduce the growth of ν that would otherwise occur (i.e. it produces density-dependent litter production). Each PFT experiences "intra-species" competition, with c ii = 1 so that the vegetation fraction is always limited to be less than one. Competition between natural PFTs is based on a tree-shrub-grass dominance hierarchy, with dominant types i limiting the expansion of sub-dominant types j (c j i = 1), but not vice-versa (c ij = 0). However, the tree types (broadleaf and needleleaf) and grass types (C 3 and C 4 ) co-compete with competition coefficients dependent on their relative heights, h i and h j : The form of this function ensures that the i-th PFT dominates when it is much taller, and the j -th PFT dominates in the opposite limit. The factor of 20 was chosen to give cocompetition over a reasonable range of height differences. Some allowance is made for agricultural regions, from which the woody types (i.e. trees and shrubs) are excluded, and only C 3 and C 4 grasses can grow. These can be interpreted as "crops" but are not simulated differently in agricultural or non-agricultural regions. The λ partitioning coefficient in Eqs. (51) and (52) is assumed to be piecewise linear in the leaf area index, with all of the NPP being used for growth for small LAI values, and all the NPP being used for "spreading" for large LAI values: where L max and L min are parameters describing the maximum and minimum leaf area index values for the given plant functional type, and L b is the "balanced" LAI which would be reached if the plant was in "full leaf". The actual LAI depends on L b and the phenological status of the vegetation type, which is updated as a function of temperature (see Sect. 4).
As the DGVM component of JULES, TRIFFID simulates the evolution of both the fractional coverage of vegetation and of terrestrial carbon storage. An option exists to disable updating of the vegetation fractions so that JULES can be run with fixed surface cover but evolving carbon storage.
Changes in vegetation carbon density, C v , are related allometrically to changes in the balanced LAI, L b . First, C v is broken down into leaf, L, root, R, and total stem carbon, W: Each components is then related to L b . Root carbon is set equal to leaf carbon, which is itself linear in LAI, and total stem carbon is related to L b by a power law (Enquist et al., 1998): where a wl and b wl are PFT-dependent parameters.
Values of canopy height, h, are found directly from W as described in Sect. 5.2.
The local litterfall rate, l , in Eq. (51), consists of contributions from leaf, root and stem carbon: where γ l , γ r and γ w are turnover rates (yr −1 ) for leaf, root and stem carbon respectively. The leaf turnover rate is calculated to be consistent with the phenological module as described in Sect. 4. There is an additional litter contribution arising from large-scale disturbance which results in loss of vegetated area at the prescribed rate γ ν , as represented by the last term on the right-hand side of Eq. (52). Using a simplified version of TRIFFID, Huntingford et al. (2000) described the response of a dominant vegetation type in different initial climatic regimes to prescribed changes in atmospheric CO 2 concentration and temperature, and noted that the model was capable of a rich range of possible behaviours of relevance to future global change.

Updating of biophysical parameters
The land-surface parameters required by JULES are recalculated directly from the LAI and canopy height of each PFT, each time the vegetation cover is updated. Values of canopy height, h, are derived by assuming a fixed ratio, a ws , of total stem carbon to respiring stem carbon: where we assume a ws = 10.0 for woody plants and a ws = 1.0 for grasses (Friend et al., 1993, and Table 7). Combining with Eqs. (58) and (46) enables canopy height to be diagnosed directly from the total stem biomass: Given the canopy height and LAI, the values of biophysical parameters are calculated as described in Part 1. Table 8. Default values of pool-specific parameters for soil carbon. The pools are decomposable and resistant plant material (DPM, RPM), biomass (BIO) and humus (HUM).

Soil carbon
Total soil carbon storage, C s , is increased by the total litterfall, c , and reduced by microbial soil respiration, R s , which returns CO 2 to the atmosphere: In each gridbox, the total litterfall is made-up of the areaweighted sum of the local litterfall from each PFT (as given by Eq. 59), along with terms due to the large-scale disturbance rate, γ ν , and PFT competition (from the TRIFFID dynamic vegetation model): The competition term (last term on the right-hand side of Eq. 63) is derived by imposing carbon conservation on the soil-vegetation system as described by Eqs. (51), (52) and (62). It implies that the NPP of each PFT will be lost entirely as litter once the PFT occupies all of the space available to it (i.e. when j c ij ν j = 1). If the TRIFFID DGVM is used, four soil pools are modelled, otherwise a single pool with a fixed (in time) value is used (Cox, 2001) to calculate soil respiration. However, in this case the soil carbon pool is not updated and net carbon fluxes may not be in balance. It is recommended that TRIF-FID is used if simulation of soil carbon fluxes is required.

Implementation of the RothC soil carbon model
The soil carbon model comprises 4 carbon pools and follows the formulation of the RothC soil carbon scheme (Jenkinson, 1990;Coleman and Jenkinson, 1999). Plant input is split between two carbon pools of decomposable (DPM) and resistant (RPM) plant material depending on the overlying vegetation type, with grasses providing a higher fraction of decomposable litter input and trees a higher fraction of resistant litter input. The other two carbon pools are microbial biomass (BIO) and long-lived humified (HUM) pools.
Respiration from each soil carbon pool is the product of a specific respiration rate, κ si , (the rate of respiration of unit Fig. 2. Comparison of alternative forms for the sensitivity of specific soil respiration to soil temperature in JULES. The dashed line shows a Q 10 form with Q 10 = 2, the solid line shows the form from the RothC model (Jenkinson, 1990). mass of that pool under standard conditions), the pool size, C i , and several rate modifying factors: Values of κ si for each pool are given in Table 8. The rate modifying factors account for the effects of soil temperature (T soil ), soil moisture (s), and vegetation fractional cover (v) on heterotrophic respiration. The form of the temperature rate modifier is controlled by a switch that allows the user to select between a Q 10 function: where the default value is Q 10 soil = 2.0, and the RothC temperature function (Jenkinson, 1990): F T (T soil ) = 47.9 1 + e 106/(T soil −254.85) The temperature of the top soil layer (typically 10 cm deep) is used in both cases. Figure 2 compares these alternative temperature functions. The effect of soil moisture is described as:  (67) where s and s o are the unfrozen soil moisture content of the top soil layer and the optimum soil moisture, both expressed as a fraction of saturation. s o = 0.5(1+s w ) and s min = 1.7s w , where s w is the soil moisture at wilting point. The general form of the moisture function from Cox (2001) has been retained in preference to the RothC moisture function because of its ability to simulate reduced respiration in very wet soils. However, it has been revised (the original set s min = s w ) so as to simulate a greater sensitivity of respiration reduction in dry soils, which gave a better fit to the observed site-level seasonal cycle of respiration. The importance of moisture controls on future soil carbon is discussed in Falloon (2009) and. The effect of vegetation cover is described as: varying linearly from 0.6 under fully vegetated soil to 1 under completely bare soil. The fraction of litter that is decomposable plant material is: where α dr controls the ratio of litter input to DPM and RPM, taking values of 0.25 for trees, 0.33 for shrubs, 0.67 for natural grass and 1.44 for crops. Carbon from decomposition of all 4 carbon pools is partly released to the atmosphere and partly feeds the BIO and HUM pools. The carbon pools are updated according to: where R s is the total respiration rate, summed over the 4 pools. β R depends on soil texture to account for the protective effect of small particle sizes. It is expected that the inclusion of multi-pool dynamics in the soil carbon model will dampen the transient response of soil carbon storage to both changes in litter input and changes in climate (Jones et al., 2005), although the long-term sensitivity will be unchanged if the same Q 10 function of sensitivity to temperature is used.

Methane emissions from wetlands
The gridbox-average methane emission from the wetland fraction of each gridbox is calculated following Gedney et al. (2004) as: where k CH 4 is a global constant, C eff is the effective substrate availability, Q 10 CH 4 is a temperature-dependent Q 10 factor and T 0 is a reference temperature. f wet is the fraction of the gridbox that is considered to be wetland (i.e. stagnant water) and is calculated using subgrid topographic information, as described in Part 1. The effective Q 10 value is temperature dependent and calculated as: When the four-pool soil carbon model is used the substrate availability is calculated by weighting the size of each pool by its specific respiration rate, otherwise C eff = C s . The default parameter values are k CH 4 = 7.4 × 10 −12 s −1 , T 0 = 273.15 K and Q 10 CH 4 (T 0 ) = 3.7.

Spin-up methodology
Soil carbon and vegetation fractions have timescales of order 100s-1000s yr to reach equilibrium which would necessitate very long spin-up simulations. Hence, TRIFFID was designed to be used in both an "equilibrium" and a "dynamic" mode. The TRIFFID equations to update the plant fractional coverage and leaf area index are written to enable both explicit and implicit timestepping so that their discretisation can be reduced to the Newton-Raphson algorithm for iteratively approaching an equilibrium given external driving conditions. In equilibrium mode TRIFFID is coupled asynchronously to the rest of the model, with accumulated carbon fluxes passed from the physiology component after an integer number of years, typically every 5 yr in order to smooth the effect of interannual variability. If JULES is being run with a single repeating year of meteorology then an equilibrium timestep of 1 yr is sufficient. On each TRIFFID call, the vegetation and soil variables are updated iteratively using an implicit scheme with a long internal timestep. This approach is very effective in producing equilibrium states for the slowest variables (e.g. large soil carbon pools or forest cover). Figure 3 shows the evolution of annual NPP and broadleaf tree fraction when run in equilibrium mode for a single point, with annually repeating meteorology, from arbitrary initial conditions of 0, 35 % and 100 % broadleaf tree cover, and calling TRIFFID once a year. Although the "real" timescale of tree growth would be decades to centuries, the spin-up mode converges to a steady state after just 5 calls to TRIF-FID. In a global simulation this rapid equilibration technique is invaluable.
This implicit equilibrium mode has also been implemented for the 4-pool soil carbon model in JULES (Sect. 6.1), but a drawback is that sub-annual timescales which are important to the small, fast DPM pool are not captured by the multi-year equilibrium approach. The mean DPM implied by a seasonally-varying input of fluxes is not the same as the DPM implied by the annual mean of those fluxes. Hence, the equilibrium mode in JULES produces only an approximation of a steady soil-carbon state. To achieve a full spin-up, either JULES can subsequently be run in dynamic mode until steady state is reached, or the soil carbon model can be run offline from the rest of JULES using carbon fluxes and rate-modifiers as inputs.

Summary
JULES is based on the MOSES land surface scheme (Cox et al., 1999) and the TRIFFID dynamic global vegetation model (Cox, 2001), but with significant improvements. In particular, JULES includes several options for scaling photosynthesis from leaf to canopy scale, with the most complex modelling the profile of light interception through the vegetation with a multi-layer scheme including representation of sunfleck penetration and variation of photosynthetic capacity through the canopy. The coupled model of leaf photosynthesis and stomatal conductance includes representation of the effect of ozone on vegetation physiology. Soil carbon processes are modelled using a 4-pool description and methane emissions from wetlands are also modelled.
The performance of JULES is quantified using a system described by Blyth et al. (2010) which includes measures of land-atmosphere fluxes of CO 2 . Note that Blyth et al. (2010) used an earlier version of JULES than that described here (version 2.1.2 rather than 2.2). This later version includes representation of the effects of ozone on leaf physiology, further options for the calculation of canopy photosynthesis (options 4 and 5 in Table 3), the ability to disable competition between vegetation types in the TRIFFID dynamic vegetation model, and options for a more advanced treatment of urban areas (see Part 1), in addition to bug fixes. The latest version of JULES is version 3.0 in which the land surface model can be coupled to the IMOGEN system (Huntingford et al., 2010), thereby allowing a first-order assessment of how the biogeochemical processes represented in JULES might respond to, and in turn feed back on, a changing climate. The details of version 3.0 are beyond the scope of this paper.
The inclusion of land-atmosphere exchanges of CO 2 , H 2 O, CH 4 and O 3 and aerosol effects in a single model framework is a significant development towards a comprehensive trace-gas enabled land surface model (Arneth et al., 2010). The development of JULES is ongoing, with revised or new representations of several key Earth system processes under consideration. These include a vegetation model based on a statistical approximation of a canopy gap model (Fisher et al., 2010b), a model of fire disturbance based on Thonicke et al. (2010), soil C and N cycles  coupled to a description of plant N uptake (Fisher et al., 2010a), and a model of biogenic isoprene emission (Pacifico et al., 2010).
For further details of JULES, including how to acquire a copy of the code, see http://www.jchmr.org/jules. Non-scattered direct radiation, normalised by the incident direct radiation above the canopy I bs Scattered direct radiation, normalised by the incident direct radiation above the canopy I c W m −2 Irradiance beneath the canopy I d Absorbed diffuse radiation, normalised by the incident direct radiation above the canopy I dif ↑ Upward diffuse flux of scattered diffuse radiation normalised by the incident diffuse flux above the canopy I dif ↓ Downward diffuse flux of scattered diffuse radiation normalised by the incident diffuse flux above the canopy I dir ↑ Upward diffuse flux of scattered direct radiation normalised by the incident direct flux above the canopy I dir ↓ Downward diffuse flux of scattered direct radiation normalised by the incident direct flux above the canopy I o W m −2 Irradiance at the top of the canopy I par mol m −2 s −1 Incident photosynthetically active radiation Symbol Units Definition I sh Fraction of incident radiation above the canopy that is absorbed by shaded leaves I sun Fraction of incident radiation above the canopy that is absorbed by sunlit leaves k Light extinction coefficient k b Direct beam canopy extinction coefficient K c Pa Michaelis-Menten constant for CO 2 k CH 4 s −1 Scaling constant for methane emission k n Nitrogen profile coefficient Ratio of canopy LAI to balanced LAI P * Pa Surface air pressure Q 10 CH 4 Temperature-dependent Q 10 factor for methane emission Q 10 Kc Q 10 value for Michaelis-Menten parameter for CO 2 Q 10 Ko Q 10 value for Michaelis-Menten parameter for O 2 Q 10 leaf Q 10 value for carboxylation of Rubisco Q 10 rs Q 10 value for Rubisco specificity for CO 2 relative to O 2 Q 10 soil Q 10 value for soil respiration R kg C m −2 Carbon content of roots r a s m −1 aerodynamic and boundary layer resistance between leaf surface and reference level R BIO kg C m −2 s −1 Ratio of nitrogen content of roots to that of leaves µ sl Ratio of nitrogen content of roots to that of respiring stem ν Fractional coverage of a vegetation type ν * Fractional coverage of a vegetation type kg C m −2 s −1 Net primary productivity G kg C m −2 s −1 Gross primary productivity σ l kg C m −2 per unit LAI Specific leaf density τ Rubisco specificity for CO 2 relative to O 2 ω Leaf scattering coefficient for PAR