Geoscientific Model Development Efficient modeling of sun / shade canopy radiation dynamics explicitly accounting for scattering

The separation of global radiation ( Rg) into its direct (Rb) and diffuse constituents ( Rd) is important when modeling plant photosynthesis because a high Rd:Rg ratio has been shown to enhance Gross Primary Production (GPP). To include this effect in vegetation models, the plant canopy must be separated into sunlit and shaded leaves. However, because such models are often too intractable and computationally expensive for theoretical or large scale studies, simpler sun-shade approaches are often preferred. A widely used and computationally efficient sun-shade model was developed by Goudriaan (1977) (GOU). However, compared to more complex models, this model’s realism is limited by its lack of explicit treatment of radiation scattering. Here we present a new model based on the GOU model, but which in contrast explicitly simulates radiation scattering by sunlit leaves and the absorption of this radiation by the canopy layers above and below (2-stream approach). Compared to the GOU model our model predicts significantly different profiles of scattered radiation that are in better agreement with measured profiles of downwelling diffuse radiation. With respect to these data our model’s performance is equal to a more complex and much slower iterative radiation model while maintaining the simplicity and computational efficiency of the GOU model.


Introduction
Realistic estimation of radiation profiles in plant canopies is important in order to correctly simulate processes occurring at the leaf-level, such as photosynthesis and evaporation.The inclusion of both diffuse and direct radiation transfer is important when modeling photosynthesis in plants.For example, an increased ratio of incoming diffuse (R 0,d ) to global radiation (R 0,g ) increases photosynthesis during clear-sky conditions (Spitters, 1986;de Pury and Farquhar, 1997;Alton et al., 2005;2007;Cai et al., 2009;Mercado et al., 2009).This increase is caused by shaded leaves, which are normally not light saturated, receiving an increased photon flux thus increasing photosynthesis (Roderick, 2001).Because the sunlit leaves are light saturated the simultaneous decrease in direct radiation on the sunlit leaves will not cause a large enough offset to decrease total canopy photosynthesis.
Existing canopy radiation models differ greatly in the level of detail used to simulate radiation profiles.More detailed models assume non-homogeneous canopies and/or calculate radiation interception by individual trees (e.g.Charles-Edwards and Thornley, 1973;Mann et al., 1979;Chen et al., 1994;North, 1996;Bartelink, 1998).These more complex models are able to more accurately represent radiative transfer in the canopy than simpler representations.However, because these models are complex, they cannot be solved analytically and are therefore not sufficiently tractable for use in theoretical studies (e.g.optimization models, e.g.Franklin, 2007).Perhaps the most important limitation of the complex models is their computational demand.For example, the new generation of large scale vegetation models include height structured competition for light as one of the most important processes shaping plant communities (Albani et al., 2006;Scheiter and Higgins, 2009).Because this process results in dynamic interactions between many plants within the canopy, the resulting computational demands are too high to allow the use of complex iterative or other computationally slow radiation interception models.
A relatively simple sun/shade model which can be used in global vegetation models was originally developed by Goudriaan (1977;13-40) (GOU) and later implemented by Spitters (1986).This model has been used in several studies (e.g.Anten, 1997;dePury and Farquhar, 1997;Wang and Leuning, 1998;Wang, 2003).However, a potential disadvantage of this model compared to similar but more computationally demanding models (e.g.Norman, 1979Norman, , 1980;;Sellers, 1985) is that it does not explicitly account for the scattering of direct radiation, which leads to as distribution of scattered radiation in the canopy that disagrees with predictions of the more realistic mechanistic radiation-scattering models (e.g.Norman, 1979Norman, , 1980;;Sellers, 1985).
In this study our aim is to derive a canopy radiation model with the simplicity and calculation speed of the GOU model combined with a more realistic representation of the scattering process.Taking the GOU model as a starting point we extend it by explicitly accounting for upstream and downstream fluxes of scattered radiation.Importantly, in contrast to iterative models (Norman, 1979(Norman, , 1980) ) we do this without sacrificing the analytical solvability.The model is then tested against a measured profile of downwelling scattered broadband radiation previously used to validate the model by Norman (1979Norman ( , 1980) ) in a study by Baldocchi et al. (1985).

Materials and methods
In this study we used the sun/shade model originally developed by Goudriaan (1977, 13-40 pp.) (GOU) as our starting point.This model has previously been described in detail (Spitters, 1986;Anten, 1997;de Pury and Farquhar, 1997;Wang and Leuning, 1998;Wang, 2003).The difference between our model (BF) and the GOU model lies in the treatment of scattered radiation.While the GOU model employs an implicitly defined one-stream flux, the BF model explicitly tracks two-stream generation and absorption of scattered radiation.However, for clarity a short description of the full GOU model is given below.The radiation fluxes in both models are explained in Fig. 1.

The GOU model
The radiation profile of diffuse radiation (R d ) is calculated as: (1) In Eq. (1) R 0,d is incoming diffuse radiation, ρ is canopy reflectance, k d is the extinction for diffuse radiation and L is cumulative (single sided) Leaf Area Index (LAI).Canopy reflectance is a function of solar elevation (β) and leaf scattering (σ ) (see Spitters, 1986).Leaf scattering (σ ) includes transmittance (t) and reflectance (r), which in the GOU model are assumed to be equal.The extinction coefficient for diffuse radiation (k d ) can be calculated as (Spitters, 1986): Following Lambert-Beer's law the direct radiation on sunlit leaves is assumed to be equal at all canopy depths but with the fraction of sunlit leaves (A sl ) decreasing with canopy depth: (3) The extinction coefficient for black leaves (k b ) is a function of solar elevation (β) where G is a function representing the leaf angle distribution in the canopy (e.g.Goudriaan, 1977;Baldocchi et al., 1985;Wang and Baldocchi, 1988) and a clumping index (Goudriaan, 1977;Baldocchi et al., 1985;Chen et al., 2008).When assuming a spherical leaf angle distribution, a common model assumption, the leaf angle distribution function G equals 0.5 (Goudriaan, 1977).
The clumping index ( , where 0 < ≤ 1) represents heterogeneity in canopy structure, which decreases with increasing so that = 1 represents a random distribution of leaves (no clumping).
In the GOU model, the profile of scattered radiation (R sc ) is calculated as the difference between "total direct radiation" (R b,tot ) and beam radiation.R b,tot includes beam radiation and scattered radiation generated by the reflectance and transmittance of beam radiation (Spitters, 1986).
While the beam radiation following Lambert-Beer's law is a mechanistically justified standard assumption in many models, the assumption of the "total direct radiation" also following Lambert-Beer's law is not readily justified mechanistically.Rather, this definition of "total direct radiation" represents an empirically based approximation, which leads to a profile for scattered radiation significantly different from the mechanistically derived profile predicted by BF model.

The BF model
In more complex and realistic radiative transfer models than the GOU model, scattered radiation is often treated as two separate streams: one upward stream generated by direct radiation being reflected by leaves, and one downward stream generated by transmittance of beam radiation through leaves.
The BF model is formed by replacing the original implicit treatment of scattering in the GOU by an explicit two stream approach for scattered beam radiation.Scattering of direct (beam) radiation gives rise to upstream (reflection) and downstream (transmission) fluxes of diffuse radiation.The fraction of incoming direct radiation that is transmitted at canopy depth ξ can be calculated by multiplying Eq. (3) with transmittance.This transmitted radiation is then treated as downwelling diffuse radiation.The fraction of transmitted radiation formed at ξ remaining at canopy depth L (which is below ξ ) is: Using Eqs. ( 3) and ( 6) the fraction of incoming direct radiation that has been transmitted at canopy depth ξ and that is available at canopy depth L can be calculated as: The downwelling scattered radiation at canopy depth L can be calculated by integrating Eq. ( 7) over all canopy depths between the canopy top and canopy depth L, multiplied by incoming beam radiation (R 0,b ): Following the same rationale the upward stream of scattered radiation becomes: The GOU model can now be converted into the BF model by replacing the equation for scattered radiation (Eq.5) in the GOU model by the sum of Eqs. ( 8) and ( 9): In addition to the upstream flux generated by leaf reflectance, an upward flux generated by surface reflectance can be added where W is ground-surface albedo.

Absorbed radiation
In order to be able to calculate canopy photosynthesis absorbed radiation needs to be calculated for sunlit and shaded leaves separately.Absorbed radiation by shaded leaves in the GOU model is calculated as (Anten, 1997) and absorbed radiation by sunlit leaves as: Absorbed radiation by shaded leaves in the BF model is calculated analogously as and absorbed radiation by sunlit leaves as:

Photosynthesis
Photosynthesis (µg C m −2 s −1 ) is calculated for sunlit and shaded leaves separately using a non-rectangular hyperbola function (Thornley, 2002) where R a is absorbed radiation, φ is the quantum yield (µg C J −1 ), P max (µg C m −2 s −1 ) is the light saturated gross photosynthetic rate (µg C s −1 ) and is the convexity of the relationship between P (L) and R a .P max is calculated as where n a is nitrogen content per leaf area (g N m −2 ), n min is the leaf nitrogen content for which P max is 0 (g N m −2 ), and a is the slope of the P max − n a relationship (µg C g N −1 s −1 ).

Canopy radiation profiles
To test a radiation scattering model against measured data requires measured radiation profiles where scattered and direct components are separable.Because such data are scarce modelers often use more detailed and more complex models as a benchmark instead of measured data.In this study we use a measured profile of downwelling scattered radiation (Baldocchi et al., 1985) which has previously been used to validate the model by Norman (1979Norman ( , 1980)).Using the same parameter values as reported in the study by Baldocchi (1985), we compare our model results against both a measured profile and a more complex iterative model by Norman (1979Norman ( , 1980) ) as well as the GOU model.For our model runs we use the same model assumptions as used in the Norman model (1979,1980), i.e. a spherical leaf angle distribution and no clumping ( = 1).In addition to this comparison, we apply a site specific estimate of = 0.84 (Baldocchi and Wilson, 2001) to test the effect of clumping on the predictions by our model (BF model) and the GOU model.In the study by Baldocchi et al. (1985) radiation was measured within an oak-hickory stand located near Oak Ridge, TN, USA (35 • 57 N 84 • 17 W), using sensors mounted on trams that traversed 30 m transects.They calculated a daily mean profile of downwelling diffuse radiation based on hourly-averaged data between 08:00 and 17:00 EST for Julian day 273, 1981.Hourly downwelling broadband diffuse radiation was modeled (on the half-hour) as the sum of Eqs. ( 1) and ( 8) for the BF model and Eqs. ( 1) and ( 5) for the GOU model using the same value for tand the fraction of diffuse to global radiation (fDif ) as reported in Baldocchi et al. (1985) (0.22 and 0.17, respectively), with the value of β calculated using a two-dimensional geometry (Appendix A).Broadband radiation was used as it is the only spectral range for which profiles for downwelling diffuse radiation was available.In the GOU model σ was set to equal t in order to only account for downwelling scattered radiation.The modeled broadband radiation profiles were then compared against the measured profile taken from Baldocchi et al. (1985).As the same model parameters and model assumptions (spherical angle distribution) were used in our model runs as in the model runs using the Norman (1979Norman ( , 1980) ) model (Baldocchi et al., 1985), we can also compare our results against this modeled radiation profile.
Modeling Efficiency (ME: Appendix B) of the simulated profiles compared to measurements was also calculated.
In addition to the testing of model results against the measured profile (above) we also compare the fluxes of scattered radiation for the respective models and spectral bands (Broadband, photosynthetically active radiation (PAR), near infrared (NIR)).For the GOU model Eq. ( 5) was used and for the BF model Eqs.( 8) and ( 9).The same values for β and fDif as above were used and for reflectance (r) and transmittance (t) we used the values reported in Baldocchi et al. (1985) (Broadband: r = 0.30 and t = 0.22; PAR: r = 0.11 and t = 0.16; NIR: r = 0.43 and t = 0.26).As the GOU model does not separate scattering into reflectance and transmittance, scattering (σ ) for this model was set to the sum of r and t.

Effects on GPP
To investigate the effect of our modifications to the GOU model on simulated GPP we perform a test where we use a range of incoming radiation between 50 and 800 W m −2 and a solar elevation ranging from 45 • to 80 • .The same LAI (5.5) and value for fDif as in Baldocchi et al. (1985) were used and no clumping was assumed ( = 1).In this test we assume φ to be 2.73 µg C J −1 ; to be 0.75 and n min 0.4 g N m −2 (Franklin and Ågren, 2002).Leaf nitrogen content (n a ) was assumed to be 2.3 g N m −2 and a to be 65.7 µg C g N −1 s −1 based on values for Q. rubra (Reich et al., 1995).

Model differences in the radiation profile
The difference between the GOU and BF models lies in the treatment of scattered radiation (R sc ), where the GOU model treats R sc as the difference between the implicitly (non-mechanistically) defined total direct radiation (R b,tot ) and beam radiation (R b,b ).Contrastingly, the BF model calculates R sc explicitly and mechanistically accounting for all processes included in the radiation interception model (Fig. 1).Consequently, the profiles of total scattered radiation differ markedly between the models (Fig. 2) at the top of the canopy.This difference is largest in the photosynthetically active radiation (PAR) spectra (Fig. 2c), and smallest in the near infrared (NIR) spectra (Fig. 2b).Total scattered radiation in the GOU model follows an exponential extinction in line with Lambert-Beer's law.In the BF model the shapes of the curves for upward and downward scattered radiation differ.Upward scattered radiation (R sc↑ ) follows a profile similar to the GOU model, whereas downwelling scattered radiation (R sc↓ ) is 0 at the canopy top and increases to a maximum at cumulative leaf area L ∼ 1.5 followed by a shallow decline.The resulting hump shaped profile for total scattered  radiation in the BF model follows from the competing effects of the accumulation of generated scattered radiation and the cumulative absorption of this radiation down the canopy.

Comparison with a measured radiation profile
Tested against measurements of downwelling diffuse radiation, the BF model performs significantly better than the original GOU model (Fig. 3) with ME = 0.46 compared to −2.45 for the GOU model (excluding the trivial point L = 0).Notably, the radiation profile for downwelling diffuse radiation predicted by the BF model is almost identical to the one using the more complex Norman model (Norman, 1979(Norman, , 1980)).By using a clumping index as reported at the same site (Baldocchi and Wilson, 2001) the model fit increased for both models (Fig. 3).For the GOU model the fit became only slightly better with ME = −2.06,whereas we get ME = 0.63 for the BF model.

Model difference in simulated GPP
The difference between the GOU and BF models in simulated GPP varied with both solar elevation and incoming radiation (R 0 ), with the latter being most important.GPP is in general larger for the GOU model; for a LAI of 5.5 the difference in total canopy GPP varies between 0.1 and 3.1 % (on average 1.4 %) (Fig. 4).For a low R 0 (below 150 W m −2 ) the difference was always larger than 1.8 %.As the difference in scattered radiation was largest at the top of the canopy we also compare simulated GPP for a LAI of 2 (Fig. 4).Here the difference between the GOU and BF model varied be- tween 1.2 and 6.2 % (on average 3.1 %), with the difference for low R 0 always being larger than 4.4 %.

Implications for canopy radiation modelling
Given the prominent role of light in shaping vegetation through responses at scales from leaf physiology to community dynamics, an accurate representation of canopy light absorption is important in vegetation modelling.Not only 1 Figure 4. Difference in simulated GPP (%) between the GOU and BF model using a LAI of 2 2.0 (thin lines) or 5.5 (thick lines) using a solar elevation of 45 ○ (solid lines) or 80 ○ (dotted 3 lines).4 Fig. 4. Difference in simulated GPP (%) between the GOU and BF model using a LAI of 2.0 (thin lines) or 5.5 (thick lines) using a solar elevation of 45 • (solid lines) or 80 • (dotted lines).
does light absorption limit photosynthesis but it also controls leaf N concentration (Franklin and Ågren, 2002) and water use (Buckley et al., 2002) with implications for ecosystem resource balances.It is therefore not surprising that much effort has been put into construction of realistic canopy radiation models.Because realism often comes at a computational cost and loss of tractability there is also a large demand for simple and computationally efficient models, such as the widely used Goudriaan model (Spitters, 1986;Anten, 1997;dePury and Farquhar, 1997;Wang and Leuning, 1998;Wang, 2003).
Our replacement of the original one-stream implicit scattered radiation flux in the Goudriaan approach by an explicit two-stream model of radiation scattering significantly and qualitatively changed the shape of the modelled scattered radiation profile compared to the original GOU model (Fig. 2).The predictions of the new BF model were in better agreement with measured canopy radiation profiles of diffuse downwelling radiation than the original GOU model (Fig. 3).Furthermore, our model results are very close to the more complex and computationally intensive Norman model (Norman, 1979, 1980) while maintaining the high computational efficiency of the GOU model.This similarity in results suggests that the BF model could be used as a computationally efficient approximation of the Norman model.
Given the importance of canopy radiation modeling, our qualitative as well as quantitative improvement of a tractable analytical canopy radiation model can be used to improve a wide range of plant, vegetation and ecosystem models.While the effect of our modification of the GOU model on simulated GPP is modest at high LAIs, at lower LAIs the effect is larger, up to 6.2 %, in our simulation of an oak/hickory forest canopy.Because GPP is the basis for most plant and ecosystem processes even small improvements in its predic-tions may be valuable.For example, the new generation of large scale vegetation models include height structured competition for light as one of the most important processes shaping plant communities but do not explicitly account for radiation scattering (e.g.Smith et al., 2001;Albani et al., 2006;Scheiter and Higgins, 2009).Such vegetation models could readily be improved at very low computational cost by adding a simple canopy scattering approach, such as the BF model.

1Figure 1 .Fig. 1 .
Figure 1.Diagram of the radiation interception in the two different models: GOU (Goudriaan, 2 1977) and our new model (BF).Both the GOU (a) and BF model (b) simulate the fraction of 3 sunlit leaves (white areas) to shaded leaves (grey areas) in the same way according to 4 Lambert-Beer's law.The absorption of incoming diffuse radiation (broken grey arrows) is 5 also simulated the same way with incoming diffuse radiation being absorbed by both sunlit 6 and shaded leaves also following Lambert-Beer's law.The difference between the models lies 7 in the fluxes of scattered radiation, which in the BF model are generated by the transmittance 8 and reflectance of beam radiation (black arrows) at all levels of the canopy.The GOU model 9 approximates the total interception of incoming direct radiation (broken black arrow) with a 10 single aggregated flux that includes both the interception of beam radiation and the scattered 11 radiation generated by the beam radiation.Scattered radiation is then calculated at each level 12 as the difference between total direct radiation (black arrow) and beam radiation (based on the 13 area of sunlit leaves).In the BF model, scattered radiation is generated by beam radiation 14 (black arrow) being intercepted by sunlit leaves and scattered into one up-and one 15 downstream (solid grey arrows) of scattered (diffuse) radiation.16 17 18 19 20 21

Figure 2 .
Figure 2. The modeled profiles of scattered radiation (R sc ) in the broadband (a), near infra red (b) and photosynthetically active radiation spectra (c) relative incoming beam radiation (R 0,b ) for the respective models.Note the difference in scale between a-b and c.For the BF model, scattering (thick solid line) is divided into upwelling (BF sc↑ : thin solid line) and downwelling scattered radiation (BF sc↓ : dashed line).For the GOU model total scattered radiation is simulated (GOU sc : dotted/dashed line) using scattering equal to the sum of transmittance (t=0.22) and reflectance (r=0.30).For the BF model scattered radiation (BF sc : thick solid line) is the sum of BF sc↓ and BF sc↑ .

Fig. 2 .
Fig. 2. The modeled profiles of scattered radiation (R sc ) in the broadband (a), near infra red (b) and photosynthetically active radiation spectra (c) relative incoming beam radiation (R 0,b ) for the respective models.Note the difference in scale between (a-b) and (c).For the BF model, scattering (thick solid line) is divided into upwelling (BF sc↑ : thin solid line) and downwelling scattered radiation (BF sc↓ : dashed line).For the GOU model total scattered radiation is simulated (GOU sc : dotted/dashed line) using scattering equal to the sum of transmittance (t = 0.22) and reflectance (r = 0.30).For the BF model scattered radiation (BF sc : thick solid line) is the sum of BF sc↓ and BF sc↑ .