Fast matrix treatment of 3-D radiative transfer in vegetation canopies: SPARTACUS-Vegetation 1.1

. A fast scheme is described to compute the 3-D interaction of solar radiation with vegetation canopies. The canopy is split in the horizontal plane into one clear region and one or more vegetated regions, and the two-stream equations are used for each, but with additional terms representing lateral exchange of radiation between regions that are proportional to the area of the interface between them. The resulting coupled set of ordinary differential equations is solved us-ing the matrix-exponential method. The scheme is compared to solar Monte Carlo calculations for idealized scenes from the “RAMI4PILPS” intercomparison project, for open forest canopies and shrublands both with and without snow on the ground. Agreement is good in both the visible and infrared: for the cases compared, the root-mean-squared difference in reﬂectance, transmittance and canopy absorptance is 0.020, 0.038 and 0.033, respectively. The technique has potential application to weather and climate modelling.


Introduction
The treatment of the interaction of vegetation with solar radiation in weather and climate models varies greatly in complexity.The simplest schemes are concerned only with surface albedo and its impact on near-surface temperature forecasts, and indeed Viterbo and Betts (1999) reported a large improvement in forecasts by the ECMWF model when the use of a fixed snow albedo was modified to account for the much lower albedo that occurs when snow falls in forested areas.Much more sophisticated treatments are used in the dynamic vegetation schemes of many climate models, which need to calculate also the fraction of absorbed photosyn-thetically active radiation (faPAR).But it was reported by Loew et al. (2014) that even state-of-the-art models, when evaluated in benchmarks for which a full physical description of the vegetation was available, had worst-case albedo errors in excess of 0.3.The challenge is to represent the complex 3-D structure of vegetation canopies with a radiative transfer algorithm that is nonetheless computationally efficient enough to use in a global model.Sellers (1985) took the two-stream equations used in atmospheric radiative transfer and applied them to a vegetation canopy.In this approach, the vegetation is treated as a single horizontally homogeneous layer, and a set of three coupled ordinary differential equations are solved for the direct downwelling irradiance and the downwelling and upwelling diffuse irradiances.If the leaves can be assumed randomly oriented then the optical depth of the layer is equal to half the leaf area index (LAI).Meador and Weaver (1980) provided an analytic solution to these equations that is still used in a number of state-of-the-art surface energy exchange schemes (e.g.Best et al., 2011).The first-order error that arises is due to the fact that vegetation canopies are not homogeneous: the heterogeneous distribution of leaves within a tree crown and crowns within a forest stand is such that leaves are more likely to be shadowed by other leaves than if they were homogeneously distributed.Typically this is treated by introducing a "clumping factor" that scales down the LAI used in the two-stream scheme.A very similar approach has previously been used in atmospheric radiation schemes to treat the clumpiness of clouds (Tiedtke, 1996).The clumping factor for vegetation is typically parameterized as an empirical function of properties of the vegetation and solar zenith angle (e.g.Ni-Meister et al., 2010), but this Published by Copernicus Publications on behalf of the European Geosciences Union.lacks a physical basis and fails to represent horizontal fluxes into and out of individual tree crowns.Pinty et al. (2006) described one of the most sophisticated yet affordable schemes to date that attempts to overcome these limitations.Their scheme sums three terms: the reflection from the vegetation assuming a black underlying surface, the reflection from the surface assuming no interaction with the vegetation, and a term representing interactions between the surface and the vegetation.Despite much improved performance compared to the Sellers (1985) scheme, their approach still uses an empirical clumping factor, and is underpinned by the Meador and Weaver (1980) solution that assumes horizontally homogeneous vegetation.
In this paper we exploit recent advances in the atmospheric literature, and adapt the "SPARTACUS" (SPeedy Algorithm for Radiative Transfer through CloUd Sides) method of Hogan et al. (2016) to the vegetation problem.As described in Sect.2, this approach employs an explicit description of the horizontal distribution of vegetation for which we can write down a modified version of the two-stream equations that includes terms for lateral radiation exchange between tree crowns and the clear regions between them.The equations are then solved exactly using the matrix-exponential method.This avoids the need for an empirical clumping factor or the Meador and Weaver (1980) solution.In Sect. 3 it is compared to Monte Carlo calculations in idealized forest and shrubland conditions.

Overview
We use a simple geometrical description of the problem, as shown in Fig. 1.Leafy vegetation is assumed to occupy a single constant-thickness "canopy layer", with the horizontal domain (corresponding to a weather-or climate-model grid box) divided into m "regions".Within an individual region, the optical properties of the atmosphere and any vegetation are assumed horizontally and vertically homogeneous.c).The use of two vegetated regions adds the flexibility to represent horizontally heterogeneous tree crowns and trees of differing leaf density, borrowing the idea of Shonk and Hogan (2008) for representing cloud heterogeneity.In Sect. 3 we compare this to a simpler tworegion approach with only one vegetated region (denoted b).While the tree crowns are depicted in Fig. 1 as cylinders, this is not explicitly assumed; rather, we assume that (1) all azimuthal orientations of the interface between the clear and vegetated regions are equally likely, and (2) the tree crowns are randomly distributed.To represent forests with a significant separation between the ground and the base of the tree crowns, an additional "sub-canopy layer" may be added, also divided into m regions (see Fig. 1).Thus we require as a min-imum just four numbers to define the geometry of the problem: the fractional area of the domain covered by vegetation, c v , the vertical depth of the canopy layer, z 1 , the vertical depth of the sub-canopy layer, z 2 (which may be zero), and the length of the interface between the clear and vegetated regions per unit area of the domain, L ab .Note that although this paper considers only up to two layers and three regions, which is an appropriate level of complexity for a weather or climate model, for other applications additional layers and regions may be added.This would enable the representation of different types of vegetation of different heights, or vegetation in the understory.
In the SPARTACUS method, the two-stream differential equations are used in each region, but with additional terms representing lateral radiation transport between regions.The formulation of these equations is given in Sect.2.2, with the coefficients to be used in the case of vegetation defined in Sect.2.3.Section 2.4 then outlines how the L ab term could be parameterized in a model.Section 2.5 describes how the equations are solved for a single layer using matrix exponentials, and Sect.2.6 describes the use of the adding method to compute the direct and diffuse albedos of the entire scene (vegetation and the surface beneath it).In the context of a weather or climate model, this could be done for the same spectral intervals as the atmospheric radiation scheme, or in the smaller number of broader spectral intervals for which optical properties of the vegetation and surface are defined.These albedos would then be used as boundary conditions for the calculation of the radiative flux profile in the atmosphere above.The downwelling direct and diffuse irradiances output from the atmospheric radiation scheme are then used in Sect.2.7 to compute the irradiance profile within the vegetation canopy, enabling the absorbed and transmitted radiation to be computed.The Appendix describes how the scheme may be made computationally faster by optimizing the treatment of the sub-canopy layer.

Differential two-stream equations in matrix form
This section summarizes the theoretical background to SPARTACUS that was introduced by Hogan et al. (2016).Solar radiation in a particular spectral interval is described by three streams: the diffuse upwelling irradiance (u), the diffuse downwelling irradiance (v) and the direct downwelling irradiance (s), where u and v are irradiances into a horizontal plane while s is into a plane oriented perpendicular to the Sun.At any given height, these are column vectors containing the irradiances in m regions; in the equations that follow we use m = 3 to match the schematic shown in Fig. 1, but it is straightforward to reduce to two regions.Thus for upwelling irradiance we have u = ( u a u b u c ) T , where each irradiance component is defined as the radiative power divided by the area of the entire grid box, such that the domain-mean irradiance is obtained by summing the elements of the vector.
The two-stream equations form a set of coupled differential equations that can be written in matrix form as where z is height measured downward from the top of the layer, and is a matrix describing the interactions between irradiance components and between different regions.It is convenient to partition it into a set of m × m component matrices as follows: where and 4 is the same as 3 but using the quantity γ 4 in place of γ 3 .Missing entries in all these matrices are taken to be zero.
The 0 and 1 matrices describe the rate at which the direct and diffuse downwelling irradiances, respectively, change along their path.They are expressed in Eqs. ( 3) and ( 4) as the sum of two matrices: the first matrix in each case represents losses due to scattering and absorption, while the second represents exchange of radiation between regions.The 2 matrix describes the rate of scattering of diffuse radiation from one direction to the other, while the 3 and 4 matrices describe the rate at which the direct solar beam is scattered into the upwelling and downwelling diffuse streams.The minus signs in front of the matrices on the top row of Eq. ( 2) are due to this line corresponding to upwelling radiation, but the vertical coordinate increasing downward.
The symbols in Eqs. ( 3)-( 6) have the following meanings.The extinction coefficient to diffuse radiation of region j is denoted σ j , and σ j 0 is the same but for direct radiation.The distinction between the two permits the flexibility to represent leaves with a preferred orientation.The cosine of the solar zenith angle is denoted µ 0 while the single-scattering albedo is ω.The coefficients γ 1 -γ 4 govern the exchange of radiation between the three streams.Finally, the coefficients f j k dir and f j k diff represent the rate at which direct and diffuse radiation, respectively, is transferred from region j to region k.All these symbols are defined in terms of physical properties of the scene in the next section.

Coefficients in the two-stream equations
The matrix form of the two-stream equations in Sect.2.2 introduced several coefficients that are themselves functions of more fundamental optical or geometric properties.The γ 1 -γ 4 coefficients may be written as (Meador and Weaver, 1980) www.geosci-model-dev.net/11/339/2018/Geosci.Model Dev., 11, 339-350, 2018 where β and β 0 are the "upscatter" fractions, the fractions of downwelling radiation (in the diffuse and direct streams respectively) that are scattered upward, and µ 1 is the cosine of the effective zenith angle of diffuse radiation.For the remainder of this paper we assume the diffuse radiation to be hemispherically isotropic, so µ 1 = 1/2.
In the simplest case, where leaves are assumed to be randomly oriented, the optical depth of a region is equal to half its LAI, and therefore for a layer of thickness z, the extinction coefficients to direct and diffuse radiation are the same and are given by Assuming the leaves to be bi-Lambertian scatterers with reflectance r and transmittance t, the single-scattering albedo is given by and the upscatter fractions by These last two formulas may be derived by equating Eqs. ( 8) and ( 9) with the definitions given in the lowest row of Table 4 of Pinty et al. (2006).Pinty et al. (2006) also provided more general expressions for leaves with a preferential alignment.The rates of lateral exchange of radiation between regions that appear in Eqs. ( 3) and (4) may be derived from geometrical arguments (Hogan and Shonk, 2013;Schäfer et al., 2016) as where θ 0 is the solar zenith angle, L ij is the length of the interface between regions i and j per unit area of the horizontal domain, and c i is the fractional area of the domain covered by region i.In the m = 3 case we have two regions to represent horizontal heterogeneity of zenith optical depth, and following the findings of Shonk and Hogan (2008) we assume them to be of equal area, i.e. c b = c c = c v /2 and c a = 1 − c v (where c v is the fractional coverage of vegetation).This leads to f bc dir = f cb dir and f bc diff = f cb diff .We may represent the effect of vertical tree trunks in region c of the sub-canopy layer (illustrated in Fig. 1) as follows.If the trunks are of a size and number such that a horizontal slice through the sub-canopy layer intercepts a normalized total trunk perimeter (per unit area of region c) of L t , then by analogy with Eqs. ( 15) and ( 16), the diffuse and direct extinction coefficients are given by For simplicity we assume the trunks to be Lambertian reflectors, in which case ω is simply the trunk albedo, and with no preference for upward or downward scattering we have β = β 0 = 1/2.Now that the problem has been formulated mathematically, we can explain how the assumption that the tree crowns are randomly distributed is implicitly encoded in the equations.At any given height in the canopy layer, the probability of direct radiation in the clear region intercepting a tree crown, per unit distance travelled vertically, is f ab dir .This factor is constant in the canopy layer.Therefore, for direct radiation emerging unscattered from the edge of a tree crown into the clear region, the fraction of that light remaining in the clear region rather than having encountered another tree varies in proportion to exp(−f ab dir z), where z is the vertical distance travelled in the clear region (assuming no absorption or scattering, and that the light remains within the canopy layer).To express this in terms of horizontal distance x, we use Eq. ( 16) and recognize that tan(θ 0 ) = x/z to obtain exp[−xL ab /(π c a )].This implies that the chord lengths between the edges of tree crowns in all possible horizontal directions also follow the same exponential distribution, which in turn defines the spatial distribution of trees as random.

Parameterizing the vegetation perimeter length
The length of the vegetation-clear boundary, L ab , is the fundamental property used by SPARTACUS to characterize the importance of lateral radiative exchange between clear and vegetated regions.It is therefore the quantity that would ideally be measured in field experiments.However, in the context of weather and climate modelling, the physiographic variable available would most likely be vegetation cover c v (e.g. from the measurements of Hansen et al., 2003), and L ab would need to be parameterized as a function of c v .This can be done by introducing an extra parameter representing the characteristic size of a tree crown that is independent of c v .We now present two possible characteristic sizes that could be used.
In the first case, we define the effective tree diameter, D, to be the diameter of identical, cylindrical and physically separated tree crowns in an idealized forest with the same L ab and c v as the real forest.The assumption that tree crowns do not touch was used by Widlowski et al. (2011) in generating the idealized scenes that we use in Sect. 3 to evaluate SPARTACUS.The phenomenon of the crowns of some tree species remaining separate even for large tree cover is known as crown shyness (e.g.Putz et al., 1984).In analogy to the concept of an effective cloud diameter by Jensen et al. (2008), this leads to the definition If region c represents the central core of the tree crowns, as depicted in Fig. 1, then this implies Geosci.Model Dev., 11, 339-350, 2018 www.geosci-model-dev.net/11/339/2018/ In the second case we assume that tree crowns can touch each other, and will do so increasingly in dense forests.This behaviour is represented by defining an effective tree scale, S, such that This form is inspired by the idealized geometrical analysis of Morcrette ( 2012): if we place idealized trees with a square footprint measuring S × S randomly on a grid, then on average the normalized perimeter length L ab will follow Eq. ( 20).
It leads to the behaviour that L ab increases with c v up to c v = 1/2, but for further increases in c v , crown touching dominates which causes L ab to reduce again.
In the field we would envisage measuring L ab and c v and then using Eqs.( 19) and ( 20) to infer D and S. The characteristic size that varies least with c v would then be the one best suited for use in a weather or climate model, and potentially a constant characteristic size could be used to characterize an entire forest on a regional scale.Within individual grid boxes of the model, it would be used to compute L ab from c v using either Eqs.( 19) or (20).

Solution to equations within one layer
We may write the solution to Eq. ( 1) in terms of a matrix exponential (Waterman, 1981;Hogan et al., 2016): the irradiances at the base of a layer of thickness z are related to the irradiances at the top of the layer via where the matrix exponential may be computed numerically using the scaling and squaring method (e.g.Higham, 2005).If 3-D radiative transfer is neglected then f diff = f dir = 0, which decouples the equations to the extent that a computationally cheaper analytical solution is possible (Meador and Weaver, 1980).Conversely, if scattering and absorption are ignored but 3-D radiative transfer is retained, a reasonable assumption in the sub-canopy layer, then σ = σ 0 = 0, which also decouples the equations and leads to the computationally cheaper solution given in the Appendix.
In order to compute the irradiance profile, we wish to work with expressions of the following form: where Eq. ( 22) states that the upwelling irradiance exiting the top of the layer is equal to transmission of the upwelling irradiance entering the base of the layer, plus reflection of the downwelling irradiance entering the top of the layer, plus scattering of the direct solar irradiance entering the top of the layer, and similarly for Eq. ( 23). Figure 1 illustrates the meaning of the elements of the diffuse reflectance matrix R for the canopy layer: where R ij is the fraction of diffuse downwelling radiation entering the top of region i that is scattered out of the top of region j without exiting the base of the layer.The other matrices have analogous definitions: T represents the transmission of diffuse radiation across the layer, and S + and S − represent the scattering of radiation from the direct downwelling stream at the top of the layer to the diffuse upwelling stream at the top of the layer and the diffuse downwelling stream at the base of the layer, respectively.These matrices may be derived from the matrix exponential, which we decompose into seven m × m matrices: It was shown by Hogan et al. (2016) that Moreover, the direct irradiance exiting the base of a layer is computed from the direct irradiance entering the top of a layer via s(z + z) = E 0 s(z).

Extension to multiple layers
To compute the irradiance profile we use the adding method (Lacis and Hansen, 1974) but in a somewhat different form to Hogan et al. (2016), in order to facilitate integration within a full atmospheric radiation scheme.This section considers the first part: stepping up through the vegetation layers computing the albedo of the scene below each layer interface.We define the matrix A i+1/2 as the albedo to diffuse downwelling radiation of the scene below interface i+1/2 (including the surface contribution), and the matrix D i+1/2 as the albedo to direct radiation.The off-diagonal terms of these matrices represent the fraction of radiation downwelling in one region that is reflected back into the other.At the surface (interface n + 1/2 for an n-layer description of the canopy), these matrices are diagonal: where for maximum flexibility we allow for separate direct and diffuse surface albedos, and separate albedos below each region to represent lower snow cover beneath trees.
We then use the adding method to compute A and D just below the interface above, accounting for the possibility of multiple scattering.In the case of the diffuse albedo matrix we have where I is the m × m identity matrix.This equation states that the albedo at interface i − 1/2 is equal to the reflection of layer i, plus the albedo at interface i + 1/2 accounting for the two-way transmission through the intervening layer.The term in square brackets accounts for multiple scattering between interface i +1/2 and layer i, and since it is a geometric series of matrices, the equation reduces to Similarly, the direct albedo matrix at the interface above is given by where D i+1/2 E 0i represents the direct radiation that passes down through layer i without being scattered and is then reflected up from interface i + 1/2, while A i+1/2 S − i represents direct radiation that is scattered into the downward diffuse stream in layer i and then reflected up from interface i+1/2.For the two-layer description of the vegetation shown in Fig. 1, Eqs. ( 33) and ( 34) are applied first at interface 1.5 (between the canopy and the sub-canopy layers) and then at interface 0.5 (the top of the canopy).It is straightforward to add additional layers.
At this point we are able to compute the scalar "scene albedos" of the surface and the vegetation.Denoting c = ( c a c b c c ) T as a column vector containing the area fractions of each region, the scene albedos to diffuse and direct radiation are When implementing the scheme described in this paper in the radiation scheme of a weather or climate model, these albedos would be used as the boundary conditions for the computation of the irradiance profile through the atmosphere.

Computing irradiances within the canopy
After running the atmospheric part of the radiation scheme, we proceed down through the vegetation to compute the direct and diffuse irradiances at each interface, ending up at the surface.The output from the atmospheric radiation calculation includes the downwelling direct and diffuse irradiances at the top of the canopy, s 1/2 and v 1/2 .These are partitioned into component irradiances at the top of each region according to the area fraction of each region: The direct irradiance is propagated down through the vegetation simply with The diffuse irradiances at the interface beneath satisfy Thus, application of Eq. ( 42) followed by Eq. ( 40) provides the irradiances at the interface below.
The horizontally averaged upwelling diffuse, downwelling diffuse and downwelling direct irradiances at interface i + 1/2, denoted u i+1/2 , v i+1/2 and s i+1/2 , respectively, are found by simply summing the elements of u i+1/2 , v i+1/2 and s i+1/2 .The total downwelling irradiance is then the sum of the direct and diffuse components: , 11, 339-350, 2018 www.geosci-model-dev.net/11/339/2018/v i+1/2 .The solar absorption by each layer is the difference in net irradiance between the interface above and below it.These definitions are used to compute normalized quantities that will be used to evaluate SPARTACUS in Sect.3.

Results
To test the application of the SPARTACUS methodology to the vegetation problem, we use two 3-D scenarios from the RAMI4PILPS1 intercomparison exercise (Widlowski et al., 2011).The first scenario is an idealized representation of an open forest canopy, and consists of spheres of leafy vegetation of diameter 10 m, while the second represents shrubland and consists of spheres of diameter 1 m.Details are provided in Table 1, including the three different area coverages of vegetation that are used.Two spectral intervals are simulated, representing the photosynthetically active visible region and   Figure 2 shows the results for the open forest canopy in the visible part of the spectrum while Fig. 3 shows the same but for the near-infrared.The corresponding results for the shrubland scenario are shown Figs. 4 and 5. Using the domain-mean irradiances defined in Sect.2.7, the quantities shown are reflectance R, transmittance T and absorptance A: It can be seen that the three-region version of SPARTACUS compares well to Monte Carlo, including all four combinations of high-and low-reflectance leaves over a high-or lowreflectance surface.In total we have 72 points of comparison with Monte Carlo calculations: two scenarios, two spectral intervals, two surface types, three vegetation covers and three solar zenith angles.Treating the Monte Carlo as "truth", we compute that the root-mean-squared error in R, T and A is 0.020, 0.038 and 0.033, respectively.Probably the worse performance occurs for low solar zenith angle in Fig. 2f (corre-sponding to visible radiation illuminating a scene with a tree cover of 0.5 over snow): A is overestimated by around 0.05, suggesting that a little too much reflected sunlight from the snow enters the tree crowns and is absorbed.
We next investigate how the results are degraded when using a more approximate description of the scene.Each panel of Figs.2-5 includes two further lines.The "homogeneous" calculation uses the same SPARTACUS code but with only one region, treating the canopy as a single horizontally homogeneous layer with the same leaf area index.This is essentially the same as the Sellers (1985) assumption and indeed with a single region the matrix-exponential method yields the same result as the Meador and Weaver (1980) solution.We see immediately that when the leaves are not clumped into trees but rather distributed uniformly, their exposure to incoming radiation is maximized and their absorptance is overestimated by up to 0.3.Conversely, both the reflectance and transmittance of the scene are underestimated, with the largest error in reflectance for overhead sun and a snow-covered surface (Fig. 2e).
The two-region SPARTACUS calculation shown in calculation in the near infrared, but absorption still tends to be overestimated in the visible.An analogous bias occurs in cloudy radiative transfer calculations in which the internal variability of clouds is neglected, which led to the proposal of Shonk and Hogan (2008) to use three regions to represent a partially cloudy scene.The success of the three-region approach suggests that it is also useful for vegetation.Having said this, the uncertainty in computing radiative transfer the vegetation canopies of weather and climate models is typically dominated by uncertainties in leaf area index.Therefore, for many applications the two-region calculation would be adequate.Since the computational cost of SPARTACUS is dominated by the matrix exponential calculation, whose cost is approximately proportional to m 3 , we would expect a tworegion SPARTACUS calculation to be at least 3 times faster than a three-region calculation.

Discussion and conclusions
This paper has demonstrated the potential for the interaction of solar radiation and complex vegetation canopies to be represented via an explicit description of the geometry, building on the SPARTACUS algorithm for representing the 3-D radiative effects of clouds (Hogan et al., 2016).The two-stream equations are written down for the tree crown and the gaps between them, but with additional terms for the horizontal exchange of radiation between regions.The equations are solved exactly using the matrix exponential method.Multiple layers are possible, although we have simplified the original SPARTACUS algorithm by assuming maximum overlap between the regions in each layer, rather than the arbitrary overlap considered by Hogan et al. (2016).Comparison against Monte Carlo calculations from the RAMI4PILPS intercomparison exercise indicates that canopy reflectance, transmittance and absorptance are computed significantly more accurately than a number of state-of-the-art models assessed by Loew et al. (2014).An advantage of the SPARTACUS approach is that in addition to LAI, only a handful of physiographic variables are required to describe the geometry of the vegetation, such as the vegetation height, coverage, and the diameter of typical tree crowns.Global estimates of the first two are now available from satellites (e.g.Simard et al., 2011;Hansen et al., 2003).
Although the testing scenarios used in this papers were simple homogeneous spheres with no woody material, the method described has the capability to represent more complex geometries.Horizontal variations in leaf density or tree crowns with different properties may be represented via two or more vegetated regions with distinct optical properties.This paper considered a two-layer description of the vegetation, with a single canopy layer overlying a sub-canopy layer, but the equations can easily be applied to a multi-layer description of the canopy, for example to compute the vertical profile of absorbed photosynthetically active radiation.The optical effects of tree trunks may also be incorporated.Moreover, the good performance with solar radiation suggests that the thermal-infrared version of SPARTACUS (Schäfer et al., 2016) could also be adapted to the vegetation problem.
A further possible extension to SPARTACUS would be to use it for remote sensing; in addition to the possibility of more accurate LAI retrievals via explicit treatment of 3-D radiative effects, this would provide a consistent framework for both remote sensing and weather/climate modelling.The challenge would be to adapt SPARTACUS to compute solar radiances rather than irradiances, which adds an extra degree of geometrical complexity.For example, trees cast shadows on the ground, but the extent to which shadows are visible to a satellite depends on the sensor zenith angle and the azimuthal separation of the sensor and the Sun.
Figure 1 considers three regions: one clear (denoted a) and two vegetated (denoted b and

Figure 1 .
Figure 1.Schematic of the idealized vegetation considered in this paper, illustrating the meanings of Layers 1 and 2 and Regions a, b and c.The diagram on the right also illustrates the interpretation of the elements of the reflectance matrix R given in Eq. (24).

Figure 2 .
Figure 2. Comparison of normalized irradiances vs. solar zenith angle for the RAMI4PILPS "Open forest canopy" scenario with optical properties appropriate for visible radiation.The two rows of panels show results for different surface albedos (α) with the top row using values appropriate for a snow-free surface and the bottom row using values for a snow-covered surface.The columns represent different areal tree fractions (c v ).The three solid lines depict the reflectance, transmittance and absorptance defined in Eqs.(43), (44) and (45), computed using the three-region version of SPARTACUS.The dashed and dot-dashed lines depict the two-region and one-region SPARTACUS calculations, respectively, where the latter involves complete horizontal homogenization of the vegetation properties through the domain.Also shown are the corresponding Monte Carlo calculations of Widlowski et al. (2011) at solar zenith angles of 27, 60 and 83 • .

Figure 3 .
Figure 3.As Fig. 2 but with optical properties appropriate for near-infrared radiation.
Figs. 2-5 treats individual trees as horizontally homogeneous cylinders, thereby neglecting the variation in zenith optical depth of the spherical trees simulated by the Monte Carlo calculations.The results are much better than those with just a single region, and virtually the same as the three-region www.geosci-model-dev.net/11/339/2018/Geosci.Model Dev., 11, 339-350, 2018

Table 1 .
Widlowski et al., 2011)e geometry of "Open forest" and "Shrubland" RAMI4PILPS scenarios simulated in this paper (seeWidlowski et al., 2011).The leaf area index of a vegetated region is defined as the total leaf surface area divided by the downward projected area of the region.

Table 2 .
Widlowski et al., 2011)e optical properties of the leaves and the surface in the visible and near-infrared in the RAMI4PILPS cases (seeWidlowski et al., 2011).= 3) version of SPARTACUS.The two vegetated regions (b and c) are of equal projected area and are configured to approximate the distribution of zenith optical depth of spheres.So for a sphere of radius r, region c represents the upper half of the optical depth distribution corresponding to a core of radius r/