The Joint UK Land Environment Simulator (JULES), model description – part 1: energy and water fluxes

. This manuscript describes the energy and water components of a new community land surface model called the Joint UK Land Environment Simulator (JULES). This is developed from the Met Ofﬁce Surface Exchange Scheme (MOSES). It can be used as a stand alone land surface model driven by observed forcing data, or coupled to an atmospheric global circulation model. The JULES model has been coupled to the Met Ofﬁce Uniﬁed Model (UM) and as such provides a unique opportunity for the research community to contribute their research to improve both world-leading operational weather forecasting and climate change prediction systems. In addition JULES, and its forerunner MOSES, have been the basis for a number of very high-proﬁle papers concerning the land-surface and climate over the last decade. JULES has a modular structure aligned to physical processes, providing the basis for a ﬂexible modelling platform.


Introduction
Traditionally Land Surface Models (LSMs) have been considered as the lower boundary condition for Global Circulation Models (GCMs) and other atmospheric modelling systems.Over the last couple of decades, the importance of the influence that the land surface has on atmospheric Correspondence to: M. Best (martin.best@metoffice.gov.uk)modelling has increased, which has led to additional focus on the complexity and accuracy of LSMs.Models have developed from a simple energy balance with a simple soil scheme (e.g., Deardorff, 1978) through to complex vegetation structures with multiple layer soil hydrology.Examples of currently used land surface schemes include the Interaction Soil-Biosphere-Atmosphere model (ISBA, Noilhan and Planton, 1989), the Canadian Land Surface Scheme (CLASS, Verseghy, 1991;Verseghy et al., 1993), the Tiled ECMWF Scheme for Surface Exchanges over Land model (TESSEL, Viterbo and Beljaars, 1995), the NOAH model (Ek et al., 2003) and the Community Land Model (CLM, Oleson et al., 2010).
The large differences in the response of the surface fluxes to various surfaces has initiated a representation of subgridscale heterogeneity, such as tile or mosaic schemes (e.g., Essery et al., 2003a).Differences at the surface can be caused by their interaction with snow (e.g., snow on top of the surface as with bare soil and short vegetation, or snow under the "surface" as with needleleaf forests), the availability of water at the surface influencing the Bowen ratio (e.g., open water, snow and ice surfaces compared to vegetation and bare soil surfaces), or in the treatment of the carbon cycle for vegetation (e.g., the difference in carbon pathways between C3 and C4 vegetation).Further increases in model resolution, particularly for regional scale operational weather forecasting, open up new challenges in the way we represent the subgridscale heterogeneity at the surface, as the nature of the heterogeneity changes.
Published by Copernicus Publications on behalf of the European Geosciences Union.As the resolution and accuracy of atmospheric modelling systems increases, there is likely to be a need for a wider diversity of land surface processes, such as river flow and flooding, groundwater, or potential crop yields.These new processes present some challenges as model developers will have to acquire new areas of expertise and integrate new science in existing modelling systems.
The development in our understanding of the interactions between the atmosphere and the biosphere for the carbon cycle has begun a new era for science in land surface modelling (e.g., Cox et al., 2000).Current research activities are not limited to the carbon cycle, but are also considering other elements such as the nitrogen cycle, methane and ozone (Gedney et al., 2004;Sitch et al., 2007;Thornton et al., 2007Thornton et al., , 2009;;Sokolov et al., 2008;Fisher et al., 2010;Zaehle et al., 2010).Again, the complexity of these new systems require additional expert knowledge that has traditionally not been held by the original LSM developers.
It is beyond most research and operational centres to have the expertise in such a diverse range of science.Therefore to develop a state of the art LSM requires an alternative perspective to the traditional isolated development of these modelling systems.The development of a community land surface model enables experts in areas of land surface science to contribute towards a leading land surface model, from which all users will benefit.This approach has been adopted with the Community Land Model (CLM) and the NOAH model, and now with the new community land surface model, the Joint UK Land Environment Simulator (JULES).JULES originated from the Met Office Surface Exchange Scheme (MOSES; Cox et al., 1999;Essery et al., 2003a), the land surface model developed at the UK Met Office for applications ranging from operational weather forecasting to Earth system modelling.The forcing data required by JULES (Table 1) are the standard information that would be exchanged when coupled to an atmospheric GCM.Hence, JULES can be linked to the UK Met Office Unified Model (Cullen, 1993) opening up the unique opportunity for the research community to contribute its science into leading operational weather forecasting and climate change prediction systems.
In addition JULES, and its forerunner MOSES, have already been the basis of a number of high-profile papers on the response of land ecosystems to climate (Cox et al., 2000;Gedney et al., 2006;Betts et al., 2007;Sitch et al., 2007;Cox et al., 2008;Mercado et al., 2009).
As well as the initialisation of the prognostic variables within the JULES model (Table 2), ancillary information is required for various soil parameters (Table 3).These data are required for both stand alone and coupled applications.In addition, information on the various parameters used within the JULES model is contained in the user documentation, which is attached as supplementary material to this paper.
JULES has been designed to be a flexible modelling system with a modular structure.This structure is illustrated in Fig. 1, where the connections between the modules show the physical processes that connect these areas.The aim of this modular structure is to make it easy to replace modules or to introduce new modules within the modelling system.For instance, whilst at present JULES can be coupled to an external river flow model via the surface and sub-surface runoff fluxes to simulate river discharge, future versions of JULES will include these processes as new modules, along with other processes such as irrigation and groundwater.
Within the modules there are also various science options (Table 4), which can be selected through a series of switches.In general the options represent subsequent developments and improvements to the physics represented in the model.The use of the scheme within an operational weather forecast model (and its evolution from the MOSES land model which was also used in the same environment) requires that such developments are not just simply replaced, but made available as options to ensure backwards compatibility between model versions.However, this presents an opportunity to analyse how developments have impacted the subsequent performance of such a land model.
In addition to the main science modules within JULES there are also three themes.These themes are not connected by physical processes to the other modules, but do impact on each of them and are critical to ensure that the JULES modelling system remains a flexible, easy to use and develop, openly validated tool that can have identifiable configurations for applied applications.These themes include the technical design of the modelling system, the validation and calibration of all aspects of the model, and setting configurations of the modelling system that are suitable for climate impact studies.The themes surround the science modules in Fig. 1 demonstrating their integrating nature.
This paper, the first of two parts that describe the JULES system, is concerned with the energy and water cycles.The second part describes the additional modules required to represent the carbon cycle (Clark et al., 2011), whilst a companion paper addresses one of the cross cutting themes with benchmarking (Blyth et al., 2011).The sections of this paper describe the modules in Fig. 1 relating to energy and water.Section 2 describes the surface exchange, covering (2.1) Geosci.Model Dev., 4, 677-699, 2011 www.geosci-model-dev.net/4/677/2011/ the energy balance equations, (2.2) the surface resistance of moisture for vegetation, (2.3) evaporation of moisture on the surface in either liquid or solid states, (2.4) how urban areas are represented, and (2.5) the treatment of surface heterogeneity.
Section 3 decribes the processes relating to snow.This includes (3.1) the interaction of snow with vegetation canopies, two methods for modelling the snow on the ground, with either (3.2) zero layer or (3.3) multi-layer models, and (3.4) the representation of snow albedo.
Section 4 deals with soil processes for temperature and moisture.This includes (4.1) the amount of water that reaches the soil surface through vegetation canopies and how this is then distributed into runoff and infiltration, (4.2) how soil moisture is extracted from the soil profile by vegetation, (4.3) the thermodynamics and water transport within the soil, (4.4) the hydraulic and (4.5) thermal characteristics of the soil, (4.6) the treatment for preventing a soil layer from becoming super-saturated, and finally (4.7) the representation of heterogeneity for soil moisture.This is done via two possible methods, the first (4.7.1) being based upon the TOP-MODEL approach (Beven and Kirkby, 1979), and the second (4.7.2) the PDM model (Moore, 1985).

Surface fluxes and energy balance
The surface fluxes of heat, moisture and momentum are calculated in JULES within the surface exchange module.To give the maximum flexibility in terms of the representation of surface heterogeneity and for the coupling of the land surface scheme to an atmospheric model, two generic types of surface are considered; vegetated and non-vegetated.The main difference between these two types of surface is the way in which the surface related parameters (e.g., albedo, roughness length) are specified.For non-vegetative surfaces they are specified by the user (with the exception of the MORUSES option for an urban surface, see Sect.2.4), whereas for vegetated surfaces these parameters are derived from the structure of the vegetation itself.This leads to an alternative set of parameters that needs to be specified (e.g., rate of change of surface albedo with leaf area index, rate of change of roughness length with canopy height).

Surface exchange equations
The standard surface energy balance equations, used to calculate the distribution of available energy between the various fluxes at the surface, have been extended to provide more flexibility to include additional physical processes.Thermal inertia is associated with the surface mass which is coupled to the underlying soil by three physical mechanisms depending upon the type of surface.The vegetation fraction is coupled to the soil using radiative exchange and atmospheric turbulence, whereas the remainder are coupled through conduction.The surface energy balance equation is then written:  where: The definitions for all symbols are given in Appendix A, along with their units.
For the longwave radiative exchange between vegetation and the soil, one reflection of the emitted radiation is modelled (hence the reason why both emissivities appear in Eq. 4).This assumes that further reflections can be neglected.
A number of options can be chosen to adjust the formulation of the surface energy balance equations.These options increase the level of complexity for the interaction between the surface and the underlying soil, but have the capability to give improvements to the representation of the surface exchange of fluxes and the surface temperature, especially at night (Best and Hopwood, 2001).The traditional surface energy balance equations can be obtained by setting the surface heat capacity to zero (i.e., setting the left-hand side of Eq. (1) to zero) and having only conductive coupling between the surface and the underlying soil (i.e., by setting the vegetation fraction variable to zero in Eq. 4).This was the original surface energy balance that was used within the MOSES model, but Best and Hopwood (2001) showed that this did not provide sufficient cooling during the night over a grass surface.A better fit to the data was given if the surface is radiativly coupled to the underlying soil rather than coupling through conduction.These improvements are provided by the second option which uses not only radiative coupling, but also turbulence between the canopy and the underlying soil for vegetation surfaces, but still retains a zero surface heat capacity (C s = 0).A third option utilises the full energy balance equations above (Eqs.1-4).This introduces a heat capacity for the surface, which not only gives further improvements for tall vegetation such as forests that have a larger heat capacity than the grass surface considered in (Best and Hopwood, 2001), but also enables other surfaces (such as urban, see section 2.4) to be easily introduced within the model framework.The surface heat capacity is specified for non-vegetation surfaces, but is determined from the leaf and woody biomass for vegetation using Larger heat capacities result in a stronger thermal inertia for the surface.
Geosci.Model Dev., 4, 677-699, 2011 www.geosci-model-dev.net/4/677/2011/Fig. 2. The two dimensional canyon geometry used in MORUSES illustrating the resistance network used in the parametrisation of the roughness length for heat.The panels depict: (a) a wide canyon geometry with both ventilation and recirculation regions; and (b) a narrow canyon with only a recirculation region (adapted from Fig. 3 of Harman et al., 2004) .
In addition to utilising the full energy balance equations, there is a fourth option which adjusts how snow is represented on vegetation by enabling the snow to exist below the canopy (see Sect. 3.1).
In order to obtain a fully implicit solution, each of the prognostic terms in the surface flux equations (apart from the soil temperature) are written in the form X i+1 = X i + X.The equations are then linearised by assuming that X X.This gives a new set of surface flux equations that can be written in the form of a fully explicit flux, an update to give an implicit solution and a further update to ensure that the atmospheric temperature and humidity satisfy implicit coupling with the atmosphere.The last update is only applied if JULES is connected to an atmospheric model with implicit coupling.So, for example, the surface moisture flux equation becomes: where α is evaluated at T i * .The implicit update to the fluxes comes from solving the surface flux equations, whilst the implicit coupling to the atmosphere comes from the coupling methodology of Best et al. (2004).
The aerodynamic resistance is calculated using standard Monin-Obukhov similarity theory (Monin and Obukhov, 1954), using the stability functions of Dyer (1974) for unstable conditions and Beljaars and Holtslag (1991) for stable conditions.The surface resistance for surfaces with potential evaporation (i.e., lake, snow and ice surfaces) is set to zero, whilst for an urban surface the conductance is set to zero unless water is available on the urban surface (i.e., the urban "canopy water").For a bare soil surface, the surface conductance (g soil , inverse of resistance) is determined by the soil moisture concentration in the top soil layer: This parametrisation was developed following problems identified with a previous scheme (Taylor and Clark, 2001).A review of bare soil evaporation (including Mahfouf and Noilhan, 1991)  For vegetation, the surface resistance is calculated using the photosynthesis model described in Sect.2.2.
For the vegetative surfaces, the latent heat flux is determined from a combination of evapotranspiration and bare soil evaporation.The relative contributions from vegetation and bare soil are a representation of the fraction of bare soil that can be seen through the vegetation canopy.Hence the fractions for each of these is determined by the density of the leaves, through the leaf area index.The combined flux represents the interaction of the atmosphere with both the canopy and the soil beneath.
Note this is different to the approach used to represent the evaporation from a sparse canopy.In this situation, due to the limitations of the tile scheme approach as used in JULES (see Sect. 2.5), the surface is distributed into a vegetation land fraction that contributes to a vegetation tile, and a bare soil land fraction that contributes to the bare soil tile.

Photosynthesis and stomatal conductance
The leaf level stomatal conductance (g s ) and net photosynthetic uptake (A) are linked via the CO 2 diffusion equation: A second equation by Jacobs (1994), which shares similarities with the simplified form of the Leuning (1995) stomatal conductance formulation, relates the ratio of internal to external CO 2 concentrations to leaf humidity deficit, where f o and D * are vegetation specific calibration parameters, which are directly related to the parameters from the Leuning (1995) model (for details, see Cox et al., 1998).This simplified formulation is convenient for large scale model applications (Cox et al., 1998).Potential (non-water stressed) leaf level photosynthesis (A P ) is calculated in JULES using the C3 and C4 photosynthesis models of Collatz et al. (1991) and Collatz et al. (1992) respectively.Photosynthesis is simulated as the minimum of three limiting rates: (i) Rubisco limited rate (W C ), (ii) light limited rate (W L ) and (iii) rate of transport of photosynthetic products (in the case of C3 plants) and PEP-Carboxylase limitation (in the case of C4 plants) W E .With both, W C and W L having a dependency on the leaf internal CO 2 concentration, C i .
Leaf photosynthesis A, is related to the potential (nonstressed) leaf photosynthesis (A P ) as follows, β is the dimensionless moisture stress factor, which is related to the mean soil moisture concentration in the root zone, and the critical and wilting point concentrations as follows: The critical point is defined by a matrix water potential of −33 kPa (Cox et al., 1999), which compares to the more commonly used field capacity that has a matrix water potential of −10 kPa.The use of the critical point enables vegetation to maintain an un-water stressed transpiration at values below field capacity.JULES uses either a big leaf or a multi-layer approach to scale photosynthesis and conductance to the canopy level.In the big leaf approach, canopy level photosynthesis and conductance are calculated using leaf level fluxes and total canopy leaf area index (Cox et al., 1998) using Beer's law (Monsi and Saeki, 1953).This is the original method used in JULES, but does not produce a realistic dirunal cylce of photosynthesis and hence evaporation (Mercado et al., 2007(Mercado et al., , 2009)).A more realistic scheme is provided by the multilayer approach, in which the radiation absorbed and photosynthesis are estimated using a user defined number of leaf area increments (canopy layers) within the canopy, with the total canopy level flux calculated as the sum of the fluxes from each individual canopy layer (Jogireddy et al., 2006;Mercado et al., 2007).A number of options are available in JULES for use with this multilayer approach.In addition to the user specifying the number of layers, a two layer approximation can also be selected.This option is not as accurate as the full multilayer scheme, but saves on computational time which can be important for weather forecasting applications.Another option also allows for the variation of leaf nitrogen within the vegetation canopy, leading to further improvements within the multilayer scheme.Equations describing the biochemistry of leaf level photosynthesis (W C , W L and W E ) and scaling up methods from leaf to canopy level are outlined in Part II, which describes the carbon cycle in JULES (Clark et al., 2011).

Freely evaporating surfaces
Evaporation from the surfaces represented within JULES comes from a number of sources.These include evapotranspiration (i.e., water extracted from the soil through vegetation) and bare soil evaporation, both of which include a surface resistance that represents the restrictions in availability of water at the surface.The other sources of evaporation come directly from moisture stores and hence have no surface resistance.These sources include evaporation from open water surfaces, evaporation from surface water held in the canopy of vegetation or ponding on urban surfaces, and sublimation from snow.
The evaporation from water held on the leaves within the vegetation canopy will deplete the canopy water store and can result in all of the water being removed within a timestep.If this occurs, then the moisture unlimited evaporation is set to the available canopy water, and any additional evaporation then comes through evapotranspiration with an associated stomatal, or surface, resistance.Such a limitation in the evaporative flux changes the surface energy balance equations, so an adjustment is made to each of the terms in the energy balance equations to ensure that the model has a consistent solution.
Each surface type within JULES can have snow on it.When snow is present, the surface resistance is set to zero to represent the fact that there is a moisture source.Within JULES there is also an option to have the snow lying underneath vegetation for the turbulent moisture flux (Sect.3.1).In this case, an additional aerodynamic resistance is added to represent the efficiency of the turbulence at transporting moisture through the canopy.Any sublimation that occurs from the snow on the surface is used to deplete the snow mass in an analogous way to the canopy water.Also like the canopy water, if the snow is removed within a timestep, then an adjustment is made to the terms in the surface energy balance equations to ensure consistency.
Within JULES, lakes can be represented in two ways through the choice of available parameters.The default setting represents lakes as a bare soil surface, except that the surface resistance for the turbulent moisture flux is set to zero, giving a freely evaporating surface.The second method makes use of the surface canopy in the energy balance equations by setting a suitably large value for the surface heat capacity (typically equivalent to water of a depth of around 1 metre, although this can be altered by the user).This option reduces the diurnal cycle of the lake surface temperature compared to the first option, giving a more realistic simulation.
For both methods, as the lake is not explicitly modelled, the evaporative flux is not removed from any moisture store within the model, since it is assumed that there is sufficient water within the lakes to ensure that they are maintained.Similarly, any precipitation that falls onto the lake surface does not contribute to any water store.This means that in order to maintain a water balance, the integrated evaporative flux from the lake surface must be determined and included in the balance equations.This is not routinely done within JULES and has to be calculated through the available diagnostics by the user.
Similarly, the permanent ice surface does not have a prognostic water store, and hence care is required to maintain water balance.To represent an ice surface in JULES, the soil temperature profile is adopted to represent the thermal structure of the ice, whilst the moisture transport used in the soil scheme is neglected.As the ice surface is taken to be one of the surface tile types, and all surface tiles share the same soil information for temperature and moisture in JULES (see Sect. 2.5), this means that it is not possible to have a fractional coverage of land ice within a gridbox or source area at present.As such, there has to be either 100 % of land ice cover or none.The specification of this fraction of land ice is therefore done through the tile fractions information.
As with snow cover (Sect.3), the surface temperature of the ice surface is prevented from rising above the melting point of water, with any resulting residual of the surface energy balance being added to the melt flux.This means that care must be taken when setting land ice within the JULES model, especially when coupled to an atmospheric model.Small areas of ice could result in large horizontal thermal gradients in the atmosphere, caused by this restriction on the surface temperature compared to ice-free land.This can result in unrealistic small scale circulations and ultimately numerical problems.Hence when coupled to an atmospheric model, this surface type should only be used to represent a large extent of permanent land ice.

Representation of urban areas
The nature and design of urban environments make their surface energy balance significantly different from natural surfaces.However, a simple bulk representation for an urban area can be obtained by introducing a suitably large thermal capacity for the surface, along with radiative coupling between the surface and the underlying soil.Best (2005) showed that such a simple representation can lead to significant improvements within numerical weather prediction models.The advantage of this approach is that it is easy to adopt within a tile scheme approach and can fit within JULES by adapting currently available parameters.
A second option to represent urban areas in JULES, is to use an additional surface tile.Best et al. (2006) showed that representing the roofs of buildings as one surface and street canyons as a second effective surface gives improvements over the one-tile approach.Also, Harman and Belcher (2006) and Porson et al. (2009) demonstrated that these two surfaces give a good approximation of more complex schemes that represent each of the facets within the urban area.The differences between the two surface types is given through the surface parameter specifications.
The third option implemented is the Met Office Reading Urban Surface Exchange Scheme (MORUSES), as described in Porson et al. (2010a,b).Again this is a two-tile scheme, but as the surface parameters are determined from the morphology and material properties of the city, this enables a distribution of surface fluxes with different structural properties.The radiative exchange within the canyon tile is formulated with an effective albedo and an effective emissivity, based upon the exchanges between the various street canyon facets.The roughness length for momentum for the urban area is determined from the formulation of Macdonald et al. (1998) The roughness length for temperature comes from a physically-based parametrisation that relates to the urban morphology and uses a resistance network to represent the transfer of heat (see Fig. 2).The canyon tile includes the effects of the recirculation jets by using two resistance pathways; one for each of the recirculation and ventilation regions.For both of these elements, three resistances are used, two representing the heat across an internal boundary layer adjacent to each facet and one representing the transfer of heat across the inertial sub-layer.The roof, which is simpler, only has two resistances representing the internal boundary layer and inertial sub-layer (see Harman et al., 2004;Porson et al., 2010a, for more details).
Effective areal heat capacities are determined to represent the roof and the canyon, which includes contributions from both the walls and the road.These are determined by considering the diurnal response using a force-restore model, whilst an adjustable roof parameter is also introduced to increase the flexibility to capture different oscillations.The canyon tile is conductively coupled through the road to the underlying soil surface, whilst the walls of the canyon and the roof tile are decoupled from the soil, by imposing a zero flux boundary condition.
The data MORUSES requires can be sourced from a variety of places, depending on availability.For example MORUSES has been used to simulate the London urban heat island as part of the LUCID project (The Development of a Local Urban Climate Model and its Application to the Intelligent Design of Cities; Bohnenstengel et al., 2011) in which the Virtual London model (Evans et al., 2005) was used for building geometry.As part of the LUCID project, empirical relationships were also formulated for urban geometry to represent areas within the study area that did not exist within the Virtual London domain.The building material properties used were typical values: clay roof, brick walls and asphalt road.Where no information of this kind is known by the user, a global dataset also exists that categorises urban areas depending on density, climatic conditions and regional culture (Jackson et al., 2010).However, the amount and quality of the data known by the user would ultimately govern the choice of urban model used.

Surface heterogeneity
The heterogeneity of the surface is modelled within JULES by using the tile, or mosaic, approach (e.g., Essery et al., 2003a).This means that a separate surface energy balance is determined for each type of surface within the domain of the gridbox or footprint, and the individual surface fluxes are then given a weighted average in order to determine the gridbox or footprint mean flux into the atmosphere.One limitation to the current structure of JULES is that although the surface exchange represents the heterogeneity through tiling, there is no representation of sub-gridscale heterogeneity within the sub-surface soil module.This will be developed in future versions of the JULES model.
In order to keep the parametrisation of surface heterogeneity as flexible as possible, the number of surface types to be considered within a model simulation is determined at run time.Hence the complexity of the heterogeneity and cost in terms of computational time have to be balanced.Thus a time-limited modelling application, such as operational weather forecasting, can run with minimal surface types to optimize cost, whereas other applications may benefit from unlimited surface types (e.g., climate applications with an interactive carbon cycle).
There are two generic types of surface in JULES having differing requirements for their surface parameters: (1) Nonvegetated surfaces with fixed parameter values (e.g., albedo and roughness length) which are specified at run time, and (2) vegetated surfaces whose parameters vary.The latter are described in the following paragraphs.
The roughness length for momentum for vegetation is determined from There are two options to determine the surface albedo (α) for vegetation.The simplest option is a bulk albedo: where α b , the soil albedo, is a spatially varying ancillary field within JULES and α ∞ is the prescribed maximum canopy albedo for dense leaf coverage.With the second option, the snow-free albedos are calculated using the two-stream model for radiative transfer through vegetation described by Sellers (1985).This scheme Geosci.Model Dev., 4, 677-699, 2011 www.geosci-model-dev.net/4/677/2011/uses separate direct-beam and diffuse albedos in the visible and near-infrared wave bands for each vegetation type.This requires four parameter values for leaf reflection coefficients and leaf scattering coefficients for both near infra-red and photosynthetically active radiation.An additional parameter for vegetation surfaces is the capacity of the canopy to hold water (C m ) through the interception of precipitation, By default nine surface types are represented; five vegetation (broadleaf trees, needleleaf trees, C3 grasses, C4 grasses and shrubs) and four non-vegetated surfaces (urban, open water, bare soil and permanent land ice).The default parameters for each of these surface types are given in Tables 5 and  6, but where possible these parameters should be calibrated for specific sites.
In addition to the surface type, each tile has an elevation above the mean gridbox height.This enables surfaces that are sensitive to the changes in atmospheric temperature and humidity, arising from displacement above the mean surface height, to experience adjusted atmospheric forcing.This is done in a simple way by adjusting the air temperature along a dry adiabat whilst keeping the specific humidity constant until the saturation point is reached.After this, the temperature is adjusted along a moist adiabat, whilst the specific humidity is then set to the saturated specific humidity at the new atmospheric temperature.To ensure consistency with the top soil level temperature, this is adjusted by the same increment as the air temperature during the calculation of the surface energy balance (note that the actual prognostic soil temperature variable is not updated by this increment).This prevents artificial warming from the soil without having to introduce heterogeneity into the soil.This assumption will be removed once soil heterogeneity is introduced into the JULES code.One impact of introducing elevation bands is to reduce spurious sublimation and melting from snow-covered surfaces.

Snow model
Two schemes are available within JULES for the representation of snow on the ground.The simplest is a zero-layer scheme that uses no explicit model layers to represent snow, instead adapting the top soil level to represent lying snow processes.The more comprehensive and physically realistic scheme takes a multi-layer approach.For vegetated surfaces, snow may additionally be partitioned between intercepted snow in the canopy and snow on the ground or held in a single effective store.The simple and multi-layer snow scheme give similar results in many conditions, but the multilayer scheme is expected to give better simulations of snow dynamics at sites with deep snow, with the possibility of midwinter melt events and better simulations of soil temperatures at sites with low winter air temperatures.

Interaction of snow with vegetation canopies
With the original scheme in MOSES, snow is held in a single store and hence sits on top of vegetation regardless of the type and height of this vegetation.The exception of this is for the albedo which does account for a darker surface when snow is under tall vegetation.However, PILPS-2e (Bowling et al., 2003) found that the models with highest winter sublimation had lowest annual runoff for a high-latitude basin.A reformulation of MOSES to distinguish between snow on and below forest canopies reduced the sublimation and improved the simulation of runoff (Essery and Clark, 2003).So an option is also available in JULES to partition the snow between snow on the canopy and the underlying ground (Essery et al., 2003a).The surface resistance for sublimation is set to zero for tiles with snow cover in the single-store option, but is for canopy snow, where I max = 4.4L is the snow interception capacity for a canopy with leaf area index L and r = 0.5 mm is a nominal grain radius for intercepted snow (Essery et al., 2003b).The change in load during a timestep with snowfall amount S f on a canopy with initial load I 0 is Snow is removed from the canopy by sublimation, and unloading of melting snow from the canopy is set equal to 40 % of the canopy snowmelt rate (Storck et al., 2002;Essery et al., 2003b).

Zero-layer snow model
The original snow scheme within JULES is a zero-layer snow model.Snow is given a constant thermal conductivity and a constant density.The heat capacity of snow is neglected, but snow decreases the bulk thermal conductivity of the surface layer due to both the increased layer thickness and the different conductivities of snow and soil.For snow depth less than half the surface soil layer thickness ( z 1 ), the thermal conductivity used in surface energy balance calculations is adjusted for insulation by snow according to The heat flux between the surface layer and the second soil layer, of thickness z 2 , is multiplied by a snow insulation factor www.geosci-model-dev.net/4/677/2011/Geosci.Model Dev., 4, 677-699, 2011  (Moore, 1985) For deeper snow, the surface conductivity is set equal to λ snow and the insulation factor is (20) (Cox et al., 1999).The surface skin temperature is not allowed to exceed 0 • C while snow remains on the ground, and the heat flux used to melt snow is diagnosed as a residual in the surface energy balance.Melt water drains immediately from the snow and is partitioned into soil infiltration and runoff; there is no storage or freezing of liquid water in snow.The snow thermal conductivity, snow density and surface layer thickness are parameters that are set by the user.Whilst the zero-layer snow scheme on the whole gives good agreement with observations, it tends to melt snow too rapidly.This is partly due to the inability of this scheme to hold liquid water within the snow that can subsequently re-freeze.In addition, the use of the top soil layer to represent the snow has a negative impact on the soil temperatures, e.g., as demonstrated in the SNOWMIP2 experiment (Essery et al., 2009).Better agreement between observations and the JULES model can be obtained by using the alternative multilayer snow scheme.

Multi-layer snow model
The maximum number of layers (N max ) that are used for deep snow and their thickness d k (k=1,...,N max ) are set by the user.
However, the number of layers actually used depends on the snow depth, which means that not all the layers exist at any one time.When a layer is at the base of the snowpack it has a variable thickness.Shallow snow is combined with the surface soil layer for snow depth d s < d 1 for numerical stability, whilst setting N max = 0 forces the use of the zero-layer option for any depth of snow.For d s ≥ d 1 , snow is represented by additional model layers on top of the soil if N max ≥ 1.As the snow depth increases, the lowest layer in the snowpack increases in thickness until it reaches twice its prescribed thickness; the layer then splits in two with the upper part staying fixed in thickness and the new lowest layer thickening as the snow accumulates.This is reversed as the snow depth decreases, with layers being progressively combined at the bottom of the snowpack.The division of a snowpack into layers is illustrated in Fig. 3.A variable snow density is used, so snow depth can decrease due to compaction as well as ablation.
Each layer in the snowpack has a thickness d k (m), a temperature T k (K), a density ρ k (kg m −3 ), an ice content I k (kg m −2 ) and a liquid water content W k (kg m −2 ).Layer thickness, density and mass are related by ρ k d k = I k + W k .The increase in layer density due to compaction over a timestep of length δt is calculated as Geosci.Model Dev., 4, 677-699, 2011 www.geosci-model-dev.net/4/677/2011/where k s = 4000 K, compactive viscosity η 0 = 10 7 Pa s, reference density ρ 0 = 50 kg m −3 , temperature T m = 273.15K, and M k = 0.5(I k +W k )+ k−1 i=1 (I i +W i ) is the mass of snow above the middle of the layer.This scheme, based on measurements by Kojima (1967), has previously been used in the snow models described by Pitman et al. (1991) and Lynch-Stieglitz (1994).The areal heat capacity of a layer is where C ice = 2100 J K −1 kg −1 and C water = 4180 J K −1 kg −1 are the specific heat capacities of ice and water, and the thermal conductivity is λ k = 2.22 ρ k ρ water 1.88 (23) where ρ water = 1000 kg m −3 is the density of water (Yen, 1981).
The structure of the multi-layer snow model is shown in Fig. 4. The conducted heat flux at the bottom of layer k is where δT k is the increment in layer temperature over a timestep, γ is the forward timestep weighting (0 for explicit and 1 for fully implicit timestepping), and is a layer thickness weighted thermal conductivity.For the lowest snow layer (k = N), T N +1 , d N+1 and λ N+1 are the temperature, thickness and conductivity of the surface soil layer.The increment in layer temperature over a timestep is Surface heat flux H 0 calculated by the surface exchange module is passed to the snow module, and ground heat flux H N calculated by the snow module is passed to the soil module; implicit timestep weighting of surface soil layer temperature T s1 is not used in calculating this flux.For a single snow layer the temperature increment is given by with solution When there are N > 1 snow layers, increments in the layer temperatures are found as the solutions of the tridiagonal set of equations for k = 2,...,N − 1, and with matrix elements and If the temperature of a layer is calculated to be above T m , the layer ice mass is reduced by an amount or the entire mass of the layer, whichever is least.The layer liquid mass is increased by the same amount and the layer www.geosci-model-dev.net/4/677/2011/Geosci.Model Dev., 4, 677-699, 2011 temperature is reset to T m .Sublimation calculated by the surface exchange module is removed from the surface layer ice mass and from deeper layers if the surface layer sublimates entirely during a timestep.A layer of depth d k entirely consisting of liquid water would have a liquid content of ρ water d k .Snow layers are allowed to retain a fraction W cap (set by the user) of this liquid content.When the liquid content of a layer exceeds its capacity, excess water is passed to the layer below.Liquid water in a layer with temperature below T m will freeze, decreasing the liquid content by an amount increasing the ice content by the same amount, and increasing the temperature by The water flux at the base of the snowpack is passed to the surface hydrology module (Sect.4.1).Fresh snow is added as an interim layer 0 with density ρ 0 and temperature equal to the surface layer temperature.After increments have been applied to the layer masses and temperatures, layers are combined or split as necessary to match the fixed layer thicknesses.The liquid water and ice contents of the revised snow layers are determined by conservation of mass and their temperatures are diagnosed from conservation of energy

Snow albedo
Diagnostic and prognostic snow albedo options are provided.The simpler diagnostic option was originally used in MOSES, but can not represent the impacts of snow ageing on the surface albedo.Hence the prognostic option provides the ability to represent the time evolution of the snow abledo, improving the physical representation of the snow.
In the diagnostic scheme, a snow-free albedo and an albedo for cold deep snow are specified for each surface type.When the surface temperature exceeds a threshold temperature T c , the snow albedo is decreased according to For a tile with snow depth d s , the albedo is a weighted average for surface masking snow depth d m .
For tall vegetation, the impact of snow lying underneath the vegetation canopy is taken into account by setting lower values for the cold deep snow albedo.
The prognostic albedo scheme uses the Wiscombe and Warren (1980) spectral snow model.The ageing of the snow surface is characterized by introducing a prognostic grain size r(t), set to r 0 = 50 µm for fresh snow and limited to a maximum value of 2000 µm.The change in grain size over a timestep is given by for snowfall rate S f .The mass of fresh snow required to refresh the albedo is set to 2.5 kg m −2 .The empirical grain area growth rate, in µm 2 s −1 , is where A s = 0.23 × 10 6 µm 2 s −1 .Snow albedos for diffuse visible and near-infrared radiation are calculated as and The zenith-angle dependence of albedos for direct-beam radiation with zenith cosine µ is represented by using an effective grain size, Geosci.Model Dev., 4, 677-699, 2011 www.geosci-model-dev.net/4/677/2011/ in place of r in the equations for diffuse albedos (Eqs.41 and 42).
For a tile with snow-free albedo α 0 , snowdepth d s and roughness length z 0 , the albedo in each band is where (45) (Oleson et al., 2004).When driving data with separate directbeam and diffuse radiation in visible and near-infrared bands are not available, the average of the diffuse albedos is simply used as an all-band snow albedo.

Hydrology and soil thermodynamics
JULES includes multi-layer, finite-difference models of subsurface heat and water fluxes, as described in Cox et al. (1999).There are options for the specification of the hydraulic and thermal characteristics, the representation of super-saturated soil moisture and the sub-surface heterogeneity of soil moisture.

Surface hydrology
To account for the size of convective storms compared to gridsize, a rainfall rate is assumed to fall on a fraction r of the grid.For large scale precipitation and point studies this fraction is set to one, whilst for convective precipitation it can take lower values, and is typically set to a value of 0.3.The amount of water that reaches the soil surface depends upon the type of surface.For non-vegetation surfaces, this is simply the precipitation rate whereas for vegetation surfaces, this becomes the throughfall and is calculated using: and the canopy water is updated by where C m is the maximum canopy water that can be held by the vegetation and t is the timestep.The canopy water can also be increased through dewfall (i.e., downward surface moisture fluxes), and is depleted by surface evaporation.Similarly, snow cover is increased through the deposition of frost (modelled as dewfall at surface temperatures below freezing), whilst the melting of snow contributes to the water available at the soil surface and updates the equivalent water within the snow pack.
The water reaching the soil surface is then split between infiltration into the soil and surface runoff.Since the throughfall can be different for various surface types within a tile scheme, whereas in JULES the infiltration into the soil is a grid-box aggregate as the soil is not tiled, the surface runoff is determined by combining equations for both the throughfall and the grid-box mean infiltration, yielding the following: where the surface infiltration rate K is equal to β s K hs and β s is an enhancement factor.These equations account for the effect of a finite model timestep on the throughfall and therefore the surface runoff.A full derivation of these equations was given by Dolman and Gregory (1992).
The infiltration into the soil is determined through the integration of the contributions for each of the surface types by using the water balance at the surface: (49)

Soil moisture extraction
The ability of vegetation to access moisture at each level in the soil is determined by root density, assumed to follow an exponential distribution with depth.The fraction of roots in soil layer k extending from depth z k−1 to z k is For transpiration E , the flux extracted from soil layer k is e 0 k E , where and is a soil moisture availability factor, defined similarly to Eq. ( 12), for a soil layer with unfrozen volumetric soil moisture concentration θ k .

Soil thermodynamics and water fluxes
The sub-surface at a gridpoint is either soil or land ice (with no water movement in the latter).Sub-surface temperatures are calculated using a finite difference form of the heat diffusion equation, including the effects of solid-liquid phase changes of water.The temperature of the k-th soil layer is incremented by the diffusive heat fluxes into and out of the layer, G k−1 and G k , and the advective flux from the layer by flowing water, J k : www.geosci-model-dev.net/4/677/2011/Geosci.Model Dev., 4, 677-699, 2011 where the fluxes are calculated as where z is the vertical coordinate.C a is the "apparent" volumetric heat capacity of the layer, including the effect of phase changes (Cox et al., 1999).For soil, the sub-surface thermal characteristics are a function of solid and liquid water contents, while land ice uses fixed characteristics.The top boundary condition for Eq. ( 53) is the surface heat flux, calculated by the surface exchange module, while at the bottom there is a zero flux boundary condition to ensure conservation of energy.
The number of soil layers is a model parameter but the default is four of thicknesses 0.1, 0.25, 0.65 and 2.0 m, giving a total soil depth of 3 m.This configuration is designed to capture the variation of soil temperature from sub-daily to annual timescales (Best et al., 2005).
Soil water contents are updated using a finite difference form of the Richards equation.The moisture content of each layer is updated as: where W k−1 and W k are the diffusive fluxes flowing in from the layer above and out to the layer below respectively, E k is the evapotranspiration extracted by plant roots in the layer (and bare soil evaporation for the top layer) and R bk is lateral runoff, which is set to zero unless the sub-surface heterogeneity of soil moisture is represented using the TOP-MODEL option (Sect.4.7.1).The vertical fluxes follow Darcy's law The top boundary condition for Eq. ( 56) is the infiltration of water at the soil surface, whilst the lower boundary condition is drainage, which contributes to sub-surface runoff.

Hydraulic characteristics
There are two options for the hydraulic characteristics.In the first the relation between soil water content, suction and hydraulic conductivity are Brooks and Corey (1964): This is the method that has traditionally been used in land surface models.It is a more simplistic formulae for the hydraulic properties of the soil than other schemes, but the required parameters can be determined from the sand, silt and clay fractions of soils, which are available in many global soil datasets.
The parameters θ s , s and b are calculated from soil texture information using the relationships of Cosby et al. (1984) or others.(Note that Cox et al. (1999) incorrectly referenced Eqs. ( 58) and ( 59) as Clapp and Hornberger (1978) rather than Brooks andCorey (1964) (T. Marthews, personal communication, 2009).) The second option uses the hydraulic relationships of van Genuchten (1980), which is a more complex formulae but more scientifically robust: where m = 1 − 1/n and S = (θ − θ r )/(θ s − θ r ).In JULES, ξ = 0.5 and the soil moisture variable is implicitly defined as θ − θ r , leaving three parameters.The specific parameters required for this formulation have not been traditionally available within soil dataset, making it difficult to use.However, more recent datasets now include these parameters within their soil information.Dharssi et al. (2009) show that with suitable parameter values, Eqs. ( 58) and ( 60) are similar over most of the soil moisture range.
The soil parameter values can vary between layers but, in the absence of suitable data with which to specify this variation, many applications ignore any variation with depth.When calculating the hydraulic characteristics using Eqs.(58-61), JULES uses θ u , the unfrozen volumetric water content, instead of θ, to capture the effects of soil freezing, following Cox et al. (1999).

Thermal characteristics
JULES has two options for calculating the effective thermal conductivity of soil λ.The first option (described by Cox et al., 1999) is a less complex scheme, but requires only a limited amount of soil information: where where and S u and S f are the unfrozen and frozen water contents as a fraction of saturation.Dharssi et al. (2009) showed that the thermal characteristics from this scheme do not agree well with those for various soil types.It was found that the formulation of Johansen (1975) gives a better fit, but this scheme requires additional soil information.Although this additional information is Geosci.Model Dev., 4, 677-699, 2011 www.geosci-model-dev.net/4/677/2011/generally available in the latest soil datasets, implementing it into the Met Office Unified Model requires significant effort.Hence Dharssi et al. (2009) derived a simplification of Johansen (1975) which gives similar response in the relationships between the thermal conductivity and soil moisture: where K e is the Kersten number with the constraint that 1•58 ≤ λ u s ≤ 2•2.This equation for the thermal conductivity of unfrozen saturated soil (λ u s ) was derived in order to give good agreement with the Johansen (1975) formulation, but without requiring knowledge of the mineral content of the soil (Dharssi et al., 2009).
This generally gives larger values for conductivity than the Cox et al. (1999) formulation, which reduces the errors in simulated air temperature when used in Numerical Weather Prediction (Dharssi et al., 2009).

Super-saturation of soil moisture
The numerical solution for the transport of soil moisture between the soil layers may result in layers which become super-saturated.JULES has two options to prevent this from occurring.With the first option, if a soil layer becomes supersaturated, then the soil moisture in this layer is limited to the saturation point and the excess water is prevented from moving into the layer from above, i.e., the drainage into the layer is restricted by the saturation.This results in the excess water being moved back up the soil layer, and if the top soil layer becomes super-saturated, then the excess water is added to the surface runoff.
The second option is to route the excess soil moisture to the soil layers below.This assumes excess soil moisture might flow laterally over land within a large gridbox, but would eventually move down through the soil layers at subgrid locations in which drainage is less impeded (e.g., where there is fractured permafrost or less compacted/faster draining soil types).This results in the excess water being moved down to lower layers, and if the bottom soil layer becomes super-saturated, then the excess water is added to the subsurface runoff.
If the total soil column is saturated, then the difference between these two options is to add the excess water to either the surface or sub-surface runoff.Whilst in both cases the water results in a runoff flux, this could impact the timing of river flow due to the delay of sub-surface runoff getting into the river network.Tests of the two options with the PILPS2d Valdai data (Schlosser et al., 2000) showed that moving the excess water in the downwards direction led to a poor surface runoff simulation and excessive soil moisture, whereas inhibiting the drainage of water into a saturated layer gave a better agreement with observations.
However, global simulations have shown that in regions of partially frozen soils, one possible result is saturated and partially frozen soil layers near the surface, with unsaturated layers below.In this situation, the option to inhibit the drainage of water into the saturated layer at the surface leads to excessive surface runoff of snowmelt, giving a dry soil during spring and hence a dry and warm bias in the atmosphere during the summer.The option to move the excess water to lower layers moistens the lower unsaturated soil layers and removes some of this dry and warm atmospheric bias whilst reducing the surface runoff of snowmelt.
These results suggest that the grid size may be important in determining the dominant physical processes that prevent the super-saturation of the soil, and further work is required to determine how this should be represented in the model.Thus care should be taken when choosing between options for controlling super-saturation, with consideration being taken for the required application.

Soil moisture heterogeneity
There are two options in JULES to introduce sub-gridscale heterogeneity into the soil moisture.One (TOPMODEL) is a more complex scheme that represents this heterogeneity throughout the soil column, including aspects such as a water table and the capability to estimate wetland fractions (Gedney et al., 2004).However, this scheme requires additional topographic ancillary information.The other option (PDM) is a much simpler scheme that does not require as much information.This scheme only considers heterogeneity in the top soil layer and thus can not be used to represent the water table depth or to determine wetland areas.However, it can still be used to increase surface runoff and has been shown to improve subsequent river discharge when fed into a river routing scheme.

TOPMODEL
JULES can optionally use a parameterisation based on TOP-MODEL (Beven and Kirkby, 1979).TOPMODEL was initially designed to include a groundwater model within a single catchment where the height of the saturated zone moves up and down and is controlled by the recharge into it and the saturated lateral flow (baseflow) out.As the water table becomes higher, more of the surface area becomes saturated (and vice versa), with the regions of higher topographic index ( λi ) flooding first.Topographic index relates to the upstream area draining into a locality and the local slope, which www.geosci-model-dev.net/4/677/2011/Geosci.Model Dev., 4, 677-699, 2011 1.00 3 ×10 −4 3×10 −4 1 × 10 −4 Canopy heat capacity (J K −1 m −2 ) 0.28×10 6 0.00 0.00 0.00 Fractional "canopy" coverage 1.00 0.00 0.00 0.00 * The snow-free albedo for soil is initialised to −1 to allow it to be set through an ancillary field instead.
is a measure of the potential to flood relative to other regions within the catchment.This distributed catchment-based model is simplified into a semi-distributed model (Sivapalan et al., 1987) for use in climate models by lumping areas of similar topographic index together from one or potentially more catchments.A gridbox mean water table depth z w is calculated, and the probability distribution function (pdf) of the topographic index within the gridbox is then used to describe the relative frequency of occurrence of the topographic indices.The gridbox fraction of the water table that is above the surface may then be calculated.This enables saturation excess runoff to occur in the model before the gridbox soil moisture is totally saturated.Runoff occurs when water is unable to permeate the fraction of the gridbox surface where the water table is above the surface.
The implementation of this approach in JULES was adapted by Gedney and Cox (2003) and Clark and Gedney (2008).With the TOPMODEL-based approach the free drainage lower boundary condition is replaced by a no flux condition, and sub-surface runoff is represented as a lateral "baseflow", described below.An extra soil layer, with simplified representation of water fluxes, is added beneath the standard soil model as a computationally efficient way in which to track the water table when it is deeper than the standard 3 m soil column.JULES assumes an exponential decrease of K s in this deeper layer with a decay constant f .The lateral sub-surface runoff, or baseflow, is calculated as where T ( z w ) is the vertical transmissivity from the bottom of the column to the z w .This transmissivity is found by summing the contributions from each layer and only considers unfrozen soil water.The transmissivity of each layer is used to partition the total baseflow between soil layers, to give the layer values R bn required in Eq. ( 56).z w is calculated by assuming that the column soil moisture is in equilibrium (Koster et al., 2000).
The "critical" value of the topographic index at which the water table reaches the surface is found as λic = ln(R b max /R b ), where R b max is the baseflow found from Eq. 68 with z w =0.The fraction of the gridbox that is saturated (f sat ) can be found by integrating the pdf of the topographic index.However, this requires numerical integration if a two-parameter gamma distribution is used for the pdf as in Gedney and Cox (2003).Instead, during the initialisation an exponential distribution is fitted to the results of the gamma distribution, and subsequently f sat is found using where a s and c s are fitted parameters for each gridbox.Saturation excess surface runoff (R se ) is then calculated as where W 0 is the rate at which water arrives at the soil surface from precipitation and snowmelt (Eq.49).
The fraction of the gridbox that is considered to be wetland (i.e., stagnant water) for the purposes of methane emissions (f wet ) is defined as that part of the gridbox at which λ ic ≤ λ i ≤ λ i max where λ i max is a global parameter.At locations with larger values of λ i (water higher above surface) the water is assumed to be flowing and not wetland.Following the procedure for f sat , an exponential relationship is fitted so that f wet can subsequently be calculated as f wet = a wt exp(−c wt f λ ic ) for parameters a wt and c wt .Gedney and Cox (2003) and Clark and Gedney (2008) showed that simulated runoff was improved by using a TOPMODEL type parameterisation, and that the global pattern of wetland is captured by this model (Gedney and Cox, 2003).

Probability Distribution Model (PDM)
An alternative to TOPMODEL is to calculate saturation excess runoff following the Probability Distributed Model (PDM, Moore, 1985).The distribution of soil storage capacity within a gridbox is modelled by a pdf, and the saturated fraction of the gridbox can be shown to be Geosci.Model Dev., 4, 677-699, 2011 www.geosci-model-dev.net/4/677/2011/ where B is a shape parameter.R se is then calculated using Eq.70.In JULES, B is kept constant across the domain, as is the depth over which W and W max are calculated (typically 1 m).The calculations of infiltration excess and sub-surface runoff are not altered if PDM is selected.Clark and Gedney (2008) showed that the use of PDM improved modelled runoff in mesoscale catchments.

Summary
The Joint UK Land Environment Simulator (JULES) is a new community land surface model, based upon the established Met Office Surface Exchange Scheme (MOSES).In addition to representing the exchange of fluxes of heat and moisture between the land surface and the atmosphere (as described here), the model also represents fluxes of carbon and some other gases such as ozone and methane (described in Clark et al., 2011).This enables JULES to be used for many applications, and results in it being a unique land surface model in the fact that it is used in both operational weather forecasting and leading climate change simulations.Unlike many land surface models, JULES has an explicit representation of the surface energy balance for vegetation, capturing the weaker coupling that exists between the canopy and underlying soil.Other models (e.g., TESSEL, Viterbo and Beljaars, 1995) represent this weaker coupling by adjusting the thermal properties of the top soil layer, but do not have the flexibility of representing radiative, turbulent and conductive exchanges that can be represented in JULES.
Like most other land surface models, JULES uses a tiled land surface scheme to represent heterogeneity in land cover.Many land models have fixed descriptions of the surface types that are designed with specific applications in mind.However, the flexible structure within JULES enables the description of the resolved surface types to be targeted for specific applications.This means, for instance, that there can be a small number of vegetation types for weather forecasting applications where computation cost is critical, but many vegetation types for climate modelling where an accurate representation of the various biomes is important.In addition, JULES introduces elevation bands to the surface types, which is not common in land surface models.The elevated surfaces enable a modified surface energy balance which can be critical for the evolution of snowmelt and sublimation.
Another feature of the snow scheme within JULES is the ability to separate snow held in vegetation canopies and the snow under the canopy, although many other models also make this distinction (e.g., CLASS, CLM).This reduces the spuriously enhanced sublimation of the snow due to an incorrectly increased surface roughness from the tall vegetation components.The new multiple layer snow scheme within JULES also impacts on the timing of snowmelt through the introduction of both solid and liquid water stores.Other land models have a range in the number of snow layers that are modelled, for instance, CLASS uses one layer (Bartlett et al., 2006) whereas CLM uses up to five (Oleson et al., 2010), whilst ISBA has both an implicit snow layer (Douville et al., 1995) or a three-layer snow model (Boone and Etchevers, 2001).However, the majority of snow schemes include both solid and liquid stores within their layers.
JULES has a multilayer approach for both radiation interception and photosynthesis for vegetation.This has been shown to give an improved diurnal cycle for photosynthesis compared to the big leaf approach to scale from leaf to canopy level, that uses only beers law for light interception through the canopy, but is used in some models (such as LPJ, Sitch et al., 2003).Other models do use a multilayer canopy scheme for photosynthesis, but still use beer's law for light interception (e.g., LPJ GUESS, Smith et al., 2001).
There is a selection of three possible options for representing urban surfaces within JULES.All three options have been shown to give a good representation of sensible and latent heat fluxes over urban surface (e.g., in the first international urban model comparison experiment, Grimmond et al., 2010Grimmond et al., , 2011)).The urban surface is integrated into the general framework of the land model, unlike some other models that have to couple an urban model to a separate land model (e.g., ISBA and TEB, Noilhan and Planton, 1989;Masson, 2000).
The heterogeneity of soil moisture can be represented with two methods of varying complexity within JULES.The simple method represents the heterogeneity in the top soil layer only, but can generate increased surface runoff, whereas the more complex scheme has a representation of the mean water table depth.Whilst some land models include the more complex scheme, many do not include soil moisture heterogeneity at all, whilst few have the simpler more computationally efficient method.
The JULES model has been designed with a flexible and modular structure, which means that new elements of science can easily be introduced as new modules into the model.The scientific developments for each module are co-ordinated by an expert in the relevant area of science, ensuring that the model will remain a state of the art land surface model for the research community.

Fig. 1 .
Fig. 1.Modular structure of the JULES model.The boxes show each of the physics modules whilst the lines between the boxes show the physical processes that connect these modules.The surrounding three boxes show the cross-cutting themes.

Fig. 4 .
Fig. 4. Structure of the numerical discretisation over the layers for the temperatures and heat fluxes within the multi-level snow scheme in JULES.

Table 1 .
Meteorological forcing data required to drive the JULES model.

Table 2 .
Prognostic variables within the JULES model.
1 Only for the multi-layer snow option.2Onlyfor the TOPMODEL soil moisture heterogeneity option.

Table 3 .
Soil ancillary data required by the JULES model.
, for a staggered array of cubes; the canyon and the roof tiles both have the same roughness.Simulation of seasonal snow depth with JULES for multilevel snow scheme, showing the division into a varying number of layer depths.The full shaded area shows the total snow depth, whilst the different shadings represent the depths of the various snow layers.Minimum layer thicknesses can be selected by the user, but in this illustration a second layer is added when the snow depth exceeds 20 cm and a third at 30 cm.

Table 4 .
Description of the various physics options within the JULES model as discussed in the identified sections.

Table 5 .
Default parameter values required by JULES for the standard vegetation surfaces.

Table 6 .
Default parameter values required by JULES for standard non-vegetation surfaces.