# The Joint UK Land Environment Simulator (JULES), model description - Part 1: Energy and water fluxes

M. J. Best, M. Pryor, D. B. Clark, G. G. Rooney, R. L. H. Essery, C. B.

Ménard, J. M. Edwards, M. A. Hendry, A. Porson, N. Gedney, et al.

► To cite this version:

M. J. Best, M. Pryor, D. B. Clark, G. G. Rooney, R. L. H. Essery, et al. The Joint UK Land Environment Simulator (JULES), model description - Part 1: Energy and water fluxes. Geoscientific Model Development, 2011, 4, pp.677-699. 10.5194/gmd-4-677-2011. hal-04113888

HAL Id: hal-04113888

https://hal.science/hal-04113888

Submitted on 2 Jun 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# The Joint UK Land Environment Simulator (JULES), model description – Part 1: Energy and water fluxes

M. J. Best<fnref n="1" />, M. Pryor<fnref n="2" />, D. B. Clark<fnref n="3" />, G. G. Rooney<fnref n="1" />, R. L. H. Essery<fnref n="4" />, C. B. Ménard<fnref n="4" />, J. M. Edwards<fnref n="1" />, M. A. Hendry<fnref n="1" />, A. Porson<fnref n="1" />, N. Gedney<fnref n="2" />, L. M. Mercado<fnref n="3" />, S. Sitch<fnref n="5" />, E. Blyth<fnref n="3" />, O. Boucher<fnref n="1" />,*, P. M. Cox<fnref n="6" />, C. S. B. Grimmond<fnref n="7" />, and R. J. Harding<fnref n="3" />

Received: 1 March 2011 – Published in Geosci. Model Dev. Discuss.: 24 March 2011

Revised: 20 July 2011 – Accepted: 22 July 2011 – Published: 1 September 2011

**Abstract.** 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 Office 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 Office Unified 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-profile 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 flexible modelling platform.

# 1 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

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.

<figure><img src="image_3.png" /></figure>

Correspondence to: M. Best
(martin.best@metoffice.gov.uk)

<footnote n="1">¹Met Office, FitzRoy Road, Exeter, EX1 3PB, UK</footnote>

<footnote n="2">²Met Office, JCHMR, Wallingford, UK</footnote>

<footnote n="3">³Centre for Ecology and Hydrology, Wallingford, UK</footnote>

<footnote n="4">⁴University of Edinburgh, Edinburgh, UK</footnote>

<footnote n="5">⁵University of Leeds, Leeds, UK</footnote>

<footnote n="6">⁶University of Exeter, Exeter, UK</footnote>

<footnote n="7">⁷King's College London, London, UK</footnote>

<footnote>* now at: Laboratoire de Météorologie Dynamique, IPSL, CNRS/UPMC, Paris, France</footnote>
<table><caption>**Table 1.** Meteorological forcing data required to drive the JULES model.</caption><thead><tr><th>Data</th><th>Units</th></tr></thead><tbody><tr><td>Downward component of shortwave radiation at the surface</td><td>W m<sup>-2</sup></td></tr><tr><td>Downward component of longwave radiation at the surface</td><td>W m<sup>-2</sup></td></tr><tr><td>Rainfall</td><td>kg m<sup>-2</sup> s<sup>-1</sup></td></tr><tr><td>Snowfall</td><td>kg m<sup>-2</sup> s<sup>-1</sup></td></tr><tr><td>U component of wind</td><td>m s<sup>-1</sup></td></tr><tr><td>V component of wind</td><td>m s<sup>-1</sup></td></tr><tr><td>Atmospheric temperature</td><td>K</td></tr><tr><td>Atmospheric specific humidity</td><td>kg kg<sup>-1</sup></td></tr><tr><td>Surface Pressure</td><td>Pa</td></tr></tbody></table>

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., 2007, 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)
<table><caption>**Table 2.** Prognostic variables within the JULES model.</caption><thead><tr><th>Data</th><th>Units</th></tr></thead><tbody><tr><td>Fractions of land surface types within gridbox</td><td></td></tr><tr><td>Surface temperature of land surface types</td><td>K</td></tr><tr><td>Temperature of each soil level</td><td>K</td></tr><tr><td>Moisture concentration of each soil layer</td><td>m<sup>3</sup> m<sup>-3</sup></td></tr><tr><td>Canopy water for vegetation surface types</td><td>kg m<sup>-2</sup></td></tr><tr><td>Canopy height for vegetation surface types</td><td>m</td></tr><tr><td>Leaf area index for vegetation surface types</td><td>m<sup>2</sup> m<sup>-2</sup></td></tr><tr><td>Soil carbon</td><td>kg C m<sup>-2</sup></td></tr><tr><td>Snow amount of each surface type</td><td>kg m<sup>-2</sup></td></tr><tr><td>Snow grain size on each surface type</td><td>µm</td></tr><tr><td>Snow on ground below vegetation surface types</td><td>kg m<sup>-2</sup></td></tr><tr><td>Snowdepth of each surface type<sup>1</sup></td><td>m</td></tr><tr><td>Temperature of each snow layer for each surface type<sup>1</sup></td><td>K</td></tr><tr><td>Ice content of each snow layer for each surface type<sup>1</sup></td><td>kg m<sup>-2</sup></td></tr><tr><td>Liquid content of each snow layer for each surface type<sup>1</sup></td><td>kg m<sup>-2</sup></td></tr><tr><td>Snow grain size of each snow layer for each surface type<sup>1</sup></td><td>µm</td></tr><tr><td>Mean water table depth<sup>2</sup></td><td>m</td></tr><tr><td>Soil moisture in deep layer as a fraction of saturation<sup>2</sup></td><td></td></tr></tbody></table>

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 describes 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 TOPMODEL approach (Beven and Kirkby, 1979), and the second (4.7.2) the PDM model (Moore, 1985).

## 2 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).

### 2.1 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:

$$C_s \frac{\delta T_*}{\delta t} = (1-\alpha) S w_{\downarrow} + \epsilon L w_{\downarrow} - \sigma \epsilon (T_*)^4 - H - L_c E - G \quad (1)$$

<footnote n="1">¹ Only for the multi-layer snow option.</footnote>

<footnote n="2">² Only for the TOPMODEL soil moisture heterogeneity option.</footnote>
<figure><img src="image_4.png" /><figcaption>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.</figcaption></figure>

where:

$$H = \frac{\rho c_p}{r_a} (T_* - T_A) \qquad (2)$$

$$E = \frac{\rho}{r_a + r_s} (Q_{\text{sat}}(T_*) - Q_1) \qquad (3)$$

$$G = \nu \left[ \sigma \epsilon \epsilon_s (T_*)^4 - \sigma \epsilon \epsilon_s (T_{s1})^4 + \frac{\rho c_p}{r_{a\text{can}}} (T_* - T_{s1}) \right] + (1 - \nu) \lambda_{\text{soil}} (T_* - T_{s1}) \qquad (4)$$

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 radiatively 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

$$C_s = C_L B_L + C_W B_W \qquad (5)$$

Larger heat capacities result in a stronger thermal inertia for the surface.
<table><caption>**Table 3.** Soil ancillary data required by the JULES model.</caption><thead><tr><th>Data</th><th>Units</th></tr></thead><tbody><tr><td>Bare soil albedo</td><td></td></tr><tr><td>Dry soil thermal conductivity</td><td>W m<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>Dry soil thermal capacity</td><td>J K<sup>-1</sup> m<sup>-3</sup></td></tr><tr><td>Volumetric saturation point for soil</td><td>m<sup>3</sup> m<sup>-3</sup> of soil</td></tr><tr><td>Critical volumetric soil moisture content</td><td>m<sup>3</sup> m<sup>-3</sup> of soil</td></tr><tr><td>Volumetric wilting point for soil</td><td>m<sup>3</sup> m<sup>-3</sup> of soil</td></tr><tr><td>Soil Saturated hydraulic conductivity of soil</td><td>kg m<sup>-2</sup> s<sup>-1</sup></td></tr><tr><td>Saturated soil water pressure (used only for the Brooks and Corey, 1964, soil hydraulics)</td><td>m</td></tr><tr><td>Clapp-Hornberger exponent (used only for the Brooks and Corey, 1964, soil hydraulics)</td><td></td></tr><tr><td>1/α (used only for the van Genuchten, 1980, soil hydraulics)</td><td></td></tr><tr><td>1/(n - 1) (used only for the van Genuchten, 1980, soil hydraulics)</td><td></td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**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)</figcaption></figure>

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 + \Delta X$. The equations are then linearised by assuming that $\Delta X \ll 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:

$$E = \frac{\rho}{r_a + r_s} (Q_{\text{sat}}(T_*^i) - Q_1^i) + \frac{\rho}{r_a + r_s} \alpha' \Delta T_* - \frac{\rho}{r_a + r_s} \Delta Q_1 \quad (6)$$

where $\alpha'$ 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_{\text{soil}}$, inverse of resistance) is determined by the soil moisture concentration in the top soil layer:

$$g_{\text{soil}} = \frac{1}{100} \left( \frac{\theta_1}{\theta_c} \right)^2 \quad (7)$$

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) along with observational studies was used to develop this relationship whilst maintaining consistency with the critical point defined for vegetation (see Sect. 2.2, C. Taylor, personal communication, 2011).
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.

## 2.2 Photosynthesis and stomatal conductance

The leaf level stomatal conductance ($g_s$) and net photosynthetic uptake (A) are linked via the CO₂ diffusion equation:

$$A = g_s(C_c - C_i)/1.6 \qquad (8)$$

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₂ concentrations to leaf humidity deficit,

$$\frac{C_i - C_*}{C_c - C_*} = f_o \left( 1 - \frac{D}{D_*} \right) \qquad (9)$$

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₂ concentration, $C_i$.

$$A_P = \min(W_C, W_L, W_E) \qquad (10)$$

Leaf photosynthesis A, is related to the potential (non-stressed) leaf photosynthesis ($A_P$) as follows,

$$A = A_P \beta \qquad (11)$$

$\beta$ 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:

$$\beta = \begin{cases} 1 & \text{for } \theta \ge \theta_c \\ \frac{\theta - \theta_w}{\theta_c - \theta_w} & \text{for } \theta_w < \theta < \theta_c \\ 0 & \text{for } \theta \le \theta_w \end{cases} \qquad (12)$$

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, 2009). A more realistic scheme is provided by the multi-layer 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).

## 2.3 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.

## 2.4 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), for a staggered array of cubes; the canyon and the roof tiles both have the same roughness.
<figure><img src="image_6.png" /><figcaption>Fig. 3. Simulation of seasonal snow depth with JULES for multi-level 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.</figcaption></figure>

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.

## 2.5 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) Non-vegetated 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

$$z_0 = \omega h \quad (13)$$

There are two options to determine the surface albedo ($\alpha$) for vegetation. The simplest option is a bulk albedo:

$$\alpha = \alpha_b \exp(-kL) + \alpha_\infty [1 - \exp(-kL)] \quad (14)$$

where $\alpha_b$, the soil albedo, is a spatially varying ancillary field within JULES and $\alpha_\infty$ 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
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,

$$C_m = A_m + B_m L \quad (15)$$

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.

## 3 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 multi-layer scheme is expected to give better simulations of snow dynamics at sites with deep snow, with the possibility of mid-winter melt events and better simulations of soil temperatures at sites with low winter air temperatures.

### 3.1 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

$$r_s = \frac{\rho_i r^2}{0.03 D_f (1.79 + 3 U^{1/2}) I} \left( \frac{I}{I_{\max}} \right)^{0.4} \quad (16)$$

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

$$\Delta I = 0.7 (I_{\max} - I_0) (1 - e^{-S_f/I_{\max}}). \quad (17)$$

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).

### 3.2 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 ($\Delta z_1$), the thermal conductivity used in surface energy balance calculations is adjusted for insulation by snow according to

$$\lambda = \lambda_{\text{soil}} \left[ 1 + \frac{2d_s}{\Delta z_1} \left( \frac{\lambda_{\text{soil}}}{\lambda_{\text{snow}}} - 1 \right) \right]^{-1} \quad (18)$$

The heat flux between the surface layer and the second soil layer, of thickness $\Delta z_2$, is multiplied by a snow insulation factor

$$\zeta = \left( 1 + \frac{2d_s}{\Delta z_1 + \Delta z_2} \right)^{-1} \quad (19)$$
<table><caption>**Table 4.** Description of the various physics options within the JULES model as discussed in the identified sections.</caption><thead><tr><th>Physics</th><th>Section</th><th>Option</th></tr></thead><tbody><tr><td>Surface Exchange</td><td>2.1</td><td>No thermal inertia and conductive coupling<br>No thermal inertia and radiative coupling<br>Thermal inertia and radiative coupling<br>Thermal inertia, radiative coupling and snow under vegetation canopy</td></tr><tr><td>Canopy radiation and scaling</td><td>2.2</td><td>Big leaf<br>Multiple canopy layers<br>2 layer approximation to multiple canopy layers<br>Multiple canopy layers with variable leaf nitrogen</td></tr><tr><td>Albedo</td><td>2.5</td><td>Bulk albedos<br>Spectral albedos and snow ageing</td></tr><tr><td>Urban model</td><td>2.4</td><td>1 tile (bulk)<br>2 tiles (roofs and canyons)<br>MORUSES; Porson et al. (2010a,b)</td></tr><tr><td>Snow</td><td>3.2<br>3.3</td><td>zero layer model<br>multi-layer model</td></tr><tr><td>Soil hydraulics</td><td>4.4</td><td>Brooks and Corey (1964)<br>van Genuchten (1980)</td></tr><tr><td>Soil thermodynamics</td><td>4.5</td><td>Cox et al. (1999)<br>Dharssi et al. (2009)</td></tr><tr><td>Soil moisture super-saturation</td><td>4.6</td><td>Restricted drainage into layer<br>Infiltration into lower layer</td></tr><tr><td>Large scale hydrology</td><td>4.7.1<br>4.7.2</td><td>TOPMODEL Gedney and Cox (2003)<br>PDM (Moore, 1985)</td></tr></tbody></table>

For deeper snow, the surface conductivity is set equal to $\lambda_{\text{snow}}$ and the insulation factor is

$$ \zeta = (\Delta z_1 + \Delta z_2) \left[ (2d_s - \Delta z_1) \frac{\lambda_{\text{soil}}}{\lambda_{\text{snow}}} + 2 \Delta z_1 + \Delta z_2 \right]^{-1} \quad (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 multi-layer snow scheme.

### 3.3 Multi-layer snow model

The maximum number of layers ($N_{\text{max}}$) that are used for deep snow and their thickness $d_k$ ($k=1,...,N_{\text{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_{\text{max}} = 0$ forces the use of the zero-layer option for any depth of snow. For $d_s \ge d_1$, snow is represented by additional model layers on top of the soil if $N_{\text{max}} \ge 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 $\rho_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 $\rho_k d_k = I_k + W_k$. The increase in layer density due to compaction over a timestep of length $\delta t$ is calculated as

$$ \frac{\delta \rho_k}{\delta t} = \frac{\rho_k g M_k}{\eta_0} \exp \left( \frac{k_s}{T_m} - \frac{k_s}{T_k} - \frac{\rho_k}{\rho_0} \right) \quad (21) $$
<table><caption>**Table 5.** Default parameter values required by JULES for the standard vegetation surfaces.</caption><thead><tr><th>Parameter</th><th>Broadleaf trees</th><th>Needleleaf trees</th><th>C3 grasses</th><th>C4 grasses</th><th>Shrubs</th></tr></thead><tbody><tr><td>Snow-covered albedo for large LAI</td><td>0.15</td><td>0.15</td><td>0.60</td><td>0.60</td><td>0.40</td></tr><tr><td>Snow-covered albedo for zero LAI</td><td>0.30</td><td>0.30</td><td>0.80</td><td>0.80</td><td>0.80</td></tr><tr><td>Snow-free albedo for large LAI</td><td>0.10</td><td>0.10</td><td>0.20</td><td>0.20</td><td>0.20</td></tr><tr><td>Rate of change of vegetation roughness length with height</td><td>0.05</td><td>0.05</td><td>0.10</td><td>0.10</td><td>0.10</td></tr><tr><td>Minimum canopy capacity (kg m<sup>-2</sup>)</td><td>0.50</td><td>0.50</td><td>0.50</td><td>0.50</td><td>0.50</td></tr><tr><td>Rate of change of canopy capacity with LAI</td><td>0.05</td><td>0.05</td><td>0.05</td><td>0.05</td><td>0.05</td></tr><tr><td>Infiltration enhancement factor</td><td>4.00</td><td>4.00</td><td>2.00</td><td>2.00</td><td>2.00</td></tr><tr><td>Light extinction coefficient</td><td>0.50</td><td>0.50</td><td>0.50</td><td>0.50</td><td>0.50</td></tr><tr><td>Rootdepth (m)</td><td>3.00</td><td>1.00</td><td>0.50</td><td>0.50</td><td>0.50</td></tr></tbody></table>

where $k_s = 4000$ K, compactive viscosity $\eta_0 = 10^7$ Pa s, reference density $\rho_0 = 50$ kg m$^{-3}$, temperature $T_m = 273.15$ K, and $M_k = 0.5(I_k + W_k) + \sum_{i=1}^{k-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

$$C_k = I_k C_{\text{ice}} + W_k C_{\text{water}} \quad (22)$$

where $C_{\text{ice}} = 2100$ J K$^{-1}$ kg$^{-1}$ and $C_{\text{water}} = 4180$ J K$^{-1}$ kg$^{-1}$ are the specific heat capacities of ice and water, and the thermal conductivity is

$$\lambda_k = 2.22 \left( \frac{\rho_k}{\rho_{\text{water}}} \right)^{1.88} \quad (23)$$

where $\rho_{\text{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

$$H_k = \Gamma_k [T_k - T_{k+1} + \gamma (\delta T_k - \delta T_{k+1})] \quad (24)$$

where $\delta T_k$ is the increment in layer temperature over a timestep, $\gamma$ is the forward timestep weighting (0 for explicit and 1 for fully implicit timestepping), and

$$\Gamma_k = \left( \frac{d_k}{2\lambda_k} + \frac{d_{k+1}}{2\lambda_{k+1}} \right)^{-1} \quad (25)$$

is a layer thickness weighted thermal conductivity. For the lowest snow layer ($k = N$), $T_{N+1}$, $d_{N+1}$ and $\lambda_{N+1}$ are the temperature, thickness and conductivity of the surface soil layer. The increment in layer temperature over a timestep is

$$\delta T_k = \frac{\delta t}{C_k} (H_{k-1} - H_k) \quad (26)$$

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

$$\delta T_1 = \frac{\delta t}{C_1} [H_0 - \Gamma_1 (T_1 - T_{s1} + \gamma \delta T_1)] \quad (27)$$

with solution

$$\delta T_1 = \frac{[H_0 + \Gamma_1 (T_{s1} - T_1)] \delta t}{\gamma \Gamma_1 \delta t + C_1} \quad (28)$$

When there are $N > 1$ snow layers, increments in the layer temperatures are found as the solutions of the tridiagonal set of equations

$$b_1 \delta T_1 + c_1 \delta T_2 = [H_0 + \Gamma_1 (T_2 - T_1)] \delta t, \quad (29)$$

$$a_k \delta T_{k-1} + b_k \delta T_k + c_k \delta T_{k+1} = [\Gamma_{k-1}(T_{k-1} - T_k) + \Gamma_k(T_{k+1} - T_k)] \delta t \quad (30)$$

for $k = 2, \dots, N-1$, and

$$a_N \delta T_{N-1} + b_N \delta T_N = [\Gamma_{N-1}(T_{N-1} - T_N) + \Gamma_N(T_1 - T_N)] \delta t \quad (31)$$

with matrix elements

$$a_k = c_{k-1} = -\gamma \Gamma_{k-1} \delta t \quad (32)$$

and

$$b_k = C_k + \gamma (\Gamma_{k-1} + \Gamma_k) \delta t. \quad (33)$$

If the temperature of a layer is calculated to be above $T_m$, the layer ice mass is reduced by an amount

$$\delta I_k = \frac{C_k}{L_f} (T_k - T_m) \quad (34)$$

or the entire mass of the layer, whichever is least. The layer liquid mass is increased by the same amount and the layer
<figure><img src="image_7.png" /><figcaption>Fig. 4. Structure of the numerical discretisation over the layers for the temperatures and heat fluxes within the multi-level snow scheme in JULES.</figcaption></figure>

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 $\rho_{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

$$\delta W_k = \frac{C_k}{L_f}(T_m - T_k), \qquad (35)$$

increasing the ice content by the same amount, and increasing the temperature by

$$\delta T_k = \frac{L_f \delta W_k}{C_k}. \qquad (36)$$

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 $\rho_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

## 3.4 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

$$\alpha_s = \alpha_{cds} + k_1(\alpha_0 - \alpha_{cds})(T_* - T_c). \qquad (37)$$

For a tile with snow depth $d_s$, the albedo is a weighted average

$$\alpha = \alpha_0 + (\alpha_s - \alpha_0)(1 - e^{-d_s/d_m}) \qquad (38)$$

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

$$r(t + \delta t) = \left[ r(t)^2 + \frac{G_r}{\pi} \delta t \right]^{1/2} - [r(t) - r_0] \frac{S_f \delta t}{d_0} \qquad (39)$$

for snowfall rate $S_f$. The mass of fresh snow required to refresh the albedo is set to $2.5 \text{ kg m}^{-2}$. The empirical grain area growth rate, in $\mu\text{m}^2 \text{ s}^{-1}$, is

$$G_r = \begin{cases} 0.6 & T_* = T_m \text{ (melting snow)} \\ 0.06 & T_* < T_m, r < 150 \mu\text{m (cold fresh snow)} \\ A_s \exp(-4550/T_*) & T_* < T_m, r > 150 \mu\text{m (cold aged snow)} \end{cases} \qquad (40)$$

where $A_s = 0.23 \times 10^6 \mu\text{m}^2 \text{ s}^{-1}$. Snow albedos for diffuse visible and near-infrared radiation are calculated as

$$\alpha_{vis} = 0.98 - 0.002(r^{1/2} - r_0^{1/2}) \qquad (41)$$

and

$$\alpha_{nir} = 0.7 - 0.09 \ln \left( \frac{r}{r_0} \right). \qquad (42)$$

The zenith-angle dependence of albedos for direct-beam radiation with zenith cosine $\mu$ is represented by using an effective grain size,

$$r_e = [1 + 0.77(\mu - 0.65)]^2 r, \qquad (43)$$
in place of *r* in the equations for diffuse albedos (Eqs. 41 and 42).

For a tile with snow-free albedo $\alpha_0$, snowdepth $d_s$ and roughness length $z_0$, the albedo in each band is

$$ \alpha = f_s \alpha_s + (1 - f_s) \alpha_0 \quad (44) $$

where

$$ f_s = \frac{d_s}{d_s + 10z_0} \quad (45) $$

(Oleson et al., 2004). When driving data with separate direct-beam 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.

# 4 Hydrology and soil thermodynamics

JULES includes multi-layer, finite-difference models of sub-surface 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.

## 4.1 Surface hydrology

To account for the size of convective storms compared to gridsize, a rainfall rate is assumed to fall on a fraction $\epsilon_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:

$$ T_F = R \left(1 - \frac{C}{C_m}\right) \exp\left(-\frac{\epsilon_r C_m}{R \Delta t}\right) + R \frac{C}{C_m} \quad (46) $$

and the canopy water is updated by

$$ C^{(n+1)} = C^{(n)} + (R - T_F) \Delta t. \quad (47) $$

where $C_m$ is the maximum canopy water that can be held by the vegetation and $\Delta 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:

$$ Y = \begin{cases} R \frac{C}{C_m} \exp\left(-\frac{\epsilon_r K C_m}{R C}\right) + R \left(1 - \frac{C}{C_m}\right) \exp\left(-\frac{\epsilon_r C_m}{R \Delta t}\right) & K \Delta t \le C \\ R \exp\left[-\frac{\epsilon_r (K \Delta t + C_m - C)}{R \Delta t}\right] & K \Delta t > C \end{cases} \quad (48) $$

where the surface infiltration rate $K$ is equal to $\beta_s K_{hs}$ and $\beta_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:

$$ W_0 = \sum_j v_j (T_{Fj} + S_{mj} - Y_j). \quad (49) $$

## 4.2 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

$$ r_k = \frac{e^{-2z_{k-1}/d_r} - e^{-2z_k/d_r}}{1 - e^{-2z_r/d_r}}, \quad (50) $$

For transpiration $E'$, the flux extracted from soil layer $k$ is $e_k^0 E'$, where

$$ e_k^0 = \frac{r_k \beta_k}{\sum_k r_k \beta_k} \quad (51) $$

and

$$ \beta_k = \begin{cases} 1 & \theta_k \ge \theta_c \\ \frac{(\theta_k - \theta_w)}{(\theta_c - \theta_w)} & \theta_w < \theta_k < \theta_c \\ 0 & \theta_k \le \theta_w \end{cases} \quad (52) $$

is a soil moisture availability factor, defined similarly to Eq. (12), for a soil layer with unfrozen volumetric soil moisture concentration $\theta_k$.

## 4.3 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$:
$$C_a \Delta z_k \frac{dT_{sk}}{dt} = G_{k-1} - G_k - J_k \Delta z_k \quad (53)$$

where the fluxes are calculated as

$$G = \lambda_{\text{soil}} \frac{\partial T_s}{\partial z} \quad (54)$$

$$J = C_{\text{water}} W' \frac{\partial T_s}{\partial z} \quad (55)$$

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:

$$\frac{d\theta_k}{dt} = W'_{k-1} - W'_{k} - E'_{k} - R_{bk} \quad (56)$$

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

$$W' = K_h \left( \frac{\partial \Psi}{\partial z} + 1 \right) \quad (57)$$

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.

## 4.4 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):

$$\theta/\theta_s = (\Psi/\Psi_s)^{-1/b} \quad (58)$$

$$K_h = K_{hs} (\theta/\theta_s)^{2b+3} \quad (59)$$

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 and Corey (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:

$$\frac{\theta - \theta_r}{\theta_s - \theta_r} = \frac{1}{[1 + (\alpha_v \Psi)^n]^m} \quad (60)$$

$$K_h = K_{hs} S^{\xi} \left[ 1 - \left(1 - S^{1/m}\right)^m \right]^2 \quad (61)$$

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).

## 4.5 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:

$$\lambda = (\lambda_s - \lambda_{\text{dry}}) \theta/\theta_s + \lambda_{\text{dry}} \quad (62)$$

where

$$\lambda_s = \lambda_{\text{water}}^{\theta_u^s} \lambda_{\text{ice}}^{\theta_f^s} \lambda_{\text{dry}}^{\theta_s^s} / \lambda_{\text{air}}^{\theta_s^s} \quad (63)$$

where θ_u^s = θ_s[S_f/(S_u + S_f)], θ_f^s = θ_s - θ_f^s, 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
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:

$$ \lambda = (\lambda_s - \lambda_{\text{dry}}) K_e + \lambda_{\text{dry}} \quad (64) $$

where $K_e$ is the Kersten number

$$ K_e = \begin{cases} \log(\theta/\theta_s) + 1.0 & (\theta/\theta_s) \ge 0.1 \\ 0 & \text{otherwise} \end{cases} \quad (65) $$

$$ \lambda_s = \frac{\lambda_{\text{water}}^{\theta_s} \lambda_{\text{ice}}^{\theta_s}}{\lambda_{\text{water}}^{\theta_s}} \lambda_s^u \quad (66) $$

$$ \lambda_s^u = 1.58 + 12.4(\lambda_{\text{dry}} - 0.25) \quad (67) $$

with the constraint that $1.58 \le \lambda_s^u \le 2.2$. This equation for the thermal conductivity of unfrozen saturated soil ($\lambda_s^u$) 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).

## 4.6 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 super-saturated, 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 sub-grid 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 sub-surface 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.

## 4.7 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.

### 4.7.1 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 ($\bar{\lambda}_i$) flooding first. Topographic index relates to the upstream area draining into a locality and the local slope, which
<table><caption>**Table 6.** Default parameter values required by JULES for standard non-vegetation surfaces.</caption><thead><tr><th>Parameter</th><th>Urban</th><th>Water</th><th>Soil</th><th>Ice</th></tr></thead><tbody><tr><td>Snow-covered albedo</td><td>0.40</td><td>0.80</td><td>0.80</td><td>0.80</td></tr><tr><td>Snow-free albedo</td><td>0.18</td><td>0.06</td><td>-1.00*</td><td>0.75</td></tr><tr><td>Canopy capacity (kg m<sup>-2</sup>)</td><td>0.50</td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>Surface conductance (m s<sup>-1</sup>)</td><td>0.00</td><td>0.00</td><td>1 &times; 10<sup>-2</sup></td><td>1 &times; 10<sup>6</sup></td></tr><tr><td>Infiltration enhancement factor</td><td>0.10</td><td>0.00</td><td>0.50</td><td>0.00</td></tr><tr><td>Roughness length (m)</td><td>1.00</td><td>3 &times; 10<sup>-4</sup></td><td>3 &times; 10<sup>-4</sup></td><td>1 &times; 10<sup>-4</sup></td></tr><tr><td>Canopy heat capacity (J K<sup>-1</sup> m<sup>-2</sup>)</td><td>0.28 &times; 10<sup>6</sup></td><td>0.00</td><td>0.00</td><td>0.00</td></tr><tr><td>Fractional “canopy” coverage</td><td>1.00</td><td>0.00</td><td>0.00</td><td>0.00</td></tr></tbody></table>

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 $\bar{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

$$R_b = \bar{T}(\bar{z}_w) \exp(-\lambda_i) \quad (68)$$

where $T(\bar{z}_w)$ is the vertical transmissivity from the bottom of the column to the $\bar{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). $\bar{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 $\lambda_{ic} = \ln(R_{b\max}/R_b)$, where $R_{b\max}$ is the baseflow found from Eq. 68 with $\bar{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

$$f_{sat} = a_s \exp(-c_s f \lambda_{ic}) \quad (69)$$

where $a_s$ and $c_s$ are fitted parameters for each gridbox.

Saturation excess surface runoff ($R_{se}$) is then calculated as

$$R_{se} = f_{sat} W_0 \quad (70)$$

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 $\lambda_{ic} \le \lambda_i \le \lambda_{i\max}$ where $\lambda_{i\max}$ is a global parameter. At locations with larger values of $\lambda_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\lambda_{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).

### 4.7.2 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

<footnote>* The snow-free albedo for soil is initialised to -1 to allow it to be set through an ancillary field instead.</footnote>
$$f_{\text{sat}} = 1 - [1 - \theta/\theta_{\text{s}}]^{B/(B+1)} \quad (71)$$

where B is a shape parameter. Rse is then calculated using Eq. 70. In JULES, B is kept constant across the domain, as is the depth over which W and Wmax 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.

## 5 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., 2010, 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.
# Appendix A

## Definitions of symbols

<table><thead><tr><th>Symbol</th><th>Units</th><th>Equation</th><th>Definition</th></tr></thead><tbody><tr><td>A</td><td>mol CO<sub>2</sub> m<sup>-2</sup> s<sup>-1</sup></td><td>8, 11</td><td>net photosynthesis uptake</td></tr><tr><td>A<sub>m</sub></td><td>kg m<sup>-2</sup></td><td>15</td><td>puddling of water on soil surface and interception by leafless vegetation</td></tr><tr><td>A<sub>p</sub></td><td>mol CO<sub>2</sub> m<sup>-2</sup> s<sup>-1</sup></td><td>10, 11</td><td>potential leaf level photosynthesis</td></tr><tr><td>b</td><td></td><td>58, 59</td><td>Clapp and Hornberger (1978) soil exponent</td></tr><tr><td>B<sub>L</sub></td><td>kg C m<sup>-2</sup></td><td>5</td><td>leaf biomass</td></tr><tr><td>B<sub>m</sub></td><td>kg C m<sup>-2</sup></td><td>15</td><td>rate of change of water holding capacity with leaf area index</td></tr><tr><td>B<sub>w</sub></td><td>kg C m<sup>-2</sup></td><td>5</td><td>woody biomass</td></tr><tr><td>C</td><td>kg m<sup>-2</sup></td><td>46, 47, 48</td><td>canopy water</td></tr><tr><td>C<sub>a</sub></td><td>J m<sup>-3</sup> K<sup>-1</sup></td><td>53</td><td>volumetric heat capacity of the soil</td></tr><tr><td>C<sub>c</sub></td><td>Pa</td><td>8, 9</td><td>leaf surface carbon dioxide concentration</td></tr><tr><td>C<sub>i</sub></td><td>Pa</td><td>8, 9</td><td>internal leaf carbon dioxide concentration</td></tr><tr><td>C<sub>ice</sub></td><td>J kg<sup>-1</sup> K<sup>-1</sup></td><td>22</td><td>specific heat capacity of ice</td></tr><tr><td>C<sub>k</sub></td><td>J m<sup>-2</sup> K<sup>-1</sup></td><td>22, 26, 27, 28</td><td>areal heat capacity of the k-th snow layer</td></tr><tr><td>C<sub>L</sub></td><td>J kg<sup>-1</sup> K<sup>-1</sup></td><td>5</td><td>specific heat capacity of leaves</td></tr><tr><td>C<sub>m</sub></td><td>kg m<sup>-2</sup></td><td>15, 46, 48</td><td>vegetation canopy water holding capacity</td></tr><tr><td>c<sub>p</sub></td><td>J kg<sup>-1</sup> K<sup>-1</sup></td><td>2, 4</td><td>specific heat capacity of air</td></tr><tr><td>C<sub>s</sub></td><td>J m<sup>-2</sup> K<sup>-1</sup></td><td>1, 5</td><td>areal heat capacity associated with the surface material</td></tr><tr><td>C<sub>w</sub></td><td>J kg<sup>-1</sup> K<sup>-1</sup></td><td>5</td><td>specific heat capacity of wood</td></tr><tr><td>C<sub>water</sub></td><td>J kg<sup>-1</sup> K<sup>-1</sup></td><td>22, 55</td><td>specific heat capacity of water</td></tr><tr><td>C<sub>*</sub></td><td>Pa</td><td>9</td><td>CO<sub>2</sub> compensation point</td></tr><tr><td>D</td><td>kg kg<sup>-1</sup></td><td>9</td><td>leaf humidity deficit</td></tr><tr><td>D<sub>f</sub></td><td>m<sup>2</sup> s<sup>-1</sup></td><td>16</td><td>diffusivity of water vapour in air</td></tr><tr><td>d<sub>k</sub></td><td>m</td><td>25</td><td>depth of k-th snow layer</td></tr><tr><td>d<sub>m</sub></td><td>m</td><td>38</td><td>surface masking snow depth</td></tr><tr><td>d<sub>r</sub></td><td>m</td><td>50</td><td>root depth</td></tr><tr><td>d<sub>s</sub></td><td>m</td><td>18, 19, 20, 38, 45</td><td>snow depth</td></tr><tr><td>d<sub>0</sub></td><td>m</td><td>39</td><td>fresh snow depth</td></tr><tr><td>E</td><td>kg m<sup>-2</sup> s<sup>-1</sup></td><td>1, 3, 6</td><td>turbulent moisture flux</td></tr><tr><td>E'</td><td>kg m<sup>-2</sup> s<sup>-1</sup></td><td>56</td><td>evapotranspiration</td></tr><tr><td>f<sub>sat</sub></td><td></td><td>69, 70, 71</td><td>fraction of gridbox with saturated soil</td></tr><tr><td>G</td><td>W m<sup>-2</sup></td><td>1, 4, 53</td><td>soil heat flux</td></tr><tr><td>g</td><td>m s<sup>-2</sup></td><td>21</td><td>acceleration due to gravity</td></tr><tr><td>G<sub>r</sub></td><td>µm<sup>2</sup> s<sup>-1</sup></td><td>39, 40</td><td>snow grain area growth rate</td></tr><tr><td>g<sub>s</sub></td><td>m s<sup>-1</sup></td><td>8</td><td>leaf level stomatal conductance</td></tr><tr><td>g<sub>soil</sub></td><td>m s<sup>-1</sup></td><td>7</td><td>bare soil surface conductance</td></tr><tr><td>H</td><td>W m<sup>-2</sup></td><td>1, 2</td><td>turbulent heat flux</td></tr><tr><td>h</td><td>m</td><td>13</td><td>height of vegetation canopy</td></tr><tr><td>H<sub>k</sub></td><td>W m<sup>-2</sup></td><td>24, 26, 27, 28</td><td>conducted heat flux at the bottom of the k-th snow layer</td></tr><tr><td>I</td><td>kg m<sup>-2</sup></td><td>16, 17</td><td>intercepted snow load</td></tr><tr><td>I<sub>0</sub></td><td>kg m<sup>-2</sup></td><td>17</td><td>initial intercepted canopy snow load</td></tr><tr><td>I<sub>k</sub></td><td>kg m<sup>-2</sup></td><td>22, 34</td><td>ice content of k-th snow layer</td></tr><tr><td>I<sub>max</sub></td><td>kg m<sup>-2</sup></td><td>16, 17</td><td>snow interception canopy capacity</td></tr><tr><td>J</td><td>W m<sup>-3</sup></td><td>53, 55</td><td>vertical advective flux for soil moisture</td></tr><tr><td>K</td><td>kg m<sup>-2</sup> s<sup>-1</sup></td><td>48</td><td>surface infiltration rate</td></tr><tr><td>k</td><td></td><td>14</td><td>light extinction coefficient</td></tr><tr><td>K<sub>e</sub></td><td></td><td>64, 65</td><td>Kersten number</td></tr><tr><td>K<sub>h</sub></td><td>m s<sup>-1</sup></td><td>57, 59, 61</td><td>hydraulic conductivity</td></tr><tr><td>K<sub>hs</sub></td><td>m s<sup>-1</sup></td><td>59, 61</td><td>hydraulic conductivity for saturated soil</td></tr><tr><td>k<sub>l</sub></td><td>K<sup>-1</sup></td><td>37</td><td>snow ageing parameter</td></tr><tr><td>L</td><td>m<sup>2</sup> m<sup>-2</sup></td><td>14, 15</td><td>leaf area index</td></tr></tbody></table>
<table><thead><tr><th>Symbol</th><th>Units</th><th>Equation</th><th>Definition</th></tr></thead><tbody><tr><td>L<sub>c</sub></td><td>J kg<sup>-1</sup></td><td>1</td><td>latent heat of condensation of water at 0 °C</td></tr><tr><td>L<sub>f</sub></td><td>J kg<sup>-1</sup></td><td>35, 36</td><td>latent heat of fusion</td></tr><tr><td>Lw<sub>↓</sub></td><td>W m<sup>-2</sup></td><td>1</td><td>downward component of the longwave radiation</td></tr><tr><td>M<sub>k</sub></td><td>kg</td><td>21</td><td>mass of snow above the middle of the k-th snow layer</td></tr><tr><td>n</td><td></td><td>21</td><td>van Genuchten (1980) soil parameter</td></tr><tr><td>Q<sub>1</sub></td><td>kg kg<sup>-1</sup></td><td>3, 6</td><td>specific humidity at the reference atmospheric level</td></tr><tr><td>Q<sub>sat</sub>(T)</td><td>kg kg<sup>-1</sup></td><td>3, 6</td><td>saturated specific humidity at the temperature T</td></tr><tr><td>R</td><td>kg m<sup>-2</sup> s<sup>-1</sup></td><td>46, 47, 48</td><td>Precipitation rate</td></tr><tr><td>r</td><td>μm</td><td>16, 39, 41, 42</td><td>snow grain size</td></tr><tr><td>r<sub>a</sub></td><td>s m<sup>-1</sup></td><td>2, 3, 6</td><td>aerodynamic resistance</td></tr><tr><td>r<sub>a,can</sub></td><td>s m<sup>-1</sup></td><td>4</td><td>aerodynamic resistance between the surface canopy of vegetation and the underlying soil</td></tr><tr><td>R<sub>b</sub></td><td>kg m<sup>-2</sup> s<sup>-1</sup></td><td>56, 68</td><td>lateral runoff</td></tr><tr><td>r<sub>e</sub></td><td>μm</td><td>43</td><td>effective snow grain size</td></tr><tr><td>R<sub>es</sub></td><td>kg m<sup>-2</sup> s<sup>-1</sup></td><td>70</td><td>saturation excess surface runoff</td></tr><tr><td>r<sub>k</sub></td><td></td><td>50, 51</td><td>fraction of roots in the k-th soil layer</td></tr><tr><td>r<sub>s</sub></td><td>s m<sup>-1</sup></td><td>3, 6, 16</td><td>stomatal or surface moisture resistance</td></tr><tr><td>r<sub>0</sub></td><td>μm</td><td>39, 41, 42</td><td>fresh snow grain size</td></tr><tr><td>S<sub>f</sub></td><td>kg m<sup>-2</sup></td><td>17, 39</td><td>snowfall amount</td></tr><tr><td>S<sub>m</sub></td><td>kg m<sup>-2</sup> s<sup>-2</sup></td><td>49</td><td>snowmelt</td></tr><tr><td>Sw<sub>↓</sub></td><td>W m<sup>-2</sup></td><td>1</td><td>downward component of the solar radiation</td></tr><tr><td>t</td><td>s</td><td>1, 21, 26, 27, 28, 29, 30, 31, 32, 33, 39, 46, 47, 48, 53</td><td>time</td></tr><tr><td>T&#773;</td><td>m<sup>2</sup> s<sup>-1</sup></td><td>68</td><td>vertical transmissivity</td></tr><tr><td>T<sub>A</sub></td><td>K</td><td>2</td><td>reference level atmospheric temperature</td></tr><tr><td>T<sub>c</sub></td><td>K</td><td>37</td><td>snow albedo threshold temperature</td></tr><tr><td>T<sub>F</sub></td><td>kg m<sup>-2</sup> s<sup>-1</sup></td><td>46, 47, 49</td><td>throughfall</td></tr><tr><td>T<sub>k</sub></td><td>K</td><td>21, 24, 26, 27, 28</td><td>temperature of the k-th snow layer</td></tr><tr><td>T<sub>m</sub></td><td>K</td><td>21, 34, 35, 40</td><td>temperature of the melting point for water</td></tr><tr><td>T<sub>s</sub></td><td>K</td><td>54, 55</td><td>soil temperature</td></tr><tr><td>T<sub>s1</sub></td><td>K</td><td>4, 27, 28</td><td>temperature of the first soil level</td></tr><tr><td>T<sub>sk</sub></td><td>K</td><td>53</td><td>temperature of the k-th soil level</td></tr><tr><td>T<sub>*</sub></td><td>K</td><td>1, 2, 3, 4, 6, 37, 40</td><td>surface temperature</td></tr><tr><td>U</td><td>m s<sup>-1</sup></td><td>16</td><td>atmospheric wind speed</td></tr><tr><td>W<sub>C</sub></td><td>mol CO<sub>2</sub> m<sup>-2</sup> s<sup>-1</sup></td><td>10</td><td>rubisco limited rate for photosynthesis</td></tr><tr><td>W<sub>E</sub></td><td>mol CO<sub>2</sub> m<sup>-2</sup> s<sup>-1</sup></td><td>10</td><td>rate of transport of photosynthetic products (for C3 plants), or PEP-Carboxylase limitation for photosynthesis (for C4 plants)</td></tr><tr><td>W<sub>L</sub></td><td>mol CO<sub>2</sub> m<sup>-2</sup> s<sup>-1</sup></td><td>10</td><td>light limited rate for photosynthesis</td></tr><tr><td>W<sub>k</sub></td><td>kg m<sup>-2</sup></td><td>22, 35, 36</td><td>water content of k-th snow layer</td></tr><tr><td>W<sub>0</sub></td><td>kg m<sup>-2</sup> s<sup>-1</sup></td><td>49, 70</td><td>infiltration rate into the soil</td></tr><tr><td>W'</td><td>kg m<sup>-2</sup> s<sup>-1</sup></td><td>56, 57</td><td>vertical flux of soil water</td></tr><tr><td>Y</td><td>kg m<sup>-2</sup> s<sup>-1</sup></td><td>48, 49</td><td>surface runoff</td></tr><tr><td>z</td><td>m</td><td>54, 55, 57</td><td>soil depth</td></tr><tr><td>z<sub>k</sub></td><td>m</td><td>50, 53</td><td>depth of the k-th soil layer</td></tr><tr><td>z<sub>t</sub></td><td>m</td><td>50</td><td>total depth of soil column</td></tr><tr><td>z&#773;<sub>w</sub></td><td>m</td><td>68</td><td>mean water table depth</td></tr><tr><td>z<sub>0</sub></td><td>m</td><td>13, 45</td><td>roughness length for momentum</td></tr><tr><td>α</td><td></td><td>1, 14, 38, 44</td><td>surface albedo</td></tr><tr><td>α<sub>b</sub></td><td></td><td>14</td><td>bare soil albedo</td></tr><tr><td>α<sub>cds</sub></td><td></td><td>37</td><td>cold deep snow albedo</td></tr><tr><td>α<sub>nir</sub></td><td></td><td>42</td><td>diffuse near-infrared snow albedo</td></tr><tr><td>α<sub>s</sub></td><td></td><td>37, 38, 44</td><td>snow albedo</td></tr></tbody></table>
<table><thead><tr><th>Symbol</th><th>Units</th><th>Equation</th><th>Definition</th></tr></thead><tbody><tr><td>α<sub>v</sub></td><td>m<sup>-1</sup></td><td>60</td><td>van Genuchten (1980) soil parameter</td></tr><tr><td>α<sub>vis</sub></td><td></td><td>41</td><td>diffuse visible snow albedo</td></tr><tr><td>α<sub>0</sub></td><td></td><td>38.44</td><td>snow free albedo</td></tr><tr><td>α<sub>∞</sub></td><td></td><td>14</td><td>maximum canopy albedo for dense leaf coverage</td></tr><tr><td>α'</td><td>kg kg<sup>-1</sup> K<sup>-1</sup></td><td>6</td><td>δQ<sub>sat</sub>/δT</td></tr><tr><td>β</td><td></td><td>11, 12, 51, 52</td><td>soil moisture factor</td></tr><tr><td>γ</td><td></td><td>24, 27, 28, 32, 33</td><td>forward timestep weighting parameter</td></tr><tr><td>Δz<sub>i</sub></td><td>m</td><td>18, 19, 20</td><td>thickness of the i-th soil layer</td></tr><tr><td>ε</td><td></td><td>1, 4</td><td>surface emissivity</td></tr><tr><td>ε<sub>r</sub></td><td></td><td>46, 48</td><td>fraction of gridcell occupied by convective precipitation</td></tr><tr><td>ε<sub>s</sub></td><td></td><td>4</td><td>emissivity of the underlying soil surface</td></tr><tr><td>ζ</td><td></td><td>19, 20</td><td>snow insulation factor</td></tr><tr><td>ν</td><td></td><td>4, 49</td><td>fraction of vegetation</td></tr><tr><td>η<sub>0</sub></td><td>Pa s</td><td>21</td><td>compactive viscosity</td></tr><tr><td>θ</td><td>m<sup>3</sup> m<sup>-3</sup></td><td>12, 56, 58, 59, 60, 62, 65, 71</td><td>soil moisture concentration</td></tr><tr><td>θ<sub>1</sub></td><td>m<sup>3</sup> m<sup>-3</sup></td><td>7</td><td>soil moisture concentration in the top soil layer</td></tr><tr><td>θ<sub>c</sub></td><td>m<sup>3</sup> m<sup>-3</sup></td><td>7, 12</td><td>soil moisture concentration at critical point</td></tr><tr><td>θ<sub>r</sub></td><td>m<sup>3</sup> m<sup>-3</sup></td><td>60</td><td>residual soil moisture concentration</td></tr><tr><td>θ<sub>s</sub></td><td>m<sup>3</sup> m<sup>-3</sup></td><td>58, 59, 60, 62, 63, 65, 71</td><td>soil moisture concentration at saturation</td></tr><tr><td>θ<sub>w</sub></td><td>m<sup>3</sup> m<sup>-3</sup></td><td>12</td><td>soil moisture concentration at wilting point</td></tr><tr><td>λ</td><td>W m<sup>-1</sup> K<sup>-1</sup></td><td>18, 62, 64</td><td>thermal conductivity</td></tr><tr><td>λ<sub>air</sub></td><td>W m<sup>-1</sup> K<sup>-1</sup></td><td>63</td><td>thermal conductivity of air</td></tr><tr><td>λ<sub>dry</sub></td><td>W m<sup>-1</sup> K<sup>-1</sup></td><td>62, 63, 64, 67</td><td>thermal conductivity of dry soil</td></tr><tr><td>λ<sub>i</sub></td><td>ln(m)</td><td>68</td><td>topographic index</td></tr><tr><td>λ<sub>ic</sub></td><td>ln(m)</td><td>69</td><td>"critical" value of topographic index</td></tr><tr><td>λ<sub>ice</sub></td><td>W m<sup>-1</sup> K<sup>-1</sup></td><td>63, 66</td><td>thermal conductivity of ice</td></tr><tr><td>λ<sub>k</sub></td><td>W m<sup>-1</sup> K<sup>-1</sup></td><td>23, 25</td><td>thermal conductivity of k'th snow layer</td></tr><tr><td>λ<sub>snow</sub></td><td>W m<sup>-1</sup> K<sup>-1</sup></td><td>18, 20</td><td>thermal conductivity of the snow</td></tr><tr><td>λ<sub>s</sub></td><td>W m<sup>-1</sup> K<sup>-1</sup></td><td>62, 63, 64</td><td>thermal conductivity for saturated soil</td></tr><tr><td>λ<sub>soil</sub></td><td>W m<sup>-1</sup> K<sup>-1</sup></td><td>4, 18, 20</td><td>thermal conductivity of the soil</td></tr><tr><td>λ<sub>water</sub></td><td>W m<sup>-1</sup> K<sup>-1</sup></td><td>63, 66</td><td>thermal conductivity of water</td></tr><tr><td>μ</td><td></td><td>43</td><td>cosine of the zenith angle</td></tr><tr><td>ν<sub>j</sub></td><td></td><td>4, 49</td><td>fraction of gridbox covered by surface type j</td></tr><tr><td>ρ</td><td>kg m<sup>-3</sup></td><td>2, 3, 4, 6</td><td>density of air</td></tr><tr><td>ρ<sub>k</sub></td><td>kg m<sup>-3</sup></td><td>21, 23</td><td>density of the k-th snow layer</td></tr><tr><td>ρ<sub>i</sub></td><td>kg m<sup>-3</sup></td><td>16</td><td>density of ice</td></tr><tr><td>ρ<sub>water</sub></td><td>kg m<sup>-3</sup></td><td>23</td><td>density of water</td></tr><tr><td>ρ<sub>0</sub></td><td>kg m<sup>-3</sup></td><td>21</td><td>reference snow density</td></tr><tr><td>σ</td><td>W m<sup>-2</sup> K<sup>-4</sup></td><td>1, 4</td><td>Stefan Boltzmann constant</td></tr><tr><td>Ψ</td><td>m</td><td>57, 58, 60</td><td>soil water suction</td></tr><tr><td>Ψ<sub>s</sub></td><td>m</td><td>58</td><td>saturated soil water suction</td></tr><tr><td>ω</td><td>m m<sup>-1</sup></td><td>13</td><td>rate of change of roughness length with vegetation canopy height</td></tr></tbody></table>
Supplementary material related to this
article is available online at:

http://www.geosci-model-dev.net/4/677/2011/gmd-4-677-2011-supplement.pdf.

*Acknowledgements.* M. Best, M. Pryor, N. Gedney and O. Boucher were supported by the Joint DECC/Defra Met Office Hadley Centre Climate Programme (GA01101).

Edited by: J. C. Hargreaves

References

Bartlett, P. A., MacKay, M. D., and Verseght, D. L.: Modified snow algorithms in the Canadian Land Surface Scheme: Model runs and sensitivity analysis at three boreal forest stands., Atmos. Ocean, 44, 207–222, 2006.

Beljaars, A. C. M. and Holtslag, A. A. M.: Flux parametrisation over land surfaces for atmospheric models, J. Appl. Meteorol., 30, 327–341, 1991.

Best, M. J.: Representing urban areas within operational numerical weather prediction models, Bound.-Lay. Meteorol., 114, 91–109, 2005.

Best, M. J. and Hopwood, W. P.: Modelling the local surface exchange over a grass field site under stable conditions, Q. J. R. Meteorol. Soc., 127, 2033–2052, 2001.

Best, M. J., Beljaars, A., Polcher, J., and Viterbo, P.: A proposed structure for coupling tiled surfaces with the planetary boundary layer, J. Hydrometeorol., 5, 1271–1278, 2004.

Best, M. J., Cox, P. M., and Warrilow, D.: Determining the optimal soil temperature scheme for atmospheric modelling applications, Bound.-Lay. Meteorol., 114, 111–142, 2005.

Best, M. J., Grimmond, C. S. B., and Villani, M. G.: Evaluation of the urban tile in MOSES using surface energy balance observations, Bound.-Lay. Meteorol., 118, 503–525, 2006.

Betts, R. A., Boucher, O., Collins, M., Cox, P. M., Falloon, P., Gedney, N., Hemming, D., Huntingford, C., Jones, C. D., Sexton, D., and Webb, M.: Projected increase in continental runoff due to plant responses to increasing carbon dioxide, Nature, 448, 1037–1041, 2007.

Beven, K. J. and Kirkby, M. J.: A physically based, variable contributing area model of basin hydrology, Hydrol. Sci. Bull., 24, 43–69, 1979.

Blyth, E., Clark, D. B., Ellis, R., Huntingford, C., Los, S., Pryor, M., Best, M., and Sitch, S.: A comprehensive set of benchmark tests for a land surface model of simultaneous fluxes of water and carbon at both the global and seasonal scale, Geosci. Model Dev., 4, 255–269, doi:10.5194/gmd-4-255-2011, 2011.

Bohnenstengel, S. I., Evans, S., Clark, P., and Belcher, S. E.: Simulations of the London urban heat island, Q. J. Roy. Meteorol. Soc., 137, doi:10.1002/qj.855, 2011.

Boone, A. and Etchevers, P.: An Intercomparison of Three Snow Schemes of Varying Complexity Coupled to the Same Land Surface Model: Local-Scale Evaluation at an Alpine Site., J. Hydromet., 2, 374–394, 2001.

Bowling, L., Lettenmaier, D., Nijssen, B., Graham, L. P., Clark, D., El Maayar, M., Essery, R., Goers, S., Gusev, Y., Habets, F., van den Hurk, B., Jin, J., Kahan, D., Lohmann, D., Ma, X., Mahanama, S., Mocko, D., Nasonova, O., Niu, G. Y., Samuels-son, P., Shmakin, A., Takata, K., Verseghy, D., Viterbo, P., Xia, Y., Xue, Y., and Yang, Z. L.: Simulation of high-latitude hydrological processes in the Torne-Kalix basin: PILPS Phase 2(e) 1: Experiment description and summary intercomparisons, Glob. Planet. Change, 38, 1–30, 2003.

Brooks, R. H. and A. T. Corey: Hydraulic properties of porous media, Hydrological Papers 3, Colorado State Univ., Fort Collins, 1964.

Clapp, R. B. and Hornberger, G. M.: Empirical equations for some soil hydraulic properties, Water Resour. Res., 14, 601–604, 1978.

Clark, D. B. and Gedney, N.: Representing the effects of subgrid variability of soil moisture on runoff generation in a land surface model, J. Geophys. Res., 113, D10111, doi:10.1029/2007JD008940, 2008.

Clark, D. B., Mercado, L. M., Sitch, S., Jones, C. D., Gedney, N., Best, M. J., Pryor, M., Rooney, G. G., Essery, R. L. H., Blyth, E., Boucher, O., Harding, R. J., and Cox, P. M.: The Joint UK Land Environment Simulator (JULES), model description – Part 2: Carbon fluxes and vegetation dynamics, Geosci. Model Dev., 4, 701–722, doi:10.5194/gmd-4-701-2011, 2011.

Collatz, G. J., Ball, J. T., Grivet, C., and Berry, J. A.: Physiological and environmental regulation of stomatal conductance, photosynthesis and transpiration: a model that includes a laminar boundary layer, Agric. Forest Meteorol., 54, 107–136, 1991.

Collatz, G. J., Ribas-Carbo, M., and Berry, J. A.: Coupled photosynthesis-stomatal conductance model for leaves of C₄ plants, Aust. J. Plant Physiol., 19, 519–538, 1992.

Cosby, B. J., Hornberger, G. M., Clapp, R. B., and Ginn, T. R.: A statistical exploration of the relationships of soil moisture characteristics to the physical properties of soils, Water Resour. Res., 20, 682–690, 1984.

Cox, P. M., Huntingford, C., and Harding, R. J.: A canopy conductance and photosynthesis model for use in a GCM land surface scheme, J. Hydrol., 213, 79–94, 1998.

Cox, P. M., Betts, R. A., Bunton, C. B., Essery, R. L. H., Rowntree, P. R., and Smith, J.: The impact of new land surface physics on the GCM simulation of climate and climate sensitivity, Clim. Dynam., 15, 183–203, 1999.

Cox, P. M., Betts, R. A., Jones, C. D., Spall, S. A., and Totterdell, I. J.: Acceleration of global warming due to carbon-cycle feedbacks in a coupled climate model, Nature, 408, 184–187, 2000.

Cox, P. M., Harris, P., Huntingford, C., Betts, R. A., Collins, M., Jones, C. D., Jupp, T. E., Marengo, J., and Nobre, C.: Increasing risk of Amazonian drought due to decreasing aerosol pollution, Nature, 453, 212–216, 2008.

Cullen, M. J. P.: The unified forecast/climate model, Meteorol. Mag., 122, 81–94, 1993.

Deardorff, J. W.: Efficient prediction of ground surface temperature and moisture with inclusion of a layer of vegetation, J. Geophys. Res., 83, 1889–1903, 1978.

Dharssi, I., Vidale, P. L., Verhoef, A., Macpherson, B., Jones, C., and Best, M.: New soil physical properties implemented in the Unified Model at PS18, Met Office Technical Report 528, Met Office, Exeter, UK, 2009.

Dolman, A. J. and Gregory, D.: The parametrization of rainfall
interception in GCMs, Q. J. R. Meteorol. Soc., 118, 455–467, 1992.

Douville, H., Royer, J.-Y., and Mahfouf, J.-F.: A new snow param- eterization for the Meteo-France climate model., Clim. Dynam., 12, 21–35, 1995.

Dyer, A. J.: A review of flux-profile relationships, Bound.-Lay. Me- teorol., 7, 363–372, 1974.

Ek, M. B., Mitchell, K. E., Lin, Y., Rogers, E., Grunmann, P., Ko- ren, V., Gayno, G., and J. D. Tarpley, J. D.: Implementation of Noah land surface model advances in the National Centers for Environmental Prediction operational mesoscale Eta model, J. Geophys. Res., 108, 163–166, doi:10.1029/2002JD003296, 2003.

Essery, R. and Clark, D.: Developments in the MOSES 2 land- surface model for PILPS 2e, Glob. Planet. Change, 38, 161–164, 2003.

Essery, R. L. H., Best, M. J., Betts, R. A., Cox, P. M., and Taylor, C. M.: Explicit representation of subgrid heterogeneity in a GCM land surface scheme, J. Hydromet., 4, 530–543, 2003a.

Essery, R. L. H., Pomeroy, J. W., Parvianen, J., and Storck, P.: Sub- limation of snow from boreal forests in a climate model, J. Cli- mate, 16, 1855–1864, 2003b.

Essery, R., Rutter, N., Pomeroy, J., Baxter, R., Stahli, M., Gustaf- son, D., Barr, A., Bartlett, P., and Elder, K.: SNPWMIP2, An evalutation of froest snow process simulations, Bull. Amer. Met. Soc., 90, 1120–1135, 2009.

Evans, S., Hudson-Smith, A., and Batty, M.: Virtual London and beyond, Cybergeo : European Journal of Geography, A Selec- tion of the best Articles (SAGEO 2005), 359, online on 27 Octo- ber 2006, modified on 04 July 2007, consulted on 16 June 2011, available at: http://cybergeo.revues.org/2871, last access: 2 Au- gust 2011, 2005.

Fisher, J. B., Sitch, S., Malhi, Y., Fisher, R. A., Huntingford, C., and Tan, S. Y.: Carbon cost of plant nitrogen acquisition: A mechanistic, globally applicable model of plant nitrogen up- take, retranslocation, and fixation, Global Biogeochem. Cy., 24, GB1014, doi:10.1029/2009GB003621, 2010.

Gedney, N. and Cox, P. M.: The sensitivity of global climate model simulations to the representation of soil moisture heterogeneity, J. Hydromet., 4, 1265–1275, 2003.

Gedney, N., Cox, P. M., and Huntingford, C.: Climate feed- back from wetland methane emissions, Geophys. Res. Lett., 31, L20503, doi:10.1029/2004GL020919, 2004.

Gedney, N., Cox, P. M., Betts, R., Boucher, O., Huntingford, C., and Stott, P. A.: Detection of a direct carbon dioxide effect in continental river runoff records, Nature, 439, 835–888, 2006.

Grimmond, C. S. B., Blackett, M., Best, M., Barlow, J., Baik, J.-J., Belcher, S. E., Bohnenstengel, S. I., Calmet, I., Chen, F., Dan- dou, A., Fortuniak, K., Gouvea, M. L., Hamdi, R., Hendry, M., Kawai, T., Kawamoto, Y., Kondo, H., Krayenhoff, E. S., Lee, S.-H., Loridan, T., Martilli, A., Masson, V., Miao, S., Oleson, K., Pigeon, G., Porson, A., Ryu, Y.-H., Salamanca, F., Shashua- Bar, L., Steeneveld, G.-J., Trombou, M., Voogt, J., Young, D., and Zhang, N.: The international urban energy balance models comparison project: first results from phase 1, J. Appl. Me- torol. Climatol., 49, 1268–1292, doi:10.1175/2010JAMC2354.1, 2010.

Grimmond, C. S. B., Blackett, M., Best, M., Baik, J.-J., Belcher, S. E., Beringer, J., Bohnenstengel, S. I., Calmet, I., Calmet, I.,

Chen, F., Coutts, A., Dandou, A., Fortuniak, K., Gouvea, M. L., Hamdi, R., Hendry, M., Kanda, M., Kawai, T., Kawamoto, Y., Kondo, H., Krayenhoff, E. S., Lee, S.-H., Loridan, T., Martilli, A., Masson, V., Miao, S., Oleson, K., Ooka, K., Pigeon, G., Por- son, A., Ryu, Y.-H., Salamanca, F., Steeneveld, G.-J., Trombou, M., Voogt, J., Young, D. T., and Zhang, N.: Initial results from phase 2 of the international urban energy balance model compar- ison, Int. J. Climatol., 30, 244–272, doi:10.1002/joc.2227, 2011.

Harman, I. N. and Belcher, S. E.: The surface energy balance and boundary layer over urban street canyons, Q. J. R. Meteorol. Soc., 132, 2749–2768, 2006.

Harman, I. N., Barlow, J., and Belcher, S. E.: Scalar Fluxes from Ur- ban Street Canyons Part II: Model, Bound.-Lay. Meteorol., 113, 387–410, doi:10.1023/B:BOUN.0000045526.07270.a3, 2004.

Jackson, T. L., Feddema, J. J., Oleson, K. W., Bonan, G. B., and Bauer, J. T.: Parameterization for global of urban characteris- tics climate modeling, Ann. Assoc. Am. Geogr., 100, 848–865, doi:10.1080/00045608.2010.497328, 2010.

Jacobs, C. M. J.: Direct impact of atmopsheric CO₂ enrichment on regional transpiration, Ph.D. thesis, Wageningen Agricultural University, 1994.

Jogireddy, V., Cox, P. M., Huntingford, C., Harding, R. J., and Mer- cado, L. M.: An improved description of canopy light intercep- tion for use in a GCM land-surface scheme: calibration and test- ing against carbon fluxes at a coniferous forest, Hadley Centre Technical Note 63, Hadley Centre, Met Office, Exeter, UK, 2006.

Johansen, O.: Thermal conductivity of soils, Ph.d. thesis, University of Trondheim, Norway, 1975.

Kojima, K.: Densification of seasonal snow cover, Proc. Interna- tional Conference on Low Temperature Science, August 14–19, 1966, Sapporo, Japan, 929–952, 1967.

Koster, R. D., Suarez, M. J., Ducharne, A., Stieglitz, M., and Ku- mar, P.: A catchment-based approach to modeling land surface processes in a general circulation model, J. Geophys. Res., 105, 24809–24822, 2000.

Leuning, R.: A critical appraisal of a combined stomatal- photosynthesis model for C₃ plants, Plant, Cell Environ., 18, 357–364, 1995.

Lynch-Stieglitz, M.: The development and validation of a simple snow model for the GISS GCM, J. Climate, 7, 1842–1855, 1994.

Macdonald, R. W., Griffiths, R. F., and Hall, D.: An improved method for the estimation of surface roughness of obstacle ar- rays, Atmos. Environ., 32, 1857–1864, 1998.

Mahfouf, J. F. and Noilhan, J.: Comparative-study of various for- mulations of evaporation from bare soil using insitu data, J. Appl. Meterol., 30, 1354–1365, 1991.

Masson, V.: A physically-based scheme for the urban energy budget in atmospheric models, Bound.-Lay. Meteorol., 41, 1011–1026, 2000.

Mercado, L. M., Huntingford, C., Gash, J. H. C., Cox, P. M., and Jogireddy, V.: Improving the representation of radiative intercep- tion and photosynthesis for climate model applications, Tellus B, 59, 553–565, 2007.

Mercado, L. M., Bellouin, N., Sitch, S., Boucher, O., Huntingford, C., Wild, M., and Cox, P. M.: Impact of changes in diffuse ra- diation on the global land carbon sink, Nature, 458, 1014–1018, doi:10.1038/nature07949, 2009.

Monin, A. S. and Obukhov, A. M.: Basic regularity in turbu- lent mixing in the surface layer of the atmosphere, Moscow,
Ak. Nauk, Geof. Inst., 24, 163–187, 1954.

Monsi, M. and Saeki, T.: Ueber den Lichtfaktor in den Pflanzengesellschaften und seine Bedeutung fuer die Stoffproduktion, Jpn. J. Bot., 14, 22–52, 1953.

Moore, R. J.: The probability-distributed principle and runoff production at point and basin scales, Hydrol. Sci. J., 30, 273–297, 1985.

Noilhan, J. and Planton, S.: A simple parameterization of land surface processes for meteorological models, Mon. Weather Rev., 117, 536–549, 1989.

Oleson, K. W., Dai, Y., Bonan, G., Bosilovich, M., Dickinson, R., Dirmeyer, P., Hoffman, F., Houser, P., Levis, S., Niu, G.-Y., Thornton, P., Vertenstein, M., Yang, Z.-L., and Zeng, X.: Technical description of the Community Land Model (CLM), 173 pp., 2004.

Oleson, K. W., Lawrence, D. M., Bonan, G. B., Flanner, M. G., Kluzek, E., Lawrence, P. J., Levis, S., Swenson, S. C., Thornton, P. E., Dai, A., Decker, M., Dickinson, R., Feddema, J., Heald, C. L., Hoffman, F., Lamarque, J. F., Mahowald, N., Niu, G. Y., Qian, T., Randerson, J., Running, S., Sakaguchi, K., Slater, A., Stockli, R., Wang, A., Yang, Z. L., Zeng, X., and Zeng, X.: Technical Description of version 4.0 of the Community Land Model (CLM), 257 pp., 2010.

Pitman, A. J., Liang, Z. L., Cogley, J. G., and Henderson-Sellers, A.: Description of bare essentials of surface transfer for the Bureau of Meteorological Research Centre AGCM, 1991.

Porson, A., Harman, I. N., Bohnenstengel, S. I., and Belcher, S. E.: How many facets are needed to represent the surface energy balance of an urban area?, Bound.-Lay. Meteorol., 132, 107–128, 2009.

Porson, A., Clark, P. A., Harman, I. N., Best, M. J., and Belcher, S. E.: Implementation of a new urban energy budget scheme in the MetUM. Part I. Description and idealized simulations, Q. J. R. Meteorol. Soc., 136, 1514–1529, 2010a.

Porson, A., Clark, P. A., Harman, I. N., Best, M. J., and Belcher, S. E.: Implementation of a new urban energy budget scheme in the MetUM, Part II. Validation against observations and model intercomparison, Q. J. R. Meteorol. Soc., 136, 1530–1542, 2010b.

Schlosser, C. A., Slater, A. G., Pitman, A. J., Vinnikov, K. Y., Henderson-Sellers, A., Speranskaya, N. A., and Mitchell, K.: Simulations of a boreal grassland hydrology at Valdai, Russia: PILPS phase 2(d), Mon. Weather Rev., 128, 301–321, 2000.

Sellers, P. J.: Canopy Reflectance, Photosynthesis, and Transpiration III, A Reanalysis Using Improved Leaf Models and a New Canopy Integration Scheme, Int. J. Remote Sens, 6, 1335–1372, 1985.

Sitch, S., Smith, B., Prentice, I. C., Arneth, A., Bondeau, A., Cramer, W., Kaplan, J. O., Levis, S., Lucht, W., Sykes, M. T., Thonicke, K., and Venevsky, S.: Evaluation of ecosystem dynamics, plant geography and terrestrial carbon cycling in the LPJ dynamic global vegetation model, Glob. Change Biol., 9, 161–185, 2003.

Sitch, S., Cox, P. M., Collins, W. J., and Huntingford, C.: Indirect radiative forcing of climate change through ozone effects on the land-carbon sink, Nature, 448, 791–794, 2007.

Sivapalan, M., Wood, E. F., and Beven, K. J.: On hydrologic similarity. 2. A scaled model of storm runoff production. 3. A dimensionless flood frequency model using a generalised geomorphologic unit hydrograph and partial area runoff generation, Water Resour. Res., 26, 43–58, 1987.

Smith, B., Prentice, I. C., and Sykes, M. T.: Representation of vegetation dynamics in the modelling of terrestrial ecosystems: comparing two contrasting approaches within European climate space., Global Ecol. Biogeogr., 10, 621–637, doi:10.1046/j.1466-822X.2001.t01-1-00256.x, 2001.

Sokolov, A. P., Kicklighter, D. W., Melillo, J. M., Felzer, B. S., Schlosser, C. A., and Cronin, T. W.: Consequences of considering carbon-nitrogen interactions on the feedbacks between climate and the terrestrial carbon cycle, J. Climate, 21, 3776–3796, doi:10.1175/2008JCLI2038.1, 2008.

Storck, P., Lettenmaier, D. P., and Bolton, S. M.: Measurement of snow interception and canopy effects on snow accumulation and melt in a mountainous maritime climate, Oregon, United States, Water Resour. Res., 38, 1223, doi:10.1029/2002WR001281, 2002.

Taylor, C. M. and Clark, D. B.: The diurnal cycle and African eastery waves: A land surface perspective, Q. J. R. Meteorol. Soc., 127, 845–867, 2001.

Thornton, P. E., Lamarque, J. F., Rosenbloom, N. A., and Mahowald, N. M.: Influence of carbon-nitrogen cycle coupling on land model response to CO2 fertilization and climate variability, Global Biogeochem. Cy., 21, GB4018, doi:10.1029/2006GB002868, 2007.

Thornton, P. E., Doney, S. C., Lindsay, K., Moore, J. K., Mahowald, N., Randerson, J. T., Fung, I., Lamarque, J. F., Feddema, J. J., and Lee, Y. H.: Carbon-nitrogen interactions regulate climate-carbon cycle feedbacks: results from an atmosphere-ocean general circulation model, Biogeosciences, 6, 2099–2120, doi:10.5194/bg-6-2099-2009, 2009.

van Genuchten, M. T.: A closed-form equation for predicting the hydraulic conductivity of unsaturated soils., Soil Sci. Soc. America J., 44, 892–898, 1980.

Verseghy, D. L.: CLASS – a Canadian land surface scheme for GSMS, I. Soil model, Int. J. Climatol., 11, 347–370, 1991.

Verseghy, D. L., Lazare, M., and McFarlane, N.: CLASS – A Canadian land surface scheme for GSMS, II. vegetation model and coupled runs, Int. J. Climatol., 13, 111–133, 1993.

Viterbo, P. and Beljaars, A. C. M.: An Improved Land-Surface Parameterization Scheme in the Ecmwf Model and Its Validation, J. Climate, 8, 2716–2748, 1995.

Wiscombe, W. J. and Warren, S. G.: A model for the spectral albedo of snow, I. Pure snow, J. Atmos. Science, 37, 2712–2733, 1980.

Yen, Y.-C.: Review of thermal properties of snow, ice and sea ice, 27 pp., 1981.

Zaehle, S., Friend, A. D., Friedlingstein, P., Dentener, F., Peylin, P., and Schulz, M.: Carbon and nitrogen cycle dynamics in the O-CN land surface model: 2. Role of the nitrogen cycle in the historical terrestrial carbon balance, Global Biogeochem. Cy., 24, GB1006, doi:10.1029/2009GB003522, 2010.