Do we (need to) care about canopy radiation schemes in DGVMs? Caveats and potential impacts

. Dynamic global vegetation models (DGVMs) are an essential part of current state-of-the-art Earth system models. In recent years, the complexity of DGVMs has increased by incorporating new important processes like, e


Introduction
Land surface models are one of the required tools for understanding land surface dynamics, land-atmosphere interactions and climate-carbon feedbacks, and are an essential part of dynamic global vegetation models (DGVMs).DGVMs are widely used to assess climate change impacts on vegetation distribution and terrestrial carbon, water and energy fluxes and feedbacks of the biosphere in the Earth system as well as climatic consequences of land cover change (Friedlingstein et al., 2006;Brovkin et al., 2013a).Given the major policy implications for climate change mitigation, much attention is placed on the performance and realism of these models, resulting in an increased overall complexity.At the same time, among land surface models, there is no consensus on important aspects of the carbon cycle in a future climate (Sitch et al., 2008).
Published by Copernicus Publications on behalf of the European Geosciences Union.
DGVMs typically represent the land surface dynamics in a simplified manner.For practical reasons they represent small-scale processes using large-scale variables.In doing so, different models make different approximations in their representation of processes.The relevance and realism of several approximations and their potential implications for the range of projections found have been discussed elsewhere as far as it concerns the ecology (Harrison et al., 2010;Van Bodegom et al., 2012) and soil carbon dynamics (Ostle et al., 2009).However, also in their representation of exchanges of energy (radiation and heat), models differ widely and occasionally have been implemented differently for albedo and fraction of absorbed photosynthetic active radiation (faPAR) calculations.Most models are confined to one-dimensional (vertical) exchange of radiation, mostly relying on solutions derived from two-stream approximations based on plane-parallel turbid media assumptions, like those of Sellers (1985).
The transfer of radiation within canopies is complicated by multiple scattering of radiation, mutual shadowing, variations in leaf orientations and crown closure, as well as variable optical properties.Sprintsin et al. (2012) have shown that neglecting, e.g., the differences between sunlit and shadowed leaves in canopy radiative transfer (RT) schemes results in a significant underestimation of the canopy-level gross primary production (GPP).A similar study analyzing the impact of a revised canopy RT scheme in the ISBA land surface model has been documented by Carrer et al. (2013).In addition, full 3-D radiative transfer models exist for canopies as well as landscapes (Gastellu-Etchegorry et al., 2004;Kobayashi and Iwabuchi, 2008), but these cannot be used directly in DGVMs due to their high computational demand and the high number of required (vegetation structural) parameters.Numerically fast 1-D models with a limited number of required input parameters are therefore still preferably used.
These 1-D models simplify radiative complexities caused by vegetation clumping, representing the concentration of vegetation and thus canopy scattering and absorption within a given area.This leads to so-called effective variables (Widlowski et al., 2005;Rochdi et al., 2006;Pinty et al., 2006).However, neither 1-D approach has been incorporated consistently within the currently state-of-the-art DGVMs.Alternative 1-D models have been developed that assume largescale canopy elements as having simple shapes (e.g., spherical crowns) and give analytical solutions for the calculation of surface fluxes using a 1-D model (Dickinson et al., 2008;Haverd et al., 2012).
Given these complications, there is increasingly a call to systematically evaluate and benchmark DGVMs (Abramowitz et al., 2008;Luo et al., 2012;Hagemann et al., 2013;Brovkin et al., 2013b).Benchmark analysis is essential for identifying uncertainties in predictions as well as for guiding priorities for further model development (Blyth et al., 2011).Despite the importance of heat and energy exchange in DGVMs as major drivers of surface temperatures and carbon productivity, current benchmark analyses do not yet account for energy budgets.Where benchmark initiatives consider albedo and faPAR, separately -and thus from an RT point of view, potentially inconsistently -Earth observation products are proposed as candidate benchmarks (Luo et al., 2012;Hagemann et al., 2013).
There is thus a need to assess the consistency and accuracy of radiative transfer schemes in DGVMs and to assess the potential impact of uncertainties in the widely used radiative transfer models on surface energy fluxes and carbon production estimates, which so far have not been quantified for widely used DGVMs.Such assessment requires physically consistent three-dimensional radiative transfer formulations as a reference.Such reference simulations are provided by, e.g., the Radiative Transfer Model Intercomparison Initiative (RAMI) (Pinty et al., 2001;Pinty, 2004;Widlowski et al., 2007).
Here, we aim at a show case to (i) assess the capabilities of canopy RT formulations, used in state-of-the-art DGVMs, to simulate consistently canopy absorption and reflectance for idealized reference cases, (ii) assess at which conditions the used canopy RT schemes (and their simplifications) might lead to major biases in faPAR and/or albedo estimates, and (iii), importantly, to assess the potential implications thereof for net irradiance and carbon productivity estimates.
The present study therefore does not aim at global-scale evaluation of DGVMs.Global-scale model benchmarking activities mainly focus on the evaluation of models in simulating the spatial and temporal patterns of land surface dynamics at climatological timescales (Brovkin et al., 2013b;Hagemann et al., 2013).The current study is complementary to these global scale evaluation studies as it explicitly assesses the accuracy of different canopy RT schemes under different, well-defined conditions.It therefore aims at an evaluation of the models on a process scale, rather than looking at large-scale, global impacts.
The study does not however provide a full assessment of the impact of deficits in canopy RT schemes on global carbon and energy fluxes in coupled Earth system models on longer timescales, as this would require a much more comprehensive analysis than the one provided in the present study.
Thus the overall objective of this paper is to raise awareness of the relevance of canopy radiation and surface albedo schemes in globally applied land surface schemes being used in a wide range of applications.

Representation of canopy RT in DGVMs
Global DGVM simulations are typically performed on coarse spatial resolutions with model grid cell sizes on the order of 10 2 to 10 8 km 2 and long time periods (decades to millenia).To represent surface processes -including radiative transfer -at the sub-grid scale, a tiling (mosaic) approach is widely used where surface processes are simulated on tiles of N individual plant functional types (PFT) and where results are combined to a model grid box average by area-weighted averaging.Plant functional types provide a means of using a finite set of model parameters to simulate plants with similar ecological behavior (Diaz and M., 1997;Prentice et al., 2007).The analysis of the present study focuses on the performance of canopy RT models and assumes that a grid cell is covered by a single land cover type (N = 1).
The applications of DGVMs demand a numerically stable and fast canopy radiative transfer scheme.Simple 1-D canopy RT schemes or parametric approaches are therefore used in DGVMs instead of more computationally demanding 3-D RT schemes.A 1-D model is understood here as a model where the surface state variables (e.g., leaf area, leaf reflectance and absorption properties) vary only along a single coordinate axis (typically the z direction).Note that this is not related to the spatial integration needed for flux calculations.When integrating over the upper and lower hemispheres, the model is called a two-stream model.However, angular integrations can also be made only for parts of spheres, resulting then in a four-stream model for half spheres and an M stream solution for M numbers of subspaces.If M is large, a bidirectional reflectance model (BRF) is obtained (Gobron et al., 1997).
A turbid medium assumption is widely used for 1-D canopy simulations, where foliage elements are assumed to be point like scatterers -leaves are assumed to have an infinitely small size -randomly and uniformly distributed, and typically stems and branches are neglected.The gaps within canopies can be expressed by the gap probability (P gap ), which is defined as the probability of a beam at the sun zenith angle θ hitting the ground without interacting with canopy elements.In other words, it is the ratio between the uncollided radiation flux and the incident radiation (Haverd et al., 2012) and represents the domain-averaged direct transmission.In the case of a horizontally homogeneous canopy with randomly distributed leaves (and no foliage), P gap follows a Poisson distribution and is given by where µ = cos(θ ), is the leaf area index (LAI) (m 2 m −2 ) and G(θ ) is the mean projection of unit leaf area in the direction perpendicular to the incoming beam.Often a spherical leaf angle distribution is assumed, which results in G(θ ) = 0.5.This representation of P gap is widely used in DGVMs and applies to all models investigated in this study.While Eq. ( 1) is a good approximation of P gap for closed canopy cases, it has been recognized that it overestimates canopy absorption for open (non-uniform) canopies (Pinty et al., 2006;Chen et al., 2008;Haverd et al., 2012).To account for non-uniform distributions, a so-called clumping factor 0 < < 1 has been introduced for use in a one-dimensional radiative transfer model.The gap probability is then calculated as where ˜ (µ) is an effective leaf area index that reproduces the correct canopy RT fluxes using a 1-D formulation instead of a 3-D model (Pinty et al., 2006).The fraction of absorbed photosynthetic radiation (faPAR) is a crucial variable in DGVM carbon flux simulations.Ignoring horizontal fluxes (Widlowski et al., 2006), faPAR is defined as the fraction of radiation absorbed by the canopy elements (foliage and woody parts) in the photosynthetic active electromagnetic spectrum (= VIS) and is defined by considering the energy balance as where PAR i is the incident down-welling PAR at the top of the canopy, PAR α is the reflected component, PAR T is the total transmission through the canopy and α s is the albedo of the soil.A first-order approximation to faPAR is given by assuming that the leaves and soil are completely black (i.e., all radiation incidents on them is absorbed).In this case PAR α = 0 and PAR T /PAR i = P gap , which yields: However, if the optical properties of the canopy elements are non-zero, then the exact solution becomes more complex as it must consider multiple scattering within the canopy, and between the crown and the ground if the soil albedo is also non-zero.Among the various currently used faPAR approximation schemes Eq. ( 4) was shown to provide the least systematic bias (Widlowski, 2010).
Ideally, a 1-D RT scheme would also consider the proportion of direct illumination and shading at different levels in the canopy, as this has potentially significant impacts on photosynthesis (Mercado et al., 2009).This entails splitting the incoming solar radiation into direct and scattered components (Spitters et al., 1986) and calculating the amount of energy intercepted at each canopy layer that is from the direct component that has not previously been scattered (i.e., the sunlit fraction), and the diffuse and multiply scattered components (i.e., the shaded fraction).Multiple scattering is typically neglected in simple canopy RT schemes employed in DGVMs.Two of the three models investigated in the present study consider however multiple scattering within the canopy and between the ground and vegetation layers.
Most DGVMs use a parametric 1-D RT scheme that is somewhere within the range of complexities described above.However, more complex 1-D schemes are available.In general the added complexity is in the form of extra detail in the description of canopy structure.from Monte Carlo ray tracing, which explicitly computes the path of a large number of photons through the vegetation canopy (Disney et al., 2000), to geometric optic techniques that treat the tree crowns as a distribution of geometric primitives and calculate the proportion of illuminated and vegetated parts of the scene (Li and Strahler, 1986).Ray tracing models are unlikely to be practical for directly embedding them in DGVMs, although they are useful tools for validating simpler models.Some models based on geometric optics are viable however.An example is the ACTS (Analytical Clumped Two-Stream) model (Ni-Meister et al., 2010), which is derived from the GORT (Geometric Optic Radiative Transfer) model (Ni et al., 1999).The approach used is to consider the canopy as a forest of randomly distributed spheroids filled with leaves.The probability of a photon entering a crown is calculated from geometric optic theory and then scattering is treated as a 1-D RT problem similar to those described elsewhere in this section.

Handling sparse vegetation (open canopies) in DGVMs
Despite the existence of concepts to account for vegetation clumping using geometric optics, more simplified approaches to handle sparse vegetated canopies have been implemented in DGVMs and are widely used.Figure 1 shows a model grid cell with unit area A, which is assumed to be fully covered by a particular PFT.The area A corresponds to the reference ground area for the leaf area index of the PFT.
In the case of sparse vegetation types like, e.g., savannas, the vegetated area is typically covered by a dominant plant functional type (e.g., trees), understory vegetation (e.g., shrubs, grass) or bare soil.In the process of converting land cover information, as is available from satellite products, a partitioning between these different components is made (Poulter et al., 2011).The total area A is thus defined as where f veg is the fractional coverage of the major PFT and f under can correspond to a different PFT or soil surface.It is important to recognize that f under is different from the gap probability P gap (θ = 0), as gaps within the dominant canopy are not considered in f under , therefore P gap (θ = 0) ≥ f under .
In other words, a common concept in DGVMs is to approximate P gap as where g( ) is a function that estimates the gaps within the canopy of the dominant vegetation type and might change over time t as a function of leaf area index , while f under is static (examples of the investigated models are given in the following section).Note that while P gap in Eq. ( 1) is a function of illumination conditions (θ), this is not the case for Pgap in Eq. ( 6).)) for a model grid cell that is assumed to be totally vegetated.

Models
The present study uses three different DGVMs, which are used in renowned Earth system models that all contribute to the Coupled Model Intercomparison Projects on a regular basis (Meehl et al., 2007;Taylor et al., 2012).These models use different canopy radiative transfer schemes that are representative of the canopy RT models widely used in DGVMs.It is therefore expected that results from this study can also be generalized to other models which are using similar canopy RT schemes.

JULES
JULES is the Joint UK Land Environment Simulator, a land surface model designed to predict the fluxes of heat, water and carbon between the land surface and the atmosphere.It originates from the Met Office Surface Exchange Scheme (MOSES) and is designed to be linked to the UK Met Office Unified Model.The fundamental equations underlying the model are common to many land surface schemes and are described in detail elsewhere (Best et al., 2011;Clark et al., 2011).In addition, JULES allows for a coupling with a general circulation model (GCM) and provides optional modules that allow longer term processes such as succession of plant functional types to be taken into account.These allow it to be used to simulate the response of the land surface to changing climatic conditions.In particular this includes the TRIFFID (Top-down Representation of Interactive Foliage and Flora Including Dynamics) dynamic vegetation model (Cox, 2001).The canopy radiative transfer scheme in JULES is based on the two-stream approximation proposed by Sellers (1985), which translates into consistent calculations of surface albedo and canopy absorption.Sellers (1985) where I ↑ (I ↓ ) are the upward and downward diffuse radiation fluxes normalized by the incoming radiative flux, τ = G(µ)/µ is the optical depth and ω = r + t the single scattering albedo.β and β 0 are upscatter parameters for the diffuse and direct beams, respectively, µ is the average inverse optical depth per unit leaf area and L corresponds to the cumulative leaf area index.For details on the solution of Eq. ( 7) see Sellers (1985).
JULES is the only one of the analyzed models that does a consistent simulation of absorbed and scattered radiation.Direct and diffuse radiation fluxes are computed individually.The leaf scattering and extinction properties are prescribed by the leaf reflectance (r), and single scattering albedo (ω) and multiple scattering within the canopy as well as between the soil and vegetation is considered.JULES uses a number of layers, as defined by the user, with equally distributed LAI density to simulate the canopy radiative transfer for the direct and diffuse flux components.Here an n = 20 layer model has been used, which is the default for JULES.A sensitivity analysis of the JULES canopy RT model revealed that the simulated surface reflectance is not sensitive to the number of layers.The simulated faPAR shows some sensitivity to the number of canopy layers.For n < 10, larger deviations in estimated faPAR occur for θ > 50 and large values of LAI ( > 4).Using n = 20 layers results in similar simulation results as using n = 100 layers and is therefore a good compromise between stability and numerical performance.
The scheme of Sellers (1985) used here is from JULES version 3.2.An earlier version (2.1) was used in the RAMI4PILPS exercise (Widlowski et al., 2011) and was found to give analogous results (not shown here).

JSBACH
The DGVM JSBACH (Raddatz et al., 2007;Brovkin et al., 2009;Reick et al., 2013) is implicitly coupled to ECHAM6, the atmospheric component of the Max Planck Institute for Meteorology Earth System model (Stevens et al., 2013).It simulates all relevant land surface water, energy and carbon fluxes.The present study uses version 2.03 of JSBACH, which is comparable to the model version that was used for the Coupled Model Intercomparison Project (CMIP5) experiments (Taylor et al., 2012).A validation of global-scale energy and water flux components of MPI-ESM CMIP5 simulations is given in Hagemann et al. (2013) and Brovkin et al. (2013b).The JSBACH model has two independent schemes for calculating surface albedo and canopy absorption.
(8) Soil albedo α s depends on soil color and the soil litter content (Vamborg et al., 2011).A PFT-specific leaf albedo (α le ) as well as the soil background albedo was derived from MODIS observations (Otto et al., 2011).
The sky-view factor ( sky ), which weights between canopy and soil albedo, is calculated assuming a random leaf-angle distribution as sky where 0 ≤ f vegmax ≤ 1 is the fraction of a model grid cell with vegetation for → ∞, which is similar to f veg in Fig. 1.Note that sky is independent of the sun zenith angle in this parameterization, which therefore does not allow us to simulate mutual shadowing in the canopy or its effect on surface albedo, neither the diurnal dependency of surface albedo nor multiple scattering effects.

Canopy absorption
FaPAR is calculated in JSBACH using the 1-D two-stream approximation based on Sellers (1985), similar to JULES.The faPAR is calculated for direct and diffuse radiation components separately using n = 3 canopy layers with equally distributed leaf densities.Further assumptions are that the leaf reflectivity and transmissivity are equal (r = t = const).The canopy single-scattering albedo is assumed to be ω = const = 0.12.Note that this is independent of the leaf albedo used for the albedo calculations.A spherical leaf angle distribution (G(θ ) = 0.5) is also assumed in JSBACH, which is similar to JULES.Multiple scattering effects are taken into account, similar to JULES.
To account for sparse vegetation, faPAR is corrected for the vegetation fraction f veg .This vegetated fraction is parametrized for each plant functional type as where max and γ are PFT-specific parameters representing the maximum leaf area index for a particular PFT and an empirical clumping parameter, respectively.The calculation of f veg ( max ) for faPAR calculations is therefore slightly different than for the calculation of the surface albedo but also independent of the solar zenith angle.While the assumptions for the canopy radiative transfer simulations and technical implementations are different between JSBACH and JULES, it can be shown that both models www.biogeosciences.net/11/1873/2014/Biogeosciences, 11, 1873-1897, 2014 provide nearly identical results for faPAR when parametrized in the same way.
The JSBACH RT scheme allows to account explicitly for gaps within the canopy.It is assumed that, due to canopy gaps, only a fraction f clump < 1 of the canopy area is actually covered with leaves.The radiative transfer simulations are then calculated with an effective LAI ( eff = /f clump ).Note that the leaf area index used for the canopy radiative transfer simulations is therefore larger than the original leaf area index ( eff > ).The simulated radiation fluxes are therefore rescaled by multiplying by f clump after the 1-D RT calculations have been performed.To assess the impact of this clumping on the simulations in the present study, two different JSBACH model versions (with/without clumping) are evaluated.

ORCHIDEE
ORCHIDEE (ORganizing Carbon and Hydrology In Dynamic EcosystEms) is a land surface model that simulates the energy and water fluxes of soil and vegetation, the terrestrial carbon cycle, and the vegetation composition and distribution (Krinner, 2005).
ORCHIDEE is used as the land surface scheme of the IPSL Earth system model.A global model validation is described in Dufresne et al. (2013) for the CMIP5 simulation experiments.We use ORCHIDEE version 1.9.5.2 (revision 816).No changes have been introduced to the calculation of albedo and faPAR since the version that was used for the CMIP5 simulations (Dufresne et al., 2013).The calculation of the light absorbed by vegetation and vegetation albedo is calculated in two independent schemes.

Albedo
The albedo calculation in ORCHIDEE is the same as in JS-BACH (Eq.8).Only the PFT-specific canopy albedo differs slightly compared to JSBACH.The soil background albedo has been derived from a database of soil colors in the standard setup of ORCHIDEE according to Wilson and Henderson-Sellers (1985).

Canopy absorption
The faPAR is not a standard output variable in ORCHIDEE.The light absorbed by the canopy is calculated by means of Beer's law assuming a constant extinction coefficient (G(θ ) = 0.5) for all PFTs and assuming an exponential profile of leaf area index ( * ) within the vertical canopy profile.FaPAR is calculated once and stored in a look-up table.The calculation is independent of the solar zenith angle The vertical distribution of LAI is calculated using n = 20 different canopy layers.The LAI for a given layer i is calculated as * (i) = max e 0.15(i−1) − 1 e 0.15n − 1 , whereas max = 12 [m 2 m −2 ], independent of the PFT used.The parameter max is empirical and the value of 12 (m 2 m −2 ) gives the value of LAI for the last index of the tabulated LAI.This value has been chosen as a maximum value that is unrealistic to occur in reality.Thus all possible LAI values can be covered.The actual LAI, as calculated dynamically by a phenology model, is compared to this precomputed profile to determine the number of canopy layers, which is thus not fixed in OR-CHIDEE.The LAI for this layer is then used to compute fa-PAR using Eq. ( 11).A maximum number of 20 canopy layers is reached at the theoretical LAI of * (20) = 12 (m 2 m −2 ).

RAMI4PILPS experiments
The RAMI4PILPS suite of virtual experiments was designed to evaluate the accuracy and consistency of shortwave RT formulations (as used in DGVMs) under perfectly controlled experimental conditions.More specifically, RAMI4PILPS prescribed a series of virtual canopy scenarios having accurately described structural, spectral and illumination-related characteristics.For these test cases, RT model simulations have been generated as a reference using a Monte Carlo approach.The Monte Carlo model in question had been extensively verified during previous RAMI phases, e.g., Widlowski et al. (2007).Models participating in the RAMI4PILPS benchmarking exercise (Widlowski et al., 2011) had to simulate the canopy albedo, transmission and absorption in both the visible and near-infrared spectral domains.The resulting data were evaluated against the 3-D Monte Carlo reference solution.
Contrary to efforts comparing model simulations against in situ observations at specific test sites, the RAMI4PILPS approach eliminates uncertainties arising from incomplete or erroneous knowledge of (1) the structural, spectral, and illumination-related characteristics of the canopy target, and (2) the uncertainties introduced into the reference solution by calibration, sampling and upscaling errors (Fig. 2).
The complexity of the RAMI4PILPS scenes was adapted to the typical capability of available shortwave RT formulations.More specifically, RAMI4PILPS proposed two homogeneous plant canopy types, i.e., grasslands (GRA) and closed forest canopy (CFC) scenes that differ only in their predominant leaf orientations as well as in the height of the canopy.This corresponds to the idealized case of a 1-D turbid medium, where 1-D canopy radiative transfer schemes are expected to perform best.
In addition, RAMI4PILPS proposed two heterogeneous canopy scenarios where tree crowns were approximated by woodless spheres, i.e., a shrubland (SHR) and an open forest canopy (OFC) scene that differed only in their size of the spheres, the height of the canopy, and the degree of mutual shading between neighboring crowns (Fig. 2).Details of the RAMI4PILPS experiments used in the present study are summarized in Table 1, while further details of the RAMI4PILPS experiments can be found in Widlowski et al. (2011).For each scenario, simulations for different leaf areas and varying soil brightness are performed, assuming direct insulation for three different sun zenith angles as well as isotropic illumination conditions.

Interfacing with RAMI4PILPS
To enable all three models to simulate the RAMI4PILPS scenes, the experiment parameters were implemented as follows: -Single scattering albedo: JULES uses the vegetation single scattering albedo (ω) as an input.It was calculated from the RAMI4PILPS leaf reflectance and transmittance as ω = r + t = 0.1301.JSBACH assumes ω to be constant for calculation of canopy absorption.The standard value used in JSBACH simulations (ω = 0.12) was also used in the present study.The impact of this assumption was however tested by comparing JSBACH simulations with ω = 0.12 and ω = 0.1301.A minor impact on the results of this study was found, indicating the minor importance of ω for the present study.The leaf albedo (α le ) and soil albedo (α s ) were both obtained from Table 1.OR-CHIDEE does not use any information about the single scattering albedo.
-Snow-covered areas: One of the RAMI4PILPS simulations assumes snow below the canopy (SNW).In these cases, the soil albedo was replaced by the snow albedo.No snow on the canopy was simulated as this was not foreseen in the RAMI4PILPS experiments.
-Open vs. closed canopies: Open canopies consist of vegetated and non-vegetated (soil) patches, as discussed before.The RAMI4PILPS experiment descrip-tion contains information on the fractional area covered by vegetation (f veg ) as well as the amount of uncollided radiation (P gap (θ = 0)).
All models rely on the assumption of a plane-parallel turbid media, which prohibits direct application to open canopies.Open canopies are challenging to simulate for 1-D RT models and require effective RT model parameters, as discussed in Sect.2.1.The applied models try to mimic open canopies by correcting by the area fraction covered by vegetation, as discussed before.
The vegetation fraction f veg was therefore obtained from the RAMI4PILPS experiment setup (Table 1).
The soil albedo α s was used for the understory.For JS-BACH we set f vegmax = f veg .The surface albedo for open canopies was calculated for JSBACH and OR-CHIDEE by using Eq. ( 8), while for JULES the surface albedo was calculated as where α le is calculated by the JULES RT model.The total faPAR was calculated in a similar way by weighting the faPAR calculated by the RT models (faPAR canopy ) by the actual area fraction covered by vegetation as The spectral definition in all models and in RAMI4PILPS is consistent.The solar spectrum was divided into two broad spectral bands (VIS: 400 to 700 nm, NIR: 700 to 3000 nm).Both bands were analyzed for albedo in this study, while only the VIS band was used for the faPAR analysis.

Impact assessment -does it matter?
This paper aims at identifying systematic deviations for albedo and faPAR for well-established DGVMs and to provide guidance for estimating the impact of these deviations A comprehensive impact assessment would, however, require detailed knowledge of the spatial distribution of fa-PAR and albedo deviations and their consideration in coupled models, with all else being the same including LAI and all other vegetation variables.This is far beyond the scope of the present paper.It is however of major importance for further model development to raise awareness of potential impacts of the identified deviations.We will detail in the following how we performed straightforward calculations for a first assessment of potential impacts caused by deviations in RT schemes.These are likely to affect net photosynthesis first by differences in absorbed radiation through the canopy at a given LAI and second through heat exchange between the land surface and the atmosphere through albedo.Whereas DGVM behavior with respect to energy, water and carbon exchange will be affected in multiple other ways, we constrained our assessment to these two prime targets.This impact assessment is expected to provide an initial estimate of direct and potential impacts for surface energy and carbon fluxes.More complex analyses using coupled climate models would however be needed for a proper impact assessment on global spatial and century time scales.

Potential impacts of RT scheme differences on the C cycle
In order to evaluate the potential implications of differences in canopy faPAR profiles, induced by different canopy RT schemes, we determined net photosynthesis rates for each of the virtual experiments of the closed forest canopy case only (CFC).We made this constraint because the photosynthesis calculations in combination with the representation of faPAR for different layers of the canopy are likely only valid for the closed forest canopy case, given the assumptions made in the photosynthesis schemes adopted in most DGVMs.To determine net photosynthesis estimates, we chose the Farquhar scheme as described by von Caemmerer (2000), which is representative of the photosynthesis schemes in current DGVMs.According to this scheme, net photosynthesis is the minimum of CO 2 -limited photosynthesis and photosynthesis as limited by the availability of ribulose bi-phosphate (RuBP), the latter being driven by absorbed radiation, corrected for dark respiration.
We evaluated this scheme for two representative closed canopy situations under typical clear sky conditions in summer, following Widlowski et al. (2011).The two selected situations coincide with a tropical forest with a solar zenith angle at local noon ( noon ) of 15 • and a boreal (presumably evergreen needle-leaved) forest with noon = 50 • .For these situations, and the given sun zenith angles of the virtual experiment, the total incident shortwave radiation at noon was calculated according to Wang et al. (2002).When combined with the faPAR profiles through the canopy, this provided the absorbed photosynthetic active radiation (aPAR) in each canopy layer.Coupling aPAR to the Farquhar scheme provided a CO 2 -limited and light-limited photosynthesis rate in each canopy layer.The minimum of these represents the net photosynthesis rate, which was summed across all layers to provide the net photosynthesis rate (µmol C m −2 s −1 ) for a given virtual experiment.Absolute differences between net photosynthesis rates between the RAMI4PILPS scenes and each of the DGVMs were compared.
The Farquhar scheme was parameterized according to Sharkey et al. (2007), which currently provides the most comprehensive analysis of the variables involved.Estimates of the maximum rates of Rubisco carboxylase activity (V c max ) and the maximum rates of photosynthetic electron transport (J max ), were however taken to represent the mean observed values at a reference temperature of 25 • C for tropical broadleaved evergreen forests and boreal (needle-leaved) evergreen forests, according to the TRY database ( (Kattge et al., 2011) and Domingues et al. (2010)), as compiled by Verheijen et al. (2013).Dark respiration was scaled to V c max .The used calculation scheme neglects effects of water, thermal or nutrient limitations.Maximum canopy productivity rates (J max , V c max ) are, however, based on observed mean values that may implicitly take some of these limitations into account.

Surface radiation budget
Surface albedo directly affects the net surface radiation budget and thus indirectly also near surface temperature.To quantify the influence of a factor that changes the balance of incoming and outgoing energy in an Earth-atmosphere system, the concept of radiative forcing (RF) has been developed (IPCC, 2007).Radiative forcing has been widely used to quantify the global impact of regional changes in surface albedo (Hall and Qu, 2006;Pongratz et al., 2009;Lenton and Vaughan, 2009;Bright et al., 2012) and can be calculated in various ways.Here, we applied a simple method to estimate RF (Bright et al., 2012;Cherubini et al., 2012).

Effect of albedo biases
The radiative forcing calculated in this study is understood as a disturbance of radiation fluxes caused by variations in surface albedo.The estimates calculated in this study must therefore not be compared directly with estimates of radiative forcing from different forcing agents like those used, e.g., in coupled climate model assessments (IPCC, 2007).The albedo-induced radiative forcing RF α p is related to a change in planetary albedo α p as where R ↓ TOA (W m −2 ) is the incoming solar radiation flux at the top of the atmosphere (TOA).The change in planetary albedo is linearly related to a change in surface albedo (Lenton and Vaughan, 2009).The corresponding radiative forcing is then given by where a is the two-way atmospheric transmissivity for solar radiation that is given as the product of the downward ( ↓ ) and upward ( ↑ ) atmospheric transmissivities.↑ was parameterized in two ways.Similar to previous studies focusing on radiative forcing (Lenton and Vaughan, 2009;Cherubini et al., 2012), we assumed a constant value of ↑ = 0.854.This global mean value is based on the assumption that the outgoing shortwave flux leaves the atmosphere given clear sky conditions.Alternatively, we assumed that under cloudy sky conditions, the upward and downward atmospheric transmissivities are equal ( ↑ = ↓ ).The TOA solar radiation flux can be calculated for a particular point on Earth using simple formulations (Widlowski et al., 2011;Bright et al., 2012).The atmospheric one-way downward transmissivity ↓ a can be derived from shortwave surface all-sky (R ↓ surf ) and TOA flux (R ↓ TOA ) estimates as available from satellite data products (Rossow and Zhang, 1995;Posselt et al., 2012;Loeb et al., 2012) as Let us now assume that we observe a surface albedo bias ( α s ) for a particular plant functional type (PFT) and that this PFT covers a fraction f pft of a model grid box; then the total radiative forcing for the Earth with surface area A = 510 × 10 6 (km 2 ) is given using Eq. ( 16) by integrating over the surface area as where "a" is a spatial index.If α s is assumed to be constant (time/space independent) for a particular PFT, then the sensitivity to surface albedo changes can be defined as which can be calculated for each combination of PFT fraction f pft and Figure A4 shows the mean fields of R ↓ surf and the estimated one-way transmissivity over land as derived from CERES EBAF (v2.6) (Loeb et al., 2012).Using Eq. ( 19) and values in Table 2 one can estimate the radiative forcing caused by an albedo change for a given distribution of plant functional types.

Impacts of neglecting the diurnal cycle on radiative forcing
The surface albedo can change substantially throughout the day due to its dependency on the solar zenith angle.The JS-BACH and ORCHIDEE surface albedo models are representative of a wide range of albedo models applied in DGVMs.These models do not simulate the surface albedo diurnal behavior, which may affect the surface radiation budget.
The surface albedo can be expressed as a weighted average between the so called black-sky albedo (BSA) and the albedo for complete isotropic diffuse radiation (white sky albedo, WSA) (Lewis and Barnsley, 1994;Schaepman-Strub et al., 2006) as where f dir (-) corresponds to the fraction of the direct surface solar radiation flux.
The surface albedo models of JSBACH and ORCHIDEE are based on the white sky albedo (α WSA ).potential impacts of neglecting the diurnal cycle in surface albedo calculations, we used the difference between the albedo given by Eq. ( 20) and α WSA , which is given as Assuming clear sky conditions, we calculated the surface solar radiation flux R ↓ surf for each latitude φ and day of the year using the MAGIC atmospheric RT code (Mueller et al., 2012;Posselt et al., 2012).Required atmospheric aerosol and water vapor content were taken from climatological mean values (Kinne et al., 2013).
The temporal average RF disturbance due to neglecting the diurnal cycle of α, averaged over a time period T , is then given as We calculated Eq. ( 22) by using either ↑ a = const = 0.854 (Lenton and Vaughan, 2009) or by assuming that ↑ (t, a) = ↓ (t, a), where the one-way downward transmissivity was obtained from Eq. ( 17) assuming clear sky conditions and as calculated by MAGIC.The fraction of direct surface radiation flux (f dir ) was calculated from monthly climatologies of direct and total radiation fluxes as obtained from satellite radiation products and provided by the NASA Atmospheric Data Center (https://eosweb.larc.nasa.gov/).The global radiative forcing RF tot was calculated for each day (T = 24 h) and two different PFT distributions (tree, grass) as derived from MODIS vegetation.Continuous fields were used (DiMiceli et al., 2011) that provide f pft for each grid cell.

Consistency of canopy absorption and reflection calculations
The canopy radiative transfer needs to fulfill the law of conservation of energy.Following Eq. ( 3) this is defined as where 0 ≤ can ≤ 1 is the one-way canopy transmission.The canopy RT schemes in the used DGVMs differ in how they simulate canopy absorption and reflection as discussed before.Two models (JSBACH, ORCHIDEE) use separate approaches for the simulation of canopy absorption and reflection, which might lead to a lack of energy conservation in these models.
To test the energy conservation, the canopy one-way transmissivity would be needed as model output.This is however not available from all of the investigated models.We therefore adopted a simpler approach to test in general for energy conservation of the used canopy RT schemes.From Eq. ( 23) it follows that can = To test if the models are in general energy conserving, we investigated if can from Eq. ( 24) is within its physically defined limits (0 ≤ can ≤ 1), which is the most conservative approach.

Canopy absorbed radiation
The investigated models generally underestimate the radiation absorbed in the canopy for the RAMI4PILPS experiments analyzed (Fig. 3).Especially for the open canopy cases (SHR, OFC), a strong underestimation of faPAR is observed with increasing negative bias for an increasing solar zenith angle.The 1-D canopy RT models without clumping (JULES, JSBACH no-clumping) perform best for the closed canopy cases (CFC, GRA), as expected.Some positive biases are observed for the isotropic (ISO) illumination conditions for JULES and JSBACH noclumping.In general, results of these models are almost identical due to the same basis for the canopy RT simulation code as discussed in Sect.2.3.

Does it matter?
Deviations in net photosynthesis due to differences in the fa-PAR canopy profile for the RAMI4PILPS reference scene vs. each of the models strongly depend on the virtual experiment (Fig. 4).On average, deviations are largest for θ = 60 • , where deviations in faPAR are considerable, while incident radiation is still high enough to support net photosynthesis (in contrast to θ = 83 • ).Moreover, the differences in net photosynthesis strongly increase with increasing LAI, when more leaf surface is available for photosynthesis.Finally, the impacts are also larger for the boreal forest ( noon = 50) and slightly increase with increasing albedo, although the latter effects are minor.
Depending on the model and even within a given model setting, differences in calculated net photosynthesis can be both negative and positive.The largest deviations occur for ORCHIDEE, with both overestimations as well as underestimations of net photosynthesis rates, although overestimations by ORCHIDEE predominated.Net photosynthesis estimates according to JSBACH are on average lower than those for the RAMI4PILPS reference profiles.The deviations in net photosynthesis for JSBACH no-clumping and JULES are minor overall, indicating a good performance of the canopy 1-D RT schemes in these cases.This is expected, because the closed canopy corresponds to the idealized case where the 1-D RT schemes are supposed to perform best.Overall, JULES performed the best of all the models, although even for this model substantial deviations occur.
These deviations in net photosynthesis rates can be considerable, up to 10 µmol C m −2 s −1 , corresponding to up to 25 % of the photosynthesis rates.These conditions (with high LAI and θ = 60 • ) commonly prevail.Therefore, the impacts for the total calculated carbon budget are likely to be similarly large, which implies that the deviations in faPAR profiles detected in the current study are of critical importance for global carbon cycle studies.At lower incident radiation conditions, e.g., for θ = 83 • , deviations even amount to more than 75 %.However, given the lower incident radiation and consequently lower net photosynthesis, their impacts on the global carbon balance will be less.

Surface albedo
Figure 5 shows deviations of simulated surface albedo compared to RAMI4PILPS reference solutions for the entire solar spectrum (0.3 to 3.0 µm).Deviations in the entire solar spectrum are of particular importance as these determine the net effect on the surface radiation budget.More details on the separate deviations in the visible and infrared spectral bands are provided in the Appendix (Figs.A1 and A2).
The closed canopy experiments (CFC, GRA) show very similar results, as expected given that they only differ in their leaf orientation (random vs. erectophile).The largest deviations from the RAMI4PILPS references solutions are found    for bright surface background (SNW) for the JSBACH and ORCHIDEE models.Deviations for these two models are between 0.18 ≤ α ≤ 0.25 for SNW experiments.In contrast, the JULES albedo simulations show much smaller deviations for all closed canopy experiments (−0.03 ≤ α ≤ 0.02).Also here, largest deviations are observed for bright backgrounds.
For soils with medium albedo (MED), the deviations are smaller for JSBACH and ORCHIDEE.Positive biases are obtained for isotropic conditions and small solar zenith angles (−0.01 ≤ α ≤ 0.06).For large solar zenith angles, JULES shows minor deviations, while JSBACH and OR-CHIDEE show a negative bias of α = −0.06.While the deviations for JSBACH and ORCHIDEE become smaller with an increasing leaf area index for bright surfaces, due to an increased masking of the bright background, α increases with increasing LAI for the medium bright soils.
Open forest canopies (OFC) and shrubland (SHR) simulations show larger deviations than the closed canopy cases.Large positive deviations (0.07 ≤ α ≤ 0.36) are observed for the snow-covered cases, with increasing deviations with increasing solar zenith angle.The larger deviations result from a lack of the representation of canopy shadowing effects in all of the canopy RT models, resulting in an overestimation of the simulated surface albedo.For medium soils, a positive bias is observed for all models.The deviations are larger (0.01 ≤ α ≤ 0.06) than for the closed canopies and an increase in α with increasing leaf area index is observed.

Effect of systematic albedo biases
It was shown in the previous section that the investigated models show albedo biases of −0.27 < α < 0.36.The extreme deviations however occur typically either at very large sun zenith angles or at open forest canopies with snow as a background.In the case of snow-free and closed canopies, the typical albedo biases for isotropic illumination conditions and small solar zenith angles are on the order of 0.02 < α < 0.05 for closed canopies.
The impact of albedo biases on the calculated radiative forcing sensitivity (λ) is summarized in Table 3.The mean radiative forcing sensitivity for trees (grassland) is −5.5 (−13.11)(W m −2 ).Thus a change in surface albedo of, e.g., α = 0.04 would correspond to a radiative forcing of −0.22 (−0.52) (W m −2 ).The use of different global radiation data sets as input has a minor impact on the radiative forcing estimates.The standard deviation of RF caused by different radiation data is 0.3 ≤ σ ≤ 0.45 (1.03 ≤ σ ≤ 1.41) for trees (grassland).However, much larger differences are observed for different assumptions for ↑ .For clear sky conditions ( ↑ = const), the radiative cooling effect is 50 to 60 % larger  than for cloudy skies ( ↑ = ↓ ).It is therefore expected that the radiative forcing for the two investigated PFT types will be somewhere in between these two extremes.

Biogeosciences
Figure 6 shows the spatial distribution of the temporal mean radiative forcing sensitivity for tree and grassland for both assumptions of ↑ .The much larger RF sensitivity for in open canopy cases with high solar zenith angles and/or snow as a background.Only in a few instances higher fa-PAR than in the RAMI4PILPS reference were found (typically in closed canopy cases with low zenith angles or isotropic insulation).Moreover, deviations were strongly 965 model-dependent.JSBACH and ORCHIDEE consistently had the highest (negative) biases.This suggests that a twostream approximation such as proposed by Sellers (1985) or other consistent treatments of radiative transfer are essential for obtaining unbiased faPAR calculations Pinty et al. (2006).The simplified surface albedo schemes of JSBACH and ORCHIDEE do not take into account the diurnal variabil-985 ity in surface albedo.Isotropic insulation conditions are assumed instead.The impact of this assumption on global scale was therefore analyzed.While the applied assessment of radiative forcing is simple, we use state-of-the art observations to get a most realistic assessment of RFtot and its uncer-990 tainties.The developed framework can be used in general to estimate radiative forcing effects of temporally invariant as well as changing surface albedo biases.
While we have identified faPAR and albedo biases by comparing to RAMI4PILPS experiments, it needs to be em-995 phasized that the RAMI4PILPS experiments provide idealized cases of canopies and that real surfaces and vegetation patches are much more complicated.For instance, none of the analyzed canopy RT schemes is capable to take into account large scale vegetation structural effects, like e.g.mu- grassland is mainly due to its larger spatial coverage globally as well as due to the higher abundance in tropical areas with large insulation.For trees, the highest radiative cooling is observed in the tropics, with a secondary maximum in boreal regions, which is mainly during the boreal summer.
Given the similarity in albedo bias for closed forest and grassland canopies, compared to RAMI4PILPS references, the radiative forcing of the forest and grassland PFTs might be summed.The global mean RF sensitivity is therefore −9.31 (W m −2 ), assuming equal probability for the two parameterizations of atmospheric transmissivity.For a typical albedo bias of α = 0.04, this would result in RF tot = −0.37 (W m −2 ).

Diurnal effects
Figure 7 shows the seasonal cycle of the changes in radiative forcing caused by neglecting the albedo diurnal cycle for grassland and tree-covered areas for different leaf areas and parameterizations of the upward atmospheric transmissivity ( ↑ ).
The maximum radiative cooling effect is −1.25 < RF tot < −0.8 (W m −2 ).The radiative cooling is mainly caused by large albedo biases in tropical areas and changes substantially throughout the year.With increasing leaf-area index, the radiative forcing effect increases, which can be explained by different reflectances in the VIS and NIR domains as a function of LAI.While the albedo decreases in the VIS with increasing leaf area index (darkening), the albedo in the NIR increases, which results in an overall increase in the surface albedo and thus the radiative effect with an increasing LAI.In general, the radiative cooling by grassland areas is higher than for tree-covered areas due to their different spatial distribution.The annual mean radiative forcing is −0.61 (−0.53) and −0.84 (−0.73) for a leaf-area index of one and four, respectively, assuming   tual shadowing or the scattering between multiple canopy elements.

Do we (need to) care?
Existing representations of shortwave radiative flux estimates in DGVMs have been developed throughout the last 1005 decades.Developing a DGVM always requires prioritization of resources between different model components.Recent DGVMs become more and more complex and incorporate complex carbon fluxes and nutrient cycling.In the past, a lot of attention has been devoted towards an improvement 1010 of especially the biochemical components in DGVMs (Goll et al., 2012;Zaehle and Friend, 2010;Brovkin et al., 2012), whereas biophysical components, like the canopy RT have been neglected because no study had analyzed systematically the major caveats of existing schemes so far.While model 1015 benchmarking and evaluation studies detail particular deviations of model results compared to observations and assess a (relative) model skill (Gleckler et al., 2008;Luo et al., 2012;Brovkin et al., 2013b;Hagemann et al., 2013), they typically don't provide an assessment of potential impacts of the ob-1020 served model deficits.
It was therefore one of the objectives of the present study to provide an assessment of the implications due to a choice for a particular canopy radiative transfer scheme.Answering the question "Does it matter?" is of particular importance 1025 when decisions on further DGVM development need to be made with limited resources.
The differences in faPAR discussed in the previous section do not translate straightforwardly into differences in net photosynthesis, and depend on total radiation and the extent 1030 to which CO2 or light is limiting photosynthesis in different canopy layers (which again depends on the maximum capacity of the light and CO2 limited photosynthesis pathways and on faPAR profiles through the canopy).Here, we only evaluated the impacts on net photosynthesis for some of the best 1035 (i.e.least deviating) cases, i.e. for closed forest canopies.We chose this setting, because the assumptions of the implementation of the Farquhar-photosynthesis scheme in DGVMs only apply to closed (forest) canopies.In that sense, it is important to recognize that some impacts on net photosynthesis 1040 Fig. 7. Difference in radiative forcing (cooling) due to neglecting the diurnal cycle for surface albedo schemes for global tree and grass cover and different leaf areas (columns).Different rows correspond to different parameterizations of the shortwave upward oneway atmospheric transmissivity ( ↑ ).Values in parentheses correspond to global mean values.
).While the LAI and seasonal cycle have the largest impacts on RF tot , a different choice in the parameterization of ↑ also results in deviations of ≈ 15 %.

Consistency of surface energy fluxes
Energy conservation was evaluated for the different experiments following Appendix D. Results are illustrated in Fig. A5.The JSBACH and ORCHIDEE models are not energy conservative for high leaf areas and bright backgrounds (SNW).The JSBACH model shows clear deviations for the closed canopy cases (CFC, GRA) for almost all experiments with bright backgrounds.
The JULES model is the only model that simulates the canopy RT fluxes in a consistent manner.No obvious violation of the energy conservation is found for JULES.
It needs to be emphasized however that the energy conservation might also be violated in the cases where the calculated canopy transmission is within its physical limits (0 ≤ can ≤ 1).It can however not be quantified in this study as the canopy transmission information is not available for the different models.JSBACH and ORCHIDEE estimates for can often differ from those of JULES by a factor of 2 or more.If we assume that JULES simulations are closest to energy conservation, it is rather likely that these models would not be energy conserving for experiments where a large deviation from the transmission estimates of JULES is observed.

Discussion
DGVMs are widely used to simulate land surface dynamics in coupled Earth system models, to assess the role of the biosphere in the Earth system and to study global water, energy and carbon fluxes.This paper aims to raise awareness of potential problems and impacts associated with different canopy radiative transfer schemes and their impacts on carbon fluxes and global climate simulations.
Large deviations in faPAR were found.In most cases, fa-PAR was underestimated with up to 60 % underestimation in open canopy cases with high solar zenith angles and/or snow as a background.Only in a few instances were higher faPAR than in the RAMI4PILPS reference found (typically in closed canopy cases with low zenith angles or isotropic insulation).Moreover, deviations were strongly model dependent.JSBACH and ORCHIDEE consistently had the highest (negative) biases.This suggests that a two-stream approximation such as proposed by Sellers (1985) or other consistent treatments of radiative transfer are essential for obtaining unbiased faPAR calculations Pinty et al. (2006).For surface albedo, the largest deviations in surface albedo were observed in the case of bright background albedo (SNW), for both open and closed canopies.High deviations ( α > 0.2) were diagnosed for the simplified surface albedo schemes implemented in JSBACH and ORCHIDEE.The albedo model of JULES only shows minor deviations (−0.03 ≤ α ≤ −0.01) for closed canopy cases, but also large deviations for the open canopy cases.For surfaces with a typical soil background and a closed canopy, the 1-D canopy RT model of JULES is superior compared to the simple albedo schemes of JSBACH and ORCHIDEE, which have an albedo bias that is twice to five times larger than the JULES albedo bias.

Biogeosciences
The simplified surface albedo schemes of JSBACH and ORCHIDEE do not take into account the diurnal variability in surface albedo.Isotropic insulation conditions are assumed instead.The impact of this assumption on a global scale was therefore analyzed.While the applied assessment of radiative forcing is simple, we use state-of-the art observations to get the most realistic assessment of RF tot and its uncertainties.The developed framework can be used in general to estimate radiative forcing effects of temporally invariant as well as changing surface albedo biases.
While we have identified faPAR and albedo biases by comparing them to RAMI4PILPS experiments, it needs to be emphasized that the RAMI4PILPS experiments provide idealized cases of canopies and that real surfaces and vegetation patches are much more complicated.For instance, none of the analyzed canopy RT schemes is capable of taking into account large-scale vegetation structural effects like, e.g., mutual shadowing or the scattering between multiple canopy elements.

Do we (need to) care?
Existing representations of shortwave radiative flux estimates in DGVMs have been developed throughout the last decades.Developing a DGVM always requires prioritization of resources between different model components.Recent DGVMs become more and more complex and incorporate complex carbon fluxes and nutrient cycling.In the past, a lot of attention has been devoted towards an improvement of especially the biochemical components in DGVMs (Goll et al., 2012;Zaehle and Friend, 2010;Brovkin et al., 2012), whereas biophysical components, like the canopy RT, have been neglected because no study had analyzed systematically the major caveats of existing schemes so far.While model benchmarking and evaluation studies detail particular deviations of model results compared to observations and assess (relative) model skill (Gleckler et al., 2008;Luo et al., 2012;Brovkin et al., 2013b;Hagemann et al., 2013), they typically do not provide an assessment of potential impacts of the observed model deficits.
It was therefore one of the objectives of the present study to provide an assessment of the implications due to a choice of a particular canopy radiative transfer scheme.Answering the question "Does it matter?" is of particular importance when decisions on further DGVM development need to be made with limited resources.
The differences in faPAR discussed in the previous section do not translate straightforwardly into differences in net photosynthesis, and depend on total radiation and the extent to which CO 2 or light limits photosynthesis in different canopy layers (which again depends on the maximum capacity of the light and CO 2 -limited photosynthesis pathways, and on fa-PAR profiles through the canopy).Here, we only assessed the impacts on net photosynthesis for some of the best (i.e., least deviating) cases, i.e., for closed forest canopies.We chose this setting because the assumptions about the implementation of the Farquhar photosynthesis scheme in DGVMs only apply to closed (forest) canopies.In that sense, it is important to recognize that some impacts on net photosynthesis (discussed below) may be even worse for open shrub or forest canopies, let alone for situations where one PFT literally grows on top of another, affecting the light regime of the one below.
Moreover, given that the differences in faPAR are situation and model specific, the impacts are also expected to differ for various regions and models.Thus, we selected two commonly closed canopy cases (boreal forests and tropical forests), which we parameterized with representative radiation at the top of canopy and maximum photosynthesis capacities.In contrast to faPAR, net photosynthesis was most affected at medium zenith angles (when radiation at the top of the canopy was still high) and high LAI (while faPAR was almost unaffected by LAI).While in this study net photosynthesis rates and their biases were not integrated to diurnal values, conditions of medium zenith angle occur commonly and large deviations in net photosynthesis estimates are thus expected.These deviations were on average stronger for the boreal forest than for the tropical forest case.
Likely, these deviations have been partly captured by tuning maximum photosynthesis rates (which may explain their deviations from the observed means (Kattge et al., 2011;Verheijen et al., 2013)).Even so, given that the magnitude and direction differ between different model experiments, tuning will not allow capturing of regional biases.While other factors like biases in precipitation also likely contribute to regional differences in GPP estimates between DGVMs and flux observation-derived estimates (Beer et al., 2010), our results suggests the importance of incorporating consistent and adequate radiative transfer schemes.
Likewise a typical surface albedo bias on the order of α = 0.04 results in a mean radiative forcing of RF tot = −0.22(−0.52)(W m −2 ) for trees (grasslands).In addition it was shown in Sect.coupled (land-atmosphere) climate model, the provided values are likely to provide a realistic range of RF tot estimates.These radiative forcing effects correspond to considerable cooling, which is mainly caused by tropical land areas.This emphasizes that the diurnal cycle cannot be neglected for albedo calculations.Especially in tropical areas, where the insulation conditions change strongly throughout the day and solar insulation is high, a large effect on RF tot is observed.It needs to be emphasized that these results for net photosynthesis and radiative forcing are only valid for some vegetation settings and are tailored to provide only the order of magnitude of the impact of faPAR and surface albedo biases on net photosynthesis and radiative forcing.They should therefore be considered as a first step towards a more thorough analysis of the impact of biases in canopy RT schemes in DGVMs.To quantify the impact on global and regional scales as well as for longer timescales it would be required to couple different canopy RT schemes in an Earth system model and perform dedicated simulations to assess differences in global energy, water and carbon fluxes.While this is beyond the scope of the present study, we outline some ways forward that could lead to an improvement in DGVM canopy radiation flux estimates.

A way forward?
A combination of different approaches could lead to a general improvement in the representation of canopy radiation fluxes in DGVMs.
-Model evaluation is essential for any model development.Only once the target is clearly defined and accurate reference solutions are available, can the potential biases in canopy RT schemes be rigorously identified.In the case that tolerance criteria are specified that allow us to judge the relevance of deviations between models and a defined reference, model evaluation studies can determine the urgency of new developments to improve relevant model formulations (Widlowski et al., 2013).Specific test cases like those provided by RAMI4PILPS are required for that purpose as well as large-scale model evaluation studies (Hagemann et al., 2013;Brovkin et al., 2013b).The current RAMI4PILPS experiments provide a very idealized setup that does not represent real canopies in a realistic manner.To approach a fully quantitative and realistic assessment of canopy RT schemes in climate models, several aspects are of particular importance: (a) snow cover strongly affects the surface albedo.Current RAMI4PILPS experiments assume snow only below the canopy, while the surface albedo of a model grid cell is largely determined by snow on top and within the canopy.It is therefore recommended to represent snow cover dynamics in a more realistic way in RAMI4PILPS.(b) To allow for an appropriate assessment of the impact of biases in surface albedo and faPAR, it is important to provide reference solutions that cover a large variety of illumination conditions (representation of diurnal cycle) and that provide information on the reflected and absorbed fluxes, also within distinct canopy layers.This is required, e.g., for the quantification of the effect of the faPAR bias on net photosynthesis as discussed.(c) RAMI4PILPS neglects the woody part of the vegetation (stems, branches).These have a major effect on the surface albedo at high latitudes, especially when the leaf area index is negligible and the solar zenith angle is large and the ground is covered by snow.In these cases, stems contribute a large part to shadowing of the bright background (snow).An accurate representation of this masking effect in models is important for the sensitivity of surface temperature to surface albedo.(d) Vegetation types within RAMI4PILPS that more closely resemble plant functional types as defined in DGVMs would allow easier translation of results.
-Model re-parameterization could be used to compensate (partly) for some of the biases observed between models and references.Structural deficits in canopy RT schemes are not likely to be overcome with a reparameterization of models, but could reduce their effect on model simulations.Some of the observed faPAR and albedo biases may be reduced by simply changing model parameters like, e.g., leaf albedo (α leaf ), canopy single scattering albedo (ω) or specific leaf area (SLA) to artificially converge model simulations towards mean reference data.This might reduce considerably potential biases and thus potential impacts on, e.g., global radiative forcing, although it will not allow alleviating of regional biases.
-Improved consistency of existing schemes: the present study has clearly shown that simplified albedo schemes result in considerable biases and radiative forcing which is not negligible.As a first step, the simplified surface albedo models should be replaced by 1-D RT schemes that are consistent with those used by the same models when computing faPAR.An adaptation of the 1-D RT model for a consistent calculation of canopy absorption and reflectance should require medium effort, would be physically consistent and could also provide more accurate calculations of the surface albedo diurnal cycle.
-Effective radiative model parameters: a replacement of 1-D canopy RT schemes by more complex 3-D canopy RT schemes in DGVMs is not foreseeable in the future as the required canopy structural information required for 3-D models is not available as prognostic variables from DGVMs and because of the much higher computational costs of these models.Effective radiative model parameters, however, aim to  (Pinty et al., 2006).Effective radiative model parameters therefore compensate for inherent differences in process descriptions between simple (1-D) and complex (3-D) radiative modeling approaches.
The key challenge, however, is to develop a deterministic or empirical relationship between model state variables and radiative model effective parameters that satisfies the consistency of canopy RT fluxes in agreement with observations.Pinty et al. (2006) have shown that the 3-D canopy radiative transfer problem can be approximated by a 1-D model using an effective leaf area index as well as effective parameters for the vegetation optical properties.Different approaches have been proposed to relate a leaf area index to an effective leaf area index by defining a structural parameter that accounts for vegetation clumping (Yang et al., 2001;Pinty, 2004;Ni-Meister et al., 2010;Haverd et al., 2012).A new globally available data set of effective RT model parameters as derived from satellite observations (Pinty et al., 2011) could facilitate RT model parametrizations of DGVMs for different plant functional types.Given the fact that it has been shown that 1-D models can provide appropriate RT fluxes when parametrized in the right way, these novel data set provides a novel source for the re-parameterization of canopy RT models at a global scale.
-As an alternative, effective vegetation model parameters or structural parameters may be obtained through re-calibration using global satellite observations.Global satellite observations provide long-term and high-resolution estimates of surface variables, like surface albedo, faPAR and canopy radiation flux partitioning (Tucker et al., 2005;Loew and Govaerts, 2010;Pinty et al., 2011).However, a direct assimilation of geophysical EO products might be complicated due to, e.g., different definitions of geophysical variables among different products as well as between a DGVM-and an EO-based observational data set.Unless the canopy RT schemes in DGVMs do not satisfy basic principles of energy conservation and the appropriate representation of canopy scattering processes, an assimilation of EO-based surface albedo or faPAR products in DGVMs might lead to erroneous results.As an example, faPAR satellite products show a wide range of variability due to different definitions and algorithms applied for the retrieval.A comprehensive review is given by Gobron and Verstraete (2009).Dahlke et al. (2012) compared different globally available faPAR satellite products and found considerable biases between the different data products, whereas the vegetation seasonality was consistently captured by the different investigated data products.Similarly, an assimilation of satellite leaf area indices or faPAR products in DGVMs might be complicated or impossible due to different canopy radiation schemes applied.While DGVMs use 1-D canopy radiation schemes, as discussed, some satellite-based products of vegetation variables (LAI, faPAR) are based on the inversion of 3-D canopy radiation models (Myneni et al., 2002).
As a consequence, it is not physically valid to ingest the retrieved values of for instance LAI into 1-D representations such as those provided by two-stream approaches.The usage of effective variables in data assimilation schemes is required instead.

Conclusions
Major conclusions of the present study are: 1. Considerable biases: currently used simplified canopy RT schemes produce considerable biases in both surface albedo and faPAR calculations.Very simple albedo schemes result in large albedo biases ( α > 0.3) in snow-covered cases and also much higher albedo biases for surfaces with a typical soil albedo.
Our first assessment of the impacts of these biases indicates that these deviations lead to a considerable radiative forcing when not corrected for by an appropriate model parameterization, while the faPAR differences can lead to biased estimates of GPP.
2. Albedo diurnal cycle: considering the albedo diurnal cycle in global surface albedo schemes is of major importance.State-of-the-art DGVM and climate models should therefore implement albedo schemes that consider appropriately the surface albedo diurnal behavior.
3. Physical inconsistency: canopy radiation schemes that use separate approaches to simulate reflected and absorbed radiation fluxes can lead to physically meaningless results.For some of the simulated experiments, it was clearly shown that the investigated models violate the conservation law of energy for the canopy RT fluxes.
4. RAMI4PILPS continued: the current RAMI4PILPS experiments (Widlowski et al., 2011) should be extended to provide more realistic settings for evaluating DGVMs.Such continued RAMI4PILPS experiments may provide the basis for a more comprehensive analysis to assess the impact of actual surface albedo changes on the climate using simulations of a coupled Earth system model.

Fig. 1 .
Fig.1.Illustration of the difference between vegetation fraction (f veg ) and gap probability (P gap (θ = 0)) for a model grid cell that is assumed to be totally vegetated.
, 11, 1873-1897, 2014 www.biogeosciences.net/11/1873/2014/ global mean estimates of the one-way atmospheric transmissivity ( ↓ ) as well as global means of R ↓ TOA and R ↓ surf based on various data sources over land areas.The one-way atmospheric transmissivity ranges from 0.41 ≤ ↓ ≤ 0.58, dependent on the data set.
RT models and RAMI4PILPS for different experiments (model-RAMI4PILPS); columns corref area or vegetation fraction and soil brightness (see

Fig. 4 .
Fig.4.Deviations in GPP for closed forest canopies as derived from differences in canopy radiation profiles between the RAMI4PILPS and canopy RT schemes.Columns correspond to different leaf areas and soil brightness, while rows correspond to different sun zenith angles for tropical ( noon = 15) and boreal ( noon = 50) conditions.Triangles correspond to different models.

Fig. 5 .
Fig. 5. Deviations of simulated surface albedo from RAMI4PILPS reference solutions (α model − α ref ) for the solar spectrum (0.3 to 4.0 µm).Details of the figure are the same as for Fig. 3.

Fig. 6 .
Fig. 6.Impact of a systematic albedo difference of ∆α = 0.04 on radiative forcing for tree cover (top) and grassland (bottom), using CERES solar radiation flux (R ↓ surf ) and different upward atmospheric one-way atmospheric transmissivity: Γ ↑ = const (left), Γ ↑ = Γ ↓ (right) 4 Discussion DGVMs are widely used to simulate land surface dynamics in coupled Earth System models, to assess the role of the biosphere in the Earth System and to study global water, energy and carbon fluxes.This paper aims to raise awareness 955

970
For surface albedo, largest deviations in surface albedo were observed in case of bright background albedo (SNW), for both open and closed canopies.High deviations (∆α > 0.2) were diagnosed for the simplified surface albedo schemes implemented in JSBACH and ORCHIDEE.The 975 albedo model of JULES only shows minor deviations (−0.03 ≤ ∆α ≤ −0.01) for closed canopy cases, but also large deviations for the open canopy cases.For surfaces with a typical soil background and a closed canopy, the 1D canopy RT model of JULES is superior compared to the simple 980 albedo schemes of JSBACH and ORCHIDEE which have an albedo bias which is twice to five times larger than the JULES albedo bias.

Fig. 7 .
Fig. 7. Difference in radiative forcing (cooling) due to neglecting the diurnal cycle for surface albedo schemes for global tree and grass cover and different leaf area (columns).Different rows correspond to different parameterizations of the shortwave upward one-way atmospheric transmissivity (Γ ↑ ).Values in parentheses correspond to global mean values.
3.2.1 that neglecting the diurnal cycle results in a radiative cooling of between −1.25 and −0.8 (W m −2 ), dependent on the leaf area index and atmospheric opacity model chosen.While an actual assessment of the radiative forcing effect would require simulations with a

Figure 5
Figure5shows surface albedo deviations for the entire solar spectrum.The following figures show the deviations of surface albedo for the visible and near infrared shortwave bands separately.Alexander Loew: DGVM canopy radiation

Fig. A1 .Fig
Fig. A1.Deviations of simulated surface albedo from RAMI4PILPS reference solutions (α model − α ref ) for the VIS spectrum (0.3 to 0.7 µm).Details of figure structure are the same as Fig. 3.

Table 2 .
Summary of atmospheric one-way or two-way transmissivity as derived from different sources.Input fluxes used for the transmissivity estimates are provided.Values are for land areas only.

Table 3 .
Radiative forcing sensitivity (∂RF/∂ α) for different radiation data sets as derived from Eq. (19) for different global radiation data sets used to estimate the one-way atmospheric transmissivity and different assumptions about the calculation of ↑ .

Surface albedo deviations for VIS and NIR
Implementing physically consistent canopy RT schemes in DGVMs would allow for the assimilation of satellite-based observations of canopy RT fluxes at the global scale and open new perspectives for DGVM model parameterization at the regional to global scale.This paper is therefore a plea for a more rigorous treatment of surface and canopy radiation fluxes in DGVMs, which has been largely neglected in the past.Collaboration among the authors was supported by the TERRABITES COST Action ES0805 and by the International Space Science Institute (no.263).The main author was supported through the "CliSAP" (EXC177) Cluster of Excellence, University of Hamburg, funded through the German Science Foundation (DFG) and the ESA Climate Change Initative Climate Modelling User Group (contract: 4000100222/10/I-AM).Juliane Otto was funded through ERC starting grant 242564.Atmospheric radiation data were obtained from the NASA Langley Research Center Atmospheric Science Data Center.Tristan Quaife was funded by the European Space Agency (contract 4000104980/11/I-LG) with additional support from the UK National Centre for Earth Observation.The authors thank Julia Pongratz for reviewing a previous version of the manuscript.