Edinburgh WRFv3.2-SPAv2: development and validation of a coupled ecosystem-atmosphere model, scaling from surface fluxes of CO2 and energy to atmospheric profiles

. The Weather Research and Forecasting meteorological (WRF) model has been coupled to the Soil–Plant–Atmosphere (SPA) terrestrial ecosystem model, to produce WRF-SPA. SPA generates realistic land– atmosphere exchanges through fully coupled hydrological, carbon and energy cycles. The addition of a land surface model (SPA) capable of modelling biospheric CO 2 exchange allows WRF-SPA to be used for investigating the feedbacks between biosphere carbon balance, meteorology, and land use and land cover change. We have extensively validated WRF-SPA using multi-annual observations of air temperature, turbulent ﬂuxes, net radiation and net ecosystem exchange of CO 2 at three sites, representing the dominant vegetation types in Scotland (forest, managed grassland and arable agriculture). For example air temperature is well simulated across all sites (forest R 2 = 0.92, RMSE = 1.7 ◦ C, bias = 0.88 ◦ C; managed grassland R 2 = 0.73, RMSE = 2.7 ◦ C, bias = − 0.30 ◦ C; arable agriculture R 2 = 0.82, RMSE = 2.2 ◦ C, bias = 0.46 ◦ C; RMSE, root mean square error). WRF-SPA generates more realistic seasonal behaviour at the site level compared to an unmodiﬁed version of WRF, such as improved simulation of seasonal transitions in latent heat ﬂux in arable systems. WRF-SPA also generates realistic seasonal CO 2 exchanges across all sites. WRF-SPA is also able to realistically model atmospheric proﬁles of CO 2 over Scotland, spanning a 3 yr (2004–2006), capturing both proﬁle structure, indicating realistic transport,


Introduction
The land surface is a key driver of climate and biogeochemical cycles at regional and global scales (Pielke et al., 1997;Cox et al., 2000;Esau and Lyons, 2002). Understanding land-atmosphere interactions for matter and energy is essential, as changes in global climate have likely led to significant but poorly understood changes in the global carbon balance (Forster et al., 2007). The direction and magnitude of future changes to the Earth system will depend on feedbacks between biogeochemistry and climate, linked to human modification of the land surface related to land use and land cover change. However, the importance of the land surface has been relatively poorly explored in terms of the climatebiogeochemical coupling .
Simulation models provide the only effective means to investigate the dynamics of land-atmosphere interactions. However, models typically have a genesis in either the atmospheric or in the terrestrial biogeochemistry/ecological communities. To overcome deficiencies arising from these alternate perspectives has required coupling the most advanced descriptions of atmospheric dynamics with ecosystem processes. Coupled land-atmosphere models are now used to investigate global-and regional-scale exchange between the land and atmosphere, and have highlighted the complex feedbacks that result (e.g. Ciais et al., 2005;Friedlingstein et al., 2006;Betts et al., 2007;Friedlingstein and Prentice, 2010).
Mesoscale models in particular are useful for studying regional-scale processes due to their ability to operate at high spatial resolutions, allowing accurate prediction of mesoscale circulation phenomena that have a significant impact on regional transport e.g. coastal winds, katabatic winds (Nicholls et al., 2004;Riley et al., 2005;Ahmadov et al., 2009), and complex forcing of heterogeneous vegetation . Mesoscale models also provide a means to upscale land surface exchanges to observations within the planetary boundary layer (PBL), of regionally integrated exchange (e.g. tall towers or aircraft observations) (Ahmadov et al., 2009;Tolk et al., 2009). Atmospheric tracers of CO 2 for respiratory and photosynthetic exchange can be transported through the atmosphere making it possible to gain information on how variations in land use and ecosystem coverage contribute to observations of regional exchange of CO 2 (Tolk et al., 2009).
Over time land surface models (LSMs), used in coupled land-atmosphere models, have been upgraded to include more processes and improved parameterisation. The current generation of LSMs commonly represent both biogeophysical and biogeochemical processes (Bonan, 2008). For example the widely used mesoscale model Weather Research and Forecasting (WRF, Skamarock et al., 2008) uses the advanced Noah-MP LSM (Niu et al., 2011). Noah-MP includes detailed parameterisation for radiative transfer of sunlit and shaded leaf area for a big leaf canopy, a semi-empirical model of stomatal conductance and a carbon model allowing dynamic ecosystem phenology. However there remains uncertainty in the predicted response of the terrestrial ecosystem, net ecosystem exchange (NEE) of CO 2 and feedbacks mediated through the close coupling between the land surface and PBL processes (Friedlingstein et al., 2006;Bonan, 2008;Sitch et al., 2008;Friedlingstein and Prentice, 2010).
Feedbacks between the land and atmosphere are highly complex and non-linear, requiring an increasingly mechanistic approach to provide realistic responses of key ecosystem processes (i.e. photosynthesis, respiration and evapotranspiration) (Tuzet et al., 2003;Bonan, 2008;Sprintsin et al., 2012). Significant improvements in the prediction of photosynthesis and stomatal conductance are made by using a multi-layer canopy in radiative transfer, where direct and diffuse radiation, and sunlit and shaded leaf areas are modelled (Wang and Leuning, 1998;Dai et al., 2003;Sprintsin et al., 2012). Moreover, stomatal conductance is a key determinant of land surface hydrological and carbon cycles, and energy balance (Avissar, 1998). Semi-empirical representations which link photosynthesis to atmospheric demand for water, such as Ball-Berry (Collatz et al., 1991) or its variants is widely used in LSMs (e.g. Sellers et al., 1996;Oleson et al., 2010;Best et al., 2011;Niu et al., 2011). Models such as Ball-Berry do not couple atmospheric demand for water to available water supply from the soil, thus lacking important ecophysiological processes (Tuzet et al., 2003).
Furthermore, land-atmosphere feedbacks are often distinct to each land cover type, particularly at longer timescales (Stoy et al., 2009). Therefore heterogeneity in the land surface needs to be realistically modelled to accurately predict regional-scale exchange (Avissar, 1998;Schomburg et al., 2012). These feedbacks are mediated through biogeochemical and biogeophysical processes (Bonan, 2008) that are related to plant phenology (e.g. canopy height and leaf area index). Realistic modelling of phenology is particularly important as most land surface systems have a distinct annual cycle. Human intervention adds further complexity to annual cycles, for example in agricultural systems. Crop modelling in LSMs has received little attention until recently (Sus et al., 2010) and crops are often modelled as a natural grassland (Osborne et al., 2007). Croplands play a significant role in global biogeochemical and biogeophysical cycles (Bondeau et al., 2007;Denman et al., 2007), significantly impacting both turbulent fluxes (Van den Hoof et al., 2011) and surface albedo . Due to the importance of cropland a number of LSMs have been modified to include developmental crop models to improve the presentation of crop phenology and management (Lokupitiya et al., 2009;Sus et al., 2010;Van den Hoof et al., 2011;Levis et al., 2012).
The Soil-Plant-Atmosphere (SPA, Williams et al., 1996) model is a mechanistic ecosystem model. SPA uses a vertically distributed canopy model allowing variation of photosynthetic parameters through the canopy, based on field measurements, and a multi-layer radiative transfer scheme that models the distribution of direct and diffuse radiation, and sunlit and shaded leaf areas (Williams et al., 1998). SPA also uses a mechanistic model of stomatal conductance linking atmospheric demand and water availability from the soil through the plant, explicitly coupling plant carbon and hydrological cycles (Williams et al., 1996(Williams et al., , 2001. Unlike a Ball-Berry model, SPA is parameterised directly from ecophysiological measurements such as rooting depth, plant hydraulic conductance and canopy structure (e.g. Wright et al., 2012). A developmental crop model has also been included in SPA allowing for inclusion of this highly important ecosystem under human management (Sus et al., 2010).
In this paper we describe a novel coupling between the WRF and SPA, forming WRF-SPA. An overview of both models, their strengths, any modifications made to either code or their forcing data is described. WRF is a state-of-theart non-hydrostatic mesoscale meteorological model (Skamarock et al., 2008), it is considered to be one of the best models of its type available Steeneveld et al., 2011). WRF-SPA uses a number of CO 2 tracer pools to upscale land surface processes, photosynthesis and respiration, specific to each land surface type.
Here WRF-SPA is validated at both local and regional scales, using observations from Scotland over multiple years, to evaluate the model against a range of meteorological conditions. At the local-scale surface observations of meteorological conditions and fluxes of CO 2 , heat and water have been used for validation. Regional-scale validation used aircraft profile measurements of atmospheric CO 2 concentrations. The unmodified WRFv3.2 and WRF-SPA have been compared, at the local scale, to investigate the impacts of the addition of a new LSM. We aim to answer a number of specific questions: i. Can WRF-SPA realistically model surface meteorological variables and fluxes across a multi-annual period?
ii. Does WRF-SPA scale realistically from surface measurements to regional-scale observations, specifically aircraft profiles?
iii. Does WRF-SPA lead to an improvement in surface fluxes compared to the unmodified WRFv3.2?

Model description: WRF
The Weather Research and Forecasting (v3.2) (http://www.mmm.ucar.edu/wrf/users/, accessed 19 October 2009) model is a well supported and rapidly developing high resolution non-hydrostatic meteorological model (Skamarock et al., 2008). WRF is designed to be highly adaptable, with a portable code for use on massively parallel systems and a modular structure to allow for tailoring to specific uses. The model has been extensively validated over a range of locations around the world (Ahamdov et al., 2007;Ahmadov et al., 2009;Borge et al., 2008;Zhang, 2008;Wang et al., 2009) and performs favourably in comparison to other commonly used regional meteorological models Steeneveld et al., 2011). Here we use the Advanced Research WRF (ARW) dynamical solver which uses non-hydrostatic equations, allowing horizontal resolutions of < 1 km.

Atmospheric CO 2 tracers
WRF-SPA has been modified with the addition of several CO 2 tracer pools (Table 1). CO 2 transport is simulated within the model domain concurrently with meteorological variables (feedback on atmospheric radiative transfer due to variable CO 2 is neglected). The CO 2 tracer scheme is a modified version of the scheme used in WRF-VPRM (Vegetation Photosynthesis and Respiration Model) (Ahamdov et al., 2007). Atmospheric CO 2 fields were provided by Carbon Tracker Europe (CTE, Peters et al., 2010) TM5 providing 1 • × 1 • resolution fields at 3 h intervals. CTE CO 2 fields were used to provide WRF-SPA CO 2 initial conditions (IC) and lateral boundary conditions (LBC) were linearly interpolated to the WRF-SPA domain. LBC for the outer domain have been set with zero inflow and zero-gradient outflow for all CO 2 fields, except total atmospheric CO 2 and "forcings only" Global flux maps of anthropogenic emissions and ocean absorption, also from CTE, at 1 • × 1 • resolution with 3 h updates were used to provide non-biospheric surface exchange. The fluxes were interpolated using 4 point weighted means based on latitude and longitude co-ordinates. Biospheric fluxes of CO 2 are simulated by the LSM (SPA, described in Sect. 3). All surface CO 2 fluxes were calculated as rates which were added to the lowest model atmospheric layer of the WRF grid at each time step.

Model description: SPA
The Soil-Plant-Atmosphere model is a high vertical resolution mechanistic point model (up to 10 canopy layers and 20 soil layers). SPA uses coupled energy, hydrological and carbon cycles to provide surface fluxes of heat, water and CO 2 to WRF. SPA provides realistic responses to meteorological drivers by coupling its hydrological and carbon cycles through ecophysiological principles (Williams et al., 1996).
SPA has been extensively validated against eddy covariance observations over several ecosystems including temperate deciduous forests (Williams et al., 1996), Arctic tundra (Williams et al., 2000), temperate evergreen forests (Williams et al., 2001) and, with the addition of a crop development model, temperate crop systems (Sus et al., 2010). SPA has been coupled to PBL models (Lee and Mahrt, 2004;Hill et al., 2008). Hill et al. (2008) successfully demonstrated that SPA could include feedbacks and drive PBL development that agreed with radiosonde observations. A brief description of the SPA model will be given here, followed by a detailed description of the modifications made to the SPA for use with the WRF; a detailed description of major SPA developments can be found in Williams et al. (1996Williams et al. ( , 1998Williams et al. ( , 2001Williams et al. ( , 2005 and Sus et al. (2010). A complete parameter list is available in the Supplement.
WRF provides SPA with meteorological drivers from the lowest atmospheric model level including air temperature, precipitation, vapour pressure deficit (VPD), wind speed, friction velocity, atmospheric CO 2 , air pressure, short and long-wave incoming radiation. SPA currently has parameters for 8 vegetation types (evergreen forest, deciduous forest, mixed forest, crops, managed grassland, grassland, upland and urban) suitable for UK application and 13 soil types. Vegetation and soil classifications are from the WRF default land cover maps (Mesoscale and Microscale Meteorology Division, 2011).
Plant phenology and carbon dynamics are described by a box carbon model (the Data Assimilation Linked Ecosystem Carbon (DALEC) model), which is fully integrated into SPA, to simulate the main ecosystem C pools (Williams et al., 2005). C pools (foliage, structural/wood carbon, fine roots, labile, soil organic matter (SOM) and surface litter) were "spun-up", in an offline SPA simulation (except for crops) using 3 yr of meteorology (1998)(1999)(2000) from Griffin Forest. These observations are broadly representative of the Scottish average and are from a period not simulated here.
The observations were obtained from the CarboEurope network (www.carboeurope.org/) and looped for a 30 yr period. A 30 yr period was found to be sufficient for carbon pools to reach steady state when SOM was initialised with realistic values for Scotland based on the soil carbon stocks from Bradley et al. (2005). No spin up of the above ground vegetation in crops was needed as arable crop systems are annual with complete clearing of the biosphere at harvest and addition of labile carbon in the form of seed. DALEC (Data Assimilation Linked Ecosystem Carbon) provides a direct coupling between the carbon cycle and plant phenology, specifically foliar and fine root C, where foliar C determines leaf area index (LAI) and root C impacts water uptake potential. Crops have two additional C pools; storage organ C (i.e. harvestable C) and dead foliar C (still standing) (Sus et al., 2010). Urban cover is assumed to be a low density evergreen forest with a reduced emissivity to be consistent with urban construction materials, all other surface properties remain unchanged; WRF-SPA used the same emissivity value for urban cover as used in the default WRF LSM, Noah. However, we expect this parameterisation to have little impact as urban cover represents < 1 % of the modelled land surface.
The Farquhar model of photosynthesis (Farquhar and von Caemmerer, 1982), the Penman-Monteith model of leaf transpiration (Jones, 1992) and the leaf energy balance are coupled via a mechanistic model of stomatal conductance. SPA maximises carbon assimilation per unit nitrogen but within a minimum leaf water potential to prevent cavitation via a series of bisection procedures. Leaf water potential links atmospheric demand for water and soil water supply, by including the effects of soil and stem hydraulic resistance on water transport to the leaves (Williams et al., 1996).
The soil surface energy balance is solved following the approach by Hinzman et al. (1998) and the soil temperature profile is updated by an implicit method of Crank-Nicolson (Farlow, 1993). Soil hydrology is calculated through determining soil surface evaporative flux and water movement within the soil profile due to gravity, root uptake and thermal distribution through the profile. Soil hydraulic parameters are calculated using equations from Saxton et al. (1986).
SPA uses a detailed radiative transfer scheme which models the absorption, transmittance and reflectance of near infra-red (NIR), direct and diffuse photosynthetically active radiation (PAR) and long-wave radiation for both sunlit and shaded fractions of each canopy level (Williams et al., 1998). Albedo is calculated from the overall reflectance and absorption of NIR and PAR from the canopy and soil surfaces.
The biological components of SPA remain unchanged in WRF-SPA. Modifications to the physical processes include updates to the canopy interception of precipitation, water storage and drainage calculations; a new aerodynamic scheme for momentum decay above and within the canopy; leaf level conductance calculating both free and forced convective exchange; an integrative procedure for calculating turbulent exchange of soil surface through the canopy; inclusion of dew formation and wet surface canopy evaporation within the canopy energy balance and addition of dead foliage LAI and post harvest litter within the radiative transfer scheme. Modifications are described in turn, below.

Canopy hydrological parameters
Canopy interception, water storage and drainage have a significant impact on potential wet canopy evaporation, dew formation (and therefore on the canopy energy balance), and soil surface water. Canopy interception of precipitation (I ; fraction) and maximum canopy water storage (C max ; mm) are related to LAI by coefficients α (0.5) and µ (0.2), respectively. α has been selected to generate canopy interception fractions which are consistent with Rutter et al. (1975). Canopy storage coefficient µ is intended to calculate values which are consistent with canopy storage values used in a previous SPA study (Williams et al., 2001). and Canopy drainage rate is calculated using an empirical relationship derived by Rutter et al. (1975), with an LAI adjustment factor.
where D is drainage rate (mm min −1 ), b is an empirical coefficient (b = 3.7), C stor is the current canopy storage of water (mm) and a accounts for the canopy water content relative to C max .
D c is the rate of drainage on a canopy where C stor = C max . D c is adjusted via a proportional relationship to C max (Rutter et al., 1975).
3.2 Aerodynamic scheme: canopy exchange SPA models both the above and within canopy momentum decay. Wind speeds are used in determining leaf level and soil surface conductance (described later). Above canopy momentum decay follows the standard log law decay with the Monin-Obukov similarity theory stability correction (Garratt, 1992).
where dU/dz is the gradient of wind speed decay above the canopy at height z (m), κ is Von Karman constant (0.41), d is the canopy zero plane displacement height (m), u * is the friction velocity (m s −1 ) and m is the Monin-Obukov stability correction coefficient. The gradient of wind speed decay is integrated over the vertical distance between the wind speed at the reference height to the canopy top. It is important to note that currently the roughness sub-layer is not included in decay calculations. The displacement height and roughness length (z o ) are calculated based on canopy structure (height z h and LAI) as described in Raupach (1994). and where z h is the canopy height (m), C d1 is an empirically fitted parameter (7.5) and h parameterises the effect of the roughness sub-layer on roughness length (0.193) (Raupach, 1994). Within canopy momentum, decay is carried out using the method described in Harman and Finnigan (2007). Decay within the canopy is assumed to be exponential, where U (z) is the wind speed (m s −1 ) at height z within the canopy. Decay is dependent on the canopy mixing length (l m ) and ratio of u * u h = u r (where u h is the wind speed at canopy top).
The canopy mixing length is described by u r and the canopy length scale L c (m), where the length scale is calculated assuming a uniform canopy.
The resulting within-canopy wind speed profile is used in the calculation of canopy layer specific boundary layer conductance (m s −1 ) of heat and water vapour. Boundary layer conductance at each canopy layer is assumed to be the maximum conductance between free and forced convection at that layer. A detailed description of the leaf level conductance is given in Nikolov et al. (1995): where g h is the leaf level conductance for heat (m s −1 ), D h is the molecular diffusivity of heat (m 2 s −1 ), S h the Sherwood number and d o is the leaf or needle (cone if free convection) diameter. For calculation of water vapour conductance, g wv , the molecular diffusivity of water, D wv , is used.

Aerodynamic scheme: soil exchange
The bare soil surface conductance (g soil ; m s −1 ) for heat and water vapour are assumed to be equal.
where U soil is the wind speed near the soil surface, and z ref is the reference height of the lowest model level to which the soil is exchanging. The soil surface roughness length (z soil ; m), is assumed to be equal to 0.01 m both when under a canopy and as bare ground. When the soil is under a canopy, the soil conductance is first calculated as a resistance. Soil resistance is integrated through the canopy based on the turbulent eddy diffusivity following Niu and Yang (2004).
where dz is the vertical step size (m) through the canopy and K h is the eddy diffusivity at z position (m) within the canopy. Eddy diffusivity (K h ; m 2 s −1 ) is assumed to have an exponential decay through the canopy (as with momentum). Eddy diffusivity at the canopy top is estimated as specified in Kaimal and Finnigan (1994).
K h is decayed through the canopy as described below. The coefficient of momentum decay f is dependant on c d the coefficient of drag for foliage (0.2), LAI, l m and soil surface m . m was calculated as described in Garratt (1992). When ζ > 0, when ζ < 0, where ζ > 0 conditions are unstable, while ζ < 0 are stable conditions, with coefficients γ = 16 and β 1 = 5. ζ = z/L is described in Qin et al. (2002).
where L is the Obukov length, g is acceleration due to gravity (9.81 m s −2 ), cp air is the specific heat capacity of dry air (1004.6 J kg −1 K −1 ), ρ is the density of the air (kg m −3 ). H soil is the sensible heat flux from the soil in the previous time step.

Leaf energy balance
SPA uses an iterative procedure to solve stomatal conductance, and as part of this procedure, the leaf energy balance is solved. Net radiation is partitioned between latent and sensible heat fluxes (W m −2 ); the metabolic storage term is assumed to be small and is neglected. Evaporation is calculated from the Penman-Monteith equation; sensible heat is based on the temperature difference between the leaf surface and surrounding air. The radiative transfer scheme initially distributes long-wave radiation assuming that the canopy and surrounding air are in isothermal net radiation balance (Rni; W m −2 ). Rni is updated to the net radiation (Rn; W m −2 ) in the first iteration.
Rn Rni + 4 σ T 3 a ( T ) The Rn correction is based on the temperature difference ( T ) between the input leaf temperature T leaf , in the first iteration, and the absolute air temperature (T a ; K). T leaf is solved by balancing the canopy energy balance Eq. (22).
E leaf is transpiration (kg m −2 s −1 ), δc w is the absolute humidity deficit (kg m −3 ) and = s/γ . s (Pa K −1 ) is the slope of curve that relates saturation vapour pressure with air temperature and γ is the psychrometer constant (Pa K −1 ). g s is stomatal conductance (m s −1 ), and λ is the latent heat of vaporisation (J kg −1 ).
H leaf is sensible heat (W m −2 ), leaf and air are potential leaf and air temperatures respectively. The factor 2 is to account for the two different sides of the leaf.
Wet canopy evaporation is calculated via a multi-stage process. The Penman-Monteith equation is used to calculate the potential evaporation or dew formation (E pot ; kg m −2 s −1 ), i.e. with aerodynamic conductance only.
δe is the vapour pressure deficit (Pa). If dew is forming this is restricted to no larger than C max . Dew mass is added to the canopy in the following time step, while the energy exchange is added to Rn for the next iteration of the canopy. Wet evaporation (E wet , kg m −2 s −1 ) is restricted by the amount of water stored on the canopy (C stor , Eq. 26) and E leaf that has occurred Eq. (26).
An iterative solution is used to model both wet canopy evaporation, (or) dew formation and update the canopy distribution of long-wave radiation. (i) Leaf temperature is calculated solving the canopy energy balance; (ii) long-wave radiation distribution is updated based on current leaf temperature values; (iii) leaf temperature is re-solved; (iv) canopy net radiation is used to calculate wet evaporation or dew; (v) the available net radiation for the canopy energy balance is adjusted based on wet evaporation or dew. The procedure is iterated up to 10 times; steady state of the canopy balance typically occurs by the 4th iteration.

Modifications to the crop model in SPA
The crop development model developed for SPA is described and validated in Sus et al. (2010). This allows SPA to model the growth of both winter wheat and winter barley. While only modelling winter cereals will introduce a bias, it is a reasonable assumption as winter wheat and barley are the dominant arable crops in the UK. The model represents crop development from sowing through vegetative growth, flowering and maturity through to harvest. Sowing is assumed to occur once the daily mean temperature has dropped below 10 • C for winter crops. Van den Hoof et al. (2011) demonstrated that this assumption produces realistic sowing times. Harvest is assumed to occur once the developmental crop model reaches the mature development stage, which typically coincides with the storage organ (the crop yield) reaching its peak value and complete foliage senescence.
The albedo of crops changes as it matures, this change leads to a significantly greater albedo in mature crops compared to during vegetative growth. Seasonal shifts in albedo significantly alter the surface energy balance, changing both surface temperature and turbulent fluxes .  To account for this, foliage which has senescenced is retained within the canopy as a non-photosynthetically active (non-transpiring) dead LAI. In addition to decoupling dead foliage from the plant hydrological cycle, dead LAI is assigned its own NIR and PAR reflectance values (Nagler et al., 2003). Assigning dead LAI reflectance allows for inclusion of their effect on the surface energy balance, impacting both leaf level processes and radiative transfer. Post harvest, litter also plays an important role in the surface energy balance by increasing surface albedo relative to the underlying soil. As with dead foliage, post harvest surface litter is prescribed a separate albedo which is weighted with that of the underlying soil, based on fractional area cover. Litter fractional cover is estimated based on the mass of surface litter (C sflit ) present. Surface litter from both wheat and barley are assumed to have the same mass-area relationship as wheat from Nagler et al. (2003).

Model domain
WRF-SPA was run over two domains with two-way nesting; the outer domain has a resolution of 18 km × 18 km and the inner 6 km × 6 km (Fig. 1). Model output from the inner domain only was used in the validation. Scotland provides a highly complex topography and land use heterogeneity, with a longitudinal gradient from dominantly forested areas in the northwest to pasture in the central and southwest and arable cropland in the east. A 5 yr period (2002)(2003)(2004)(2005)(2006) was simulated for use in a multi-annual validation of the model at different spatial scales, from surface measurements to vertical aircraft profiles of CO 2 atmospheric concentrations. The first 2 yr were considered to be a spin up period to allow for differentiation of the vegetation phenology. The main features of the model set-up are presented in Table 2.

Validation data
The surface validation used surface observations of net radiation, turbulent fluxes (latent and sensible), net ecosystem exchange of CO 2 (NEE) and air temperature from three sites in Scotland. The sign convention used for NEE is negative fluxes represent net sequestration of carbon and positive fluxes represent a net source. Observations used were averaged to an hourly time step. The three sites are important as they are representative of the dominant land cover types in Scotland outside of the central northern mountain ranges (Fig. 1).
The sites are (i) Griffin Forest, an intensively managed Sitka spruce plantation (LAI ∼ 6 m 2 m −2 ) in central Scotland (56.61 • N, 3.80 • W, 340 m a.s.l.). Established in 1982, the site has an mean annual air temperature of ∼ 6.6 • C and precipitation of ∼ 1126 mm yr −1 (Clement et al., 2012). Air samples were collected at a range of altitudes above sea level (m a.s.l.), 800, 1100, 1600, 2100, 2600 and 3100. Sampling (all daytime) occurred throughout each year covering the whole seasonal cycle. This provides information at all development stages of the vegetation and includes a range of meteorological conditions. The profiles provide integrated regional-scale information allowing for validation at the regional-scale of WRF-SPA's carbon balance.

Surface validation
Statistical validation (analysed in R v2.15.2; R Core Team, 2012) against hourly observations are presented for WRF-SPA and WRF (Table 3). Seasonal behaviour is explored by comparing monthly mean (and standard error) values for observations and the models. To ensure comparability, the model means are calculated using only values where a corresponding observation is available. Both models demonstrate good skill in predicting surface observations, particularly air temperature. R 2 and biases are similar (Table 3); on average the differences between the modelled biases are 0.16 • C and 4.2 W m −2 for temperature and turbulent fluxes, respectively. There is a tendency for latent and sensible fluxes to be positively biased, except at Griffin Forest where WRF underestimates latent heat. The overestimation of turbulent fluxes is in part due to overestimation of net radiation, typically 10-15 W m −2 for both models (Table 3). However errors in partitioning to ground heat flux are also likely to play a significant role. In WRF-SPA, ground heat flux is relatively insensitive seasonally; underestimating the magnitude of ground heat flux may lead to a bias in energy partitioning to turbulent fluxes during summer (data not shown). WRF-SPA tends to have larger RMSEs for latent (5-20 %) and sensible heat fluxes (5-92 %) than WRF, while for temperature RMSEs are similar with a maximum difference of 0.2 • C. The largest RMSEs are for East Saltoun and are discussed later.
Air temperature is well predicted by both models with a typical annual absolute bias of < 1 • C and RMSE of ∼ 2 • C. This is broadly consistent across all sites at the annual timescale, representing skill in modelling a range of highly distinct vegetative systems. At seasonal scales, more significant differences are apparent between the two models. WRF produces a consistent overestimation of monthly mean temperatures at Griffin and East Saltoun, while during winter both models overestimate air temperature by a similar amount (Fig. 2). WRF-SPA more closely predicts monthly mean temperatures during the majority of the year, except at Easter Bush during the spring and summer months where WRF-SPA underestimates temperature. At Griffin Forest both models realistically predict seasonal behaviour, while at East Saltoun and Easter Bush WRF-SPA more accurately predicts the observed seasonality than WRF. Including accurate prediction of peak summer monthly mean temperature in August for both East Saltoun and Easter Bush during 2004. NEE is reasonably well predicted by WRF-SPA at each site, particularly at the seasonal timescales (Fig. 2). However hourly RMSEs are high, typically ∼ 6 µmol CO 2 m −2 s −1 while biases remain relatively small, usually < 1 µmol CO 2 m −2 s −1 (Table 3). Seasonally, each of the sites show different biases throughout the year from each other. WRF-SPA underestimates net carbon sequestration at both Griffin Forest and Easter Bush over the observation period, including during winter indicating an overestimation of respiration in the model. The WRF-SPA model-observation mismatch at Griffin Forest varies between years. WRF-SPA overestimates peak sequestration (NEE is more negative) during 2005 and fails to capture a peak summer reduction in carbon sequestration shown in the observations (Fig. 2), whereas in 2006 WRF-SPA captures the seasonal behaviour of the forest but there is some underestimation in fluxes (Fig. 2).
At East Saltoun, NEE appears to be poorly modelled at the hourly level with an R 2 of 0.32. While at seasonal timescales we can see that WRF-SPA models NEE reasonably well, however the phenology is out of phase (Fig. 2). Despite this the growing-season peak sink strength and post-harvest respiration peaks are of appropriate magnitude. In both years WRF-SPA predicts peak sequestration two months early and overestimates early season sequestration. The model observation mismatch is consistent with WRF-SPA modelling winter barley at East Saltoun, while the crop planted during 2004 and 2005 is spring barley.
WRF-SPA overestimates latent heat fluxes at all three sites, particularly during the summer when more energy is available. WRF overestimates at both East Saltoun and Easter Bush, and underestimates latent heat at Griffin Forest (Fig. 3). Seasonally we can see that WRF-SPA performs better, particularly in modelling the transitions into and out of summer. This is most evident during 2004 at East Saltoun where, while the peak latent heat flux is too high, the peak period timing and duration is considerably better represented than by WRF. A similar situation is observed for sensible heat (Fig. 3) with particularly good agreement at Griffin Forest by WRF-SPA while less so by WRF. Seasonal transition periods are well captured (Fig. 3), which is expected due to the coupling between latent and sensible heat fluxes. WRF-SPA captures the bimodal peaks during both years, which is at least partly driven by harvest. When the crop is harvested this removes much of the high albedo mature vegetation and exposes the low albedo soil, resulting in a higher net radiation. Removal of crop vegetation also restricts evaporation to water available in the upper soil layers. This has a significant impact on partitioning of net radiation into sensible and ground heat fluxes.

Scaling to aircraft profiles
Aircraft profiles were compared to WRF-SPA modelled profiles of CO 2 to provide validation of regionally integrated measurement of CO 2 exchange (Fig. 4). Profiles include both summer and winter flights throughout the simulation period, to investigate model performance both seasonally and multiannually. Figures show two modelled profiles, (i) total atmospheric concentration and (ii) "forcings only" profiles. The "forcings only" pool, contains CO 2 which originates from external forcings only, including LBC nudging, oceanic fluxes and anthropogenic emissions (i.e. total CO 2 minus the modelled biospheric fluxes). This allows us to visually show the impact of the "simulated biosphere" on atmospheric CO 2 concentrations within the domain.
The modelled profiles of total CO 2 compare well with observations. All predicted profiles are typically within Fig. 4. Observed and WRF-SPA modelled profiles of CO 2 above Griffin Forest. Modelled profiles include the total atmosphere profile CO 2 and "forcings only" CO 2 to show the impact of the modelled biosphere.
2-4 ppm of observations, a selection of which are shown in Fig. 4. The profile structure and PBL heights of modelled total CO 2 profiles compare well with observations indicating that regional transport and distribution of biospheric sources and sinks is broadly realistic. "Forcing only" CO 2 profiles show distinctly different structures and CO 2 concentrations within the PBL, indicating that the good agreement we find in the total CO 2 pool is due to the modelling of the biosphere by WRF-SPA.
There are differences in model performance across seasonal scales. During peak growing season, WRF-SPA overestimates regional biosphere sequestration, leading to up to 4 ppm underestimation in the modelled profile (e.g. July 2004 and June 2006; Fig. 4). This is different to the observed underestimation of carbon sink strength at Griffin Forest (Table 3). The likely reason for this is that forest cover is overrepresented within the model domain. At the resolution used in WRF-SPA, the MODIS land cover map used in WRF-SPA estimates that forest-type land covers make up 43 % of the total land cover. The MODIS forest cover estimate is significantly greater than the Scottish estimate of ∼ 17.8 % forest cover (National Forest Inventory, 2011). This issue could be rectified through replacement of the MODIS map with a more realistic representation of UK land cover (e.g. LCM 2007, http://www.ceh.ac. uk/LandCoverMap2007.html).

Discussion and conclusions
The purpose of this study is to present and validate the novel coupled model WRF-SPA. Validation of WRF-SPA was conducted at both local (site level) and regional (aircraft profiles) scales, against a range of observations relevant to CO 2 exchange and validation of the meteorological capability of the WRF-SPA.
WRF-SPA demonstrated that it can produce comparable statistics to those of WRF (Table 3) at the site level against hourly observations. WRF-SPA tends to have lower annual bias, however WRF-SPA has higher RMSEs. Higher RMSEs are not unexpected due to the greater level of complexity of SPA compared to the default WRF LSM, Noah. WRF is often considered to be one of the best mesoscale models available in terms of simulating vertical profiles of temperature (Steeneveld et al., 2011) and surface meteorological variables . Therefore, we can infer that WRF-SPA is also comparable to many other models that are currently used in regional-scale research.
At seasonal timescales WRF-SPA shows realistic behaviour across each of the land cover types presented here. Air temperature was consistently better modelled by WRF-SPA at both Griffin Forest and East Saltoun, while summer peak temperatures were better captured at Easter Bush by WRF-SPA. NEE is well modelled by WRF-SPA, particularly seasonality (Fig. 2). There remain, however, several issues including an overestimation in winter respiration at Griffin Forest and Easter Bush, which may be linked to an overestimation of soil carbon stocks or an overestimation of soil organic matter turnover. Given Scotland's high soil organic matter content it is likely that WRF-SPA's modelling of soil processes as inorganic is a significant component of this error.
There remain several model-observation mismatches at Griffin Forest during the growing season (Fig. 2). Where WRF-SPA overestimates carbon sequestration in 2005 but underestimates sequestration in 2006. Further, there is presence of a lag in simulated NEE, particularly evident in 2006. Underestimating sequestration in 2006 is consistent with WRF-SPA underestimation of LAI at Griffin Forest (modelled ∼ 4 m 2 m −2 , observed ∼ 6 m 2 m −2 ). Griffin Forest underwent selective logging in 2004 (∼ 37 % reduction in above ground biomass), therefore during 2005 the forest was likely in the process of recovery before returning to preharvest sequestration levels in 2006. The early season lag is a known issue with the evergreen carbon model in SPA, were the evergreen model lacks a labile carbon pool. As a result needle growth is limited in springtime to carbon available from photosynthesis at that time, rather than rapid release of carbon stored from the previous year (Williams et al., 2005).
WRF-SPA overestimates latent heat at all three sites compared to observations (Fig. 3). However the model-data mismatch may be due to an underestimate of turbulent fluxes in observations due to non-closure of the surface energy balance (Stoy et al., 2013). Stoy et al. (2013) also found evidence that non-closure of the energy balance is greater in crop systems, which is consistent with the greatest model-data mismatch seen as East Saltoun.
Seasonal changes in latent and sensible heat fluxes are well predicted by WRF-SPA compared to WRF. Transition periods in particular from winter to spring and summer to autumn are well captured by WRF-SPA compared to WRF, e.g. latent heat fluxes observed at East Saltoun during 2004 (Fig. 3). As the bimodal peak in sensible heat flux at East Saltoun appears to be driven by human intervention on the ecosystem, demonstrating it is important to include the impacts of human management of ecosystems. This is consistent with other studies which have previously highlighted agricultural land as a significant component of the European and global energy and carbon balance Denman et al., 2007). WRF also displays a bimodal peak during 2006, the cause of this is unclear but appears to be linked to a larger amount of incoming short-wave radiation in WRF than WRF-SPA during the mid-to late-summer. The differences in incoming radiation between the models is likely as a result of combination of feedbacks where in WRF-SPA the higher evaporation rate and sensible heat flux lead to increased cloud cover and atmospheric albedo.
Comparison with aircraft profiles show that WRF-SPA is capable of upscaling to regional measurements of atmospheric CO 2 . Indicating that the modelled CO 2 sources/sink T. L. Smallman et al.: WRF-SPA: development and validation distribution is well modelled as is transport within the atmosphere. The absolute errors of the modelled profiles are comparable to other studies (e.g. Ahamdov et al., 2007;Ter Maat et al., 2010). What is unique about our study is that we carried out this analysis over a several year period showing multi-annual consistency. WRF-SPA successfully predicts PBL height for the observed profiles which is critical to achieve realistic regional-scale mixing and meteorological variables (Steeneveld et al., 2011;Xie et al., 2012). Steeneveld et al. (2011) showed in a comparison between different PBL parameterisations, using WRF and RAMS, that correct modelling of PBL structure and temperature profiles was only achieved with significant overestimation of sensible heat fluxes of up to 50 %, which is consistent with the overestimate of sensible heat fluxes predicted by WRF-SPA.
WRF-SPA has demonstrated that it is capable of realistically predicting exchange of CO 2 at the regional scale. Further work will look at using CO 2 tracers to de-construct observations of regional exchange to investigate how each vegetative ecosystem contributes to the regional signal. In particular a currently unexplored multi-annual dataset of continuous measurement of CO 2 from the tall tower Angus site, Scotland.
In conclusion, three specific questions were asked of this validation. (i) Can WRF-SPA realistically model surface meteorological variables and fluxes across a multi-annual period? We have shown that WRF-SPA can realistically model surface observations at hourly timescales which are comparable to WRF. (ii) Does WRF-SPA scale realistically from the surface measurements to regional-scale observations, specifically aircraft profiles? We have shown across multiple years and across seasons that WRF-SPA upscales to observed atmospheric CO 2 concentration profiles. (iii) Does WRF-SPA lead to an improvement in surface fluxes compared to the unmodified WRFv3.2? We have demonstrated that at monthly means WRF-SPA predicts more realistic seasonal behaviour than WRF.
Acknowledgements. The authors would like to thank the PhD project funding body, the National Centre for Earth Observation, a UK Natural Environment Research Council research centre. Christoph Gerbig of the Max Plank Institute is thanked for providing the original CO 2 tracer modifications for WRF. John Finnigan of CSIRO is thanked for help in implementing a new canopy momentum decay parameterisation within SPA. The aircraft observations were funded by the AEROCARB project of the EU Framework Programme 5.
Edited by: R. Sander