Chapter 2: Radiation

Planck’s Law and the Stefan-Boltzmann equation

All objects with a temperature above absolute zero emit radiation and one of the significant advances in modern physics was the discovery by the German physicist, Max Planck, of a model for predicting how this works. Planck’s model predicts how the radiant spectral flux density (\(\small E_{b}(λ,T)\) in W·m-3), from a black body radiator changes as a function of absolute temperature (\(\small T\)) in Kelvin and its wavelength (\(\small \lambda\)) in m

\[\begin{equation} E_{b}(λ,T)=\frac{2\pi h_p~c^2}{\lambda^5\lgroup\exp\lgroup\frac{h_pc}{k\lambda T}\rgroup-1\rgroup} \end{equation}\]

where \(\small h_p\) is Planck’s constant (6.6256 × 10^-34 J·s), \(\small k\) is the Boltzmann constant (1.38054 × 10^-23 J·K^-1) and \(\small c\) is the speed of light in a vacuum (2.997925 × 10⁸ m·s^-1).

This is useful for a number of reasons. Firstly, if one knows the temperature of an object one can predict the spectral properties of radiation emitted. So for example, it is used to describe the spectral properties of radiation emitted by the sun (c. 5778K) or the earth (c. 288K) as produced by the example code below using the function Planck (which takes wavelength in nm and temperature in degrees C as inputs). It can be seen that the majority of radiation emitted by the sun is in the visible range, though a significant amount is also emitted in the ultraviolet and near-infrared ranges. Collectively, radiation emitted by the sun is referred to as shortwave radiation. By contrast radiation emitted by the earth, referred to as thermal or longwave radiation, has much longer wavelengths (note the different scales used on the figure axes).

library(microclimbook)
w<-c(0:3000)
r <- Planck(w) / (10^3 * 10^9)
par(mar=c(6, 6, 2, 2), mfrow = c(1, 2))
plot(r~w, type = "l", xlab = "Wavelength (nm)", ylab = 
    expression(paste("Spectral emittance (", kW/mm^3, ")")))
abline(v = 380, lwd = 2, lty = 2)
abline(v = 750, lwd = 2, lty = 2)
text(565, 15,"Visible", srt = 90)
text(900, 15,"Near-infrared", srt = 90)
w<-c(0:800)*100
r <- Planck(w, 15) / 10^6
w<- w / 1000 # convert to micro
plot(r~w, type = "l", xlab=expression(paste("Wavelength (",mu,m,")")), 
    ylab = expression(paste("Spectral emittance (", MW/m^3, ")")))

Emitted radiation

A second useful property of Planck’s Law is that it allows the total radiant energy emitted by a unit area of surface of a blackbody radiator to be calculated. This equation, referred to as the Stefan-Boltzmann equation and derived by the Slovene physicist Josef Stefan and the Austrian physicist Ludwig Eduard Boltzmann, essentially integrates the Planck equation over all wavelengths to give

\[\begin{equation} R_{em}=\sigma T^4 \end{equation}\]

where \(\small R_{em}\) is emitted radiation (W·m^-2) and \(\small \sigma\) is the Stefan-Boltzmann constant (5.6697 × 10^-8 W·m^-2·K^-4).

It is worth stating briefly that a blackbody is an idealized physical body that absorbs all incident electromagnetic radiation, regardless of frequency or angle of incidence. Both the sun and the earth behave nearly like a black body, but a general version of the Stefan-Boltzmann equation is given by

\[\begin{equation} \tag{2.1} R_{em}=\varepsilon \sigma T^4 \end{equation}\]

where \(\small \varepsilon\) is the emissivity of the object (with values in the range 0 to 1). The ground surface and vegetation typically have emissivity values of around 0.97.

In addition to emitted radiation, all objects absorb, reflect and/or transmit radiation and thus incoming radiation from the sun, for example, can always be portioned between these three components. It is thus worth describing these components in more detail. Before doing so, it is worth remarking briefly that from Kirchoff’s Law, emission and absorption of radiation are governed by the same process – that of changing the energy status of the emitted or absorbing atoms are molecules – and are thus equivalent. In general, however, when the term emissivity is used, it is usually to describe emissivity of longwave radiation. For opaque surfaces, owing to the conservation law, reflectivity in a given wavelength is one minus absorptivity. For translucent surfaces that are against a black background, reflectance, transmittance and absorption sum to one. The relevance of the black background becomes clear when one considers the distinction between reflectance and reflectivity.

Reflectance, reflectivity and albedo

Reflectance is the ratio of the radiant flux reflected from a material to the incident radiant flux and is therefore unitless. Reflectivity is the property at the boundary between two particular types of media, and for opaque surfaces the distinction is immaterial as the values are the same. However, when surfaces are translucent as is the case for vegetation for example, internal reflections and scattering become important and the two terms cannot be used interchangeably: reflectivity can be thought of as the very first intersection of the radiation with the surface. Over a black background, the various internal and external radiative fluxes must all balance and hence the conservation law applies.

Objects also reflect different amounts of radiation in different wavelengths. In the real world, humans experience this by perceiving the colour of objects. For example, leafs are generally green because chlorophyll is good at absorbing radiation in the visible spectrum except in green wavelengths – leaves have a fairly high reflectance of green light. They also have fairly high reflectance in near-infrared range, but the human eye cannot perceive this. In practical terms it is often convenient to average these wavelength specific effects across the full shortwave spectrum and the term albedo is used. The albedo of a surface is formally defined as the ratio of radiosity (the radiant flux density leaving a surface) to the irradiance (the radiant flux density received at the surface). Since radiative fluxes are measurements in units of energy per time interval, it is implicit that the spectral density (i.e. the amount of energy concentrated in a given wavelength) is accounted for. However, it is important to recognise that an albedo is not an intrinsic property of a surface. Instead, albedo depends on the spectral and angular distributions of the incident radiation, which in turn are influenced by atmospheric composition and the direction of the beam of radiation from the sun.

It is a consequence of the wavelength-dependent reflectance of most surfaces that albedo changes in response to the spectral distribution of solar radiation reaching the object of interest. Atmospheric constituents such as ozone, water and carbon dioxide absorb a significant amount of solar radiation in specific wavelengths, and in consequence, the spectral distribution of solar radiation changes in response to atmospheric conditions. Similarly, below canopy, the spectral distribution of solar radiation is quite different to that in open areas. However, for many practical purposes, it is adequate to define the albedo of a surface without explicitly considering how this changes in response to atmospheric absorption. The same is not necessarily true of anguler dependence.

Insight into the angular dependence of albedo can be gained by considering a few fundamental principles. When radiation is reflected at interfaces between media that are spatially uniform over scales large enough for individual surface elements to be treated as planar, it is reflected at an angle equal to the angle of incidence. However, many natural surfaces are composed of elements that have scales comparable to the wavelength of the incident light, and resultantly the reflected radiation is scattered in multiple directions. Such scattering can also occur when larger surface elements, such as leafs, are orientated differently with respect to the incident light. In consequence, light reflected by one surface element may be reflected by another. This is the reason why the calculation of radiative fluxes through vegetated canopies is not entirely straightforward: a topic to which we return later. It is also the reason that upward fluxes of reflected radiation are generally treated as diffuse (see below).

Because surface albedo depends on the angle of incidence of radiation (and hence also on the ratio of direct and diffuse radiation reaching the surface), it is common adopt a simplification in the representation of the surface albedo. The most popular approach is one in which the surface albedo is approximated by a simple weighting of two distinct surface albedo types, each associated with an extreme incident radiation field. The ‘black sky’ albedo, \(\small \alpha_{b}\), (more technically known as Directional Hemispherical Reflectance), is the albedo when radiation reaching the surface is purely collimated (i.e. entirely comprising direct radiation, which reaches the surface in one direction only). The ‘white sky’ albedo, \(\small \alpha_{d}\), (more technically known as Bi-Hemispherical Reflectance), is the albedo when radiation reaching the surface is completely isotropic (i.e. entirely comprising diffuse radiation, which when isotropic is equal in all directions). They can then be combined to approximate the surface ‘blue sky’ albedo \(\small \alpha\) as follows

\[\begin{equation} \tag{2.2} \alpha=F_d \alpha_d+F_b \alpha_b \end{equation}\]

Where \(\small F_d\) and \(\small F_b\) sum to 1 and correspond to the fractions of diffuse and direct (beam) to the total downward flux density of radiation respectively. We describe exactly what is meant by direct and diffuse radiation in the next section. In the section on radiative transfer within canopies, we provide a more detailed overview of the how leaf and ground reflectance contribute to overall canopy albedo.

Absorptance

A key feature of radiation absorption is that it depends on the angle of incidence of radiation relative to the surface of interest. It is therefore worth elaborating on the directional properties of radiation. Radiation absorbed by an object in the environment has several sources. The most obvious source of radiation is from the sun, but this can reach the canopy either as direct (beam) radiation or in the form of diffuse radiation, which is that scattered by particles and clouds in the atmosphere. Direct radiation is directional (referred to as anisotropic or collimated) and therefore its absorption depends on the angle of the surface relative to that of the solar beam.

Diffuse radiation is isotropic and therefore angle and orientation of the surface relative to the sun has no effect on the amount of diffuse radiation absorbed. Another source is solar radiation is that reflected from surrounding surfaces. In general this is anisotropic, but as reflected radiation in complex real environments may emanate from multiple sources, which themselves potential receive reflected radiation from multiple sources, it is often more convenient and computationally tractable to treat reflected radiation as isotropic. This negates the need to use a finite elements approach and perform complex simulations of photon scattering.

The final source is longwave (thermal) radiation comprising both that emitted downward from the sky and that emitted from surfaces surrounding the object. In practical terms, longwave radiation is usually treated as isotropic. However, an important caveat is needed to qualify the assumption of diffuse, reflected and longwave radiation being isotropic. They are generally only isotropic over a particular hemisphere that is in view. So, for example, a convex shaped animal on a bare ground surface receives diffuse radiation from above and the flux density of diffuse radiation received on the upper, illuminated surface of the animal is equal all over the upper surface. However, the underside of animal receives shortwave radiation that is reflected from the ground surface and this will not have the same flux density as the upper surface.

When a surface is exposed to direct radiation, the rays of light are nearly parallel and the flux density of radiation absorbed by the surface depends on its orientation with the respect to the sun. Intuitively one can grasp this concept by considering a torch beam shone oblique on a surface. The torch beam is far less concentrated than if it were shone directly. Where the term direct beam radiation is used, as in this book, it is taken to mean the flux density of radiation incident on a surface perpendicular to the solar beam.

The effects of slope angle are described mathematically by Lambert’s cosine law, which relates the flux density of direct radiation absorbed by an opaque surface (\(\small R_{abs}^b\)) to the flux density normal to the beam (\(\small R_0^b\)) as follows

\[\begin{equation} R_{abs}=(1-\alpha_b)(R_{0}^b \cos \theta)⁡ \end{equation}\]

where \(\small \theta\) is the angle of the radiant beam and normal to the surface. If the surface is horizontal, this is simply the solar zenith angle (the angle between the solar beam and vertical). If the surface itself is inclined, then the both the azimuth and zenith of the sun must be known (as well as the aspect of the surface). An extension to Lambert’s cosine law is thus

\[\begin{equation} \tag{2.3a} R_{abs}=(1-\alpha) R_0 \cos Z \cos s+\sin Z\sin s\cos(\Lambda_s-\Lambda) \end{equation}\]

where \(\small Z\) the solar zenith, \(\small s\) the slope angle of the surface measured from horizontal, \(\small \Lambda_s\) the solar azimuth (direction from north) and \(\small \Lambda\) the aspect of the surface (relative to north).

It is thus evident that to calculate radiation absorption in real environments one must know both the direct and diffuse flux density, but most meteorological datasets do not distinguish between the two, providing instead, an estimate of the downward flux of both direct and diffuse combined. Fortuitously, for a given location and solar zenith angle the diffuse fraction can usually be relative accurately estimated from the total shortwave flux, as both relate predictably to cloud cover. The function difprop, included with this book, implements the Skartveit et al (1998) method as in the example below.

rad <- c(1:1352) # typical values of radiation in W/m^2
jd <- julday(2022, 6, 21) # julian day
dfr <- difprop(rad, jd, 12, 49.97, -5.22) # Expected diffuse fraction and midday
plot(dfr ~ rad, type = "l", lwd = 2, 
     xlab = expression(paste("Incoming shortwave radiation (", W*M^-2, ")")),
     ylab = "Diffuse fraction")

It is useful to define a coefficient (\(s_c\)) given by

\[\begin{equation} \tag{2.3b} s_c=\cos Z \cos s+\sin Z\sin s\cos(\Lambda_s-\Lambda) \end{equation}\]

This is, in effect, the proportion of direct radiation intercepted by the surface relative to that received by a surface normal to the solar beam. Thus, \(\small s_c=1\) when the slope is normal to the solar beam and \(\small s_c\to 0\) as sunlight strikes the surface obliquely. The effects of slope and aspect on \(\small s_c\) in a real landscape can be seen pictorially by running the code below, which applies function solarcoef, which calculates \(\small s_c\). The solarcoef calculates the solar position automatically using the method detailed below.

require(plotly)
jd <- julday(2022, 6, 21)
sc <- solarcoef(dtm, 9, jd) # solar coefficient at 9am on 21st Jun 2022
# Interactive 3D plot
plot_ly(z = ~.isr(dtm)) %>%
        add_surface(surfacecolor = ~.isr(sc)) %>%
        layout(scene = list(zaxis = list(range = c(0, 200))))

Two further convenient metric are useful to define and calculate. Firstly, the ratio of the shadow area of an object on a surface perpendicular to the solar beam (\(\small A_p\)) relative to the total surface of the object (\(\small A\)). This is equivalent to the average flux density over the entire surface of the object relative to that normal to the surface. It is also equivalent to the silhouette area of the object normal to the solar beam (i.e. the area of the silhouette if the object and the sun were directly in line with one another) and would also be area of shadow cast by that object on a surface perpendicular to the solar beam. For example, a horizontal leaf with total two-sided area \(\small A\) would have a silhouette area of \(\small A_p⁄A=0.5\cos Z\) (0.5 because only one side of the leaf is illuminated). This is a particularly useful concept for more complex shapes such as spheres, cylinders, ellipsoids that approximate the shape of an animal, as it allows a mathematical calculation of the average flux density for objects where the angle of incidence to the solar beam varies across the object. For a variety of shapes, the silhouette area can be calculated using function silhouette. We return to this concept and demonstrate use of the function in the section on radiation absorption by animals.

The second useful metric to define is the area of shadow cast on a horizontal surface by an object, relative to the illuminated area of that surface (\(\small K\)). This is useful as it is equivalent to the ratio of radiation absorbed by a bumpy surface such as a vegetation canopy to that which would be absorbed by a flat surface with equivalent albedo (Monteith and Unsworth 2013). It also has other uses, for example, in the computation of transmission of radiation, to which we return below. The three metrics relate to each other as \(\small A_p/A=0.5s_c=0.5K/\cos Z\)

The sun’s position

To calculate radiation absorption in real environments, it is necessary to calculate the solar zenith and azimuth. The solar zenith (\(\small Z\)) depends on latitude \(\small \phi_y\) and time as follows

\[\begin{equation} \cos Z=\sin\delta\sin \phi_y+\cos\delta\cos\phi_y\cos h_s \end{equation}\]

where \(\small h_s\) is the hour angle in solar time (\(\small t_s\)) given by

\[\begin{equation} h_s=(\pi⁄12)(t_s-12) \end{equation}\]

and \(\small \delta\) is the current declination angle of the sun calculated from the Astronomical Julian day (\(\small j\)) as

\[\begin{equation} \delta=\frac{\pi 23.5}{180} \cos\lgroup 2\pi\frac{j-159.5}{365.25}\rgroup \end{equation}\]

The Astronomical Julian day is the continuous count of days since the beginning of the Julian period. The date 1^st Jan 2022 has an Astronomical Julian day of 2459581. The solar time is calculated from longitude (\(\small \phi_x\)) and local time (\(\small t_l\)) as

\[\begin{equation} t_s=t_l+\frac{4(\phi_x-\phi_m)-\Delta_t}{60} \end{equation}\]

where \(\small \phi_m\) is the longitude of the local time zone meridian (i.e. 0 for Greenwich Mean Time) and \(\small \Delta_t\) is the equation of time – a correction applied to account obliquity due to tilt of the Earth’s rotational axis and for the east-west component of the analemma, namely the angular offset of the Sun from its mean position on the celestial sphere due the eccentricity of the Earth’s orbit. These two factors have different wavelengths, amplitudes and phases that vary over geological timescales. From Milne (1921)

\[\begin{equation} \Delta_t=-7.659\sin(M+9.683)\sin(2M-3.5932) \end{equation}\]

where M is the mean anomaly given by \(\small 6.24004077+0.01720197(j-2451545)\).

The solar azimuth (\(\small \Lambda_s\)) is calculated as

\[\begin{equation} \sin\Lambda_s=-\frac{\sin h_s\cos\delta}{\sin Z} \end{equation}\]

A series of functions have been included with the book for calculating solar positions. The are julday for calculating \(\small j\), solartime for calculating \(\small t_s\), solarazi for calculating \(\small \Lambda_s\) and solarzen for calculating \(\small Z\). A demonstration of their application is provided in the example below

# Time series for 2022:
tme <- as.POSIXlt(c(0:364) * 3600 * 24, origin = "2022-01-01 00:00", tz = "UTC")
jd <- julday(tme$year + 1900, tme$mon + 1, tme$mday) # Julian day
Z<-solarzen(12, 49.97, -5.2, jd)  # Midday solar zenith at Caerthillian, UK
plot(Z ~ as.POSIXct(tme), type="l", xlab = "Month", ylab = "Midday solar zenith", lwd = 2)

Transmittance

Transmittance is relevant to consider when radiation is attenuated as it passes through turbid mediums such as air. It is also often convenient to consider canopies as translucent mediums, in which radiation is attenuated as one moves down through the canopy. Transmission of beam radiation has generally been described using an equation similar to Beer’s law. In its simplest form, consider an entirely translucent gas containing small and randomly placed particles that are shaped like horizontal flat plates with a total one-sided area of \(\small 0.5A\). If the sun is at its zenith, radiation passing through this gas will be attenuated by the presence of the particles, but since the particles are partially overlapping, the transmission of radiation is not directly and linearly related to the total surface area of particles, but rather tends to zero as the total surface area increases. The equation describing this relationship is

\[\begin{equation} \tau= \exp(-0.5A) \end{equation}\]

where \(\tau\) is the fraction of radiation transmitted through the medium - i.e. its transmittance. In this particular instance, the solar zenith angle makes no difference. Because the plates are horizontal, the chances of encountering particles does not change with zenith angle, but the same is not true if the particles are spherical as is more commonly the case. With spherical particles, a further relevance of \(\small K\) can be seen as it also describes the relationship between the total surface area of particles (\(\small A\)) and the transmission of radiation such that

\[\begin{equation} \tag{2.4} \tau=\exp(-KA) \end{equation}\]

For the sake of simplicity, let us instead define the sphere in terms of its diameter. In this instances, the shadow cast on the horizontal with the sun is at its zenith is equivalent to the diameter of the sphere and \(\small \tau_{(Z=0)}=exp(-d)\) where \(\small d\) is the summed diameter of all particles. A more general expression for any solar zenith angle is then given by

\[\begin{equation} \tau=\exp(-K_dd) \end{equation}\]

Where \(\small K_d=1/\cos Z\) which is in effect, the ratio of the path length through the medium relative to that when the sun is at its zenith. This concept is useful when considering attenuation of radiation through the atmosphere, though because of the curvature of the earth, calculation of the path length ratio is a little more complex. In the literature, this is usually referred to as the air mass coefficient. A further complication is that atmospheric constituents, notably water vapor, affect the attenuation of radiation through the atmosphere even in the absence of clouds. However, it is nevertheless useful to apply these principles to calculate the expected ‘clear-sky’ radiation for a given zenith angle as doing so provides a reasonable approximation of the diurnal cycle of direct radiation and provides a means to reproduce broadly realistic hourly values of radiation when, for example, only daily values are known. The function clearskyrad applies the Crawford & Duchon (1999) method for doing so and its application in the context of diurnal radiation cycles is demonstrated below. The function takes pressure, temperature and humidity as inputs, but is not greatly sensitive to these variables. Approximate values can therefore be provided, or the defaults chosen, if these data are unavailable.

tme<-as.POSIXlt(c(0:95) * 3600, origin = "2022-06-19", tz = "UTC")
jd <- julday(tme$year + 1900, tme$mon + 1, tme$mday) # julian days
radd <- c(180, 310, 180, 145) # average daily radiation over 4 days (W/m^2)
# hourly clear sky radiation at Caerthilliean: 
csr <- clearskyrad(tme$hour, 49.97, -5.22, jd) 
# Calculate daily clear sky radiation
dcsr <- matrix(csr, ncol = 24, byrow = TRUE) 
dcsr <- apply(dcsr, 1, mean)
# Calculate hourly clear sky fraction:
radf <- rep(radd / dcsr, each = 24)
# Calculate hourly radiation:
radh <- radf * csr
# Plot results:
plot(radh~as.POSIXct(tme), type = "l", lwd = 2, xlab = "Day", 
     ylab = expression(paste("Radiation (",W*m^-2,")")))

Direct radiative transfer through vegetated canopies

Though real canopies are structurally complex comprising numerous small surfaces with different angles and orientations, it is convenient to consider the way in which solar radiation is attenuated as it passes down through the canopy as being quite similar to the way in which it is attenuated by a turbid medium. For illustration, consider a canopy comprising randomly distributed, small, horizontal leaves that are completely opaque and black. Consider also a completely clear sky with the sun directly overhead (such that \(\small Z=0\) and there is no diffuse radiation). If one divides the canopy into multiple horizontal thin layers then the fraction of radiation intercepted by each layer is simply the total one-sided leaf area per unit ground area within each layer.

To describe this mathematically, let us envisage a series of nodes \(\small i\) situation at the top of each layer and let us number the nodes sequentially downward from the top of the canopy. Let us combine the effects of leafs and other canopy elements such as branches and tree tracks and let \(\small P_i\) be the cumulative area of these canopy constituents per unit ground area (from the top of the canopy) to node \(\small i\) and let \(\small P_{i+1}\) be the cumulative area at the node below. Let \(\small \Delta P=-(P_i-P_{i+1})\) (i.e. the canopy element area per unit ground area within a layer between \(\small i\) and \(\small i+1\) where the sign is negative as \(\small P_{i+1}>P_i\)). Let \(\small R_i^b\) and \(\small R_{i+1}^b\) be the flux densities of radiation at node \(\small i\) and \(\small i+1\) and let \(\small \Delta R^b=R_i^b-R_{i+1}^b\). It then follows that \(\small R_{i+1}^b=R_i^b-\Delta PR_i^b\), which as the layers become thinner such that \(\small \Delta\to d\), describes the rate of change in the flux density with canopy element area:

\[\begin{equation} \tag{2.5a} \frac{dR^b}{dP}=-R_i^b \end{equation}\]

The notation used in the equation above and the description of how it is derived, is not the most intuitive. It is presented in this way primarily because it is similar to the format used to capture the more complex two-stream radiative transfer model describe below. It is therefore worth elaborating on its derivation with a simple example. However, all that really needs to be understood from this is that the fraction of radiation intercepted relative the layer above is determine by the leaf (and other canopy element) area within each layer, and hence, radiation is attenuated exponentially through the canopy such that

\[\begin{equation} \tag{2.6a} R_i^b=R_0^b\exp(-P_i) \end{equation}\]

where \(\small R_0^b\) is the flux density of direct radiation at the top of the canopy. As with the turbid medium where small particles are responsible for the radiation attenuation, the reason for this exponential decay is due to the the horizontal overlap of canopy elements owing to their random distribution – i.e. some of the radiation in any given layer has already been intercepted by canopy elements in layers above.

The same concepts can now be extended to consider situations where the canopy elements are not horizontal. Consider a single leaf at an inclined angle – the flux density of direct radiation intercepted by that leaf per unit one sided leaf area is the same as the area of shadow cast on the horizontal per unit illuminated leaf area, which has already previously been defined as \(\small K\). If the leaf is at a \(\small 45^\circ\) degree angle relative to the solar beam, for example, then \(\small K=\cos(45^\circ)\approx 0.707\). Equations (2.5a) and 2.6a) can then be written in a more general form as

Thus, when the sun is close to its zenith, planophile canopy elements will cast a greater shadow area and intercept more radiation than erectophiles with equivalent area. Conversely, when the sun is low above the horizon, erectophiles will intercept more radiation, and in real canopies, the attenuation of radiation is thus dependent on the zenith angle.

This is all relatively straightforward for a single elaf with a known inclination angle, but in real canopies a complexity is that constituents of in the canopy are at different angles. Thus, rather than defining each canopy element individually, it is more convenient to characterise the canopy as comprising surfaces that have a continuous distribution of inclinations. From Campbell (1986) and Campbell (1990), it is reasonable to assume that the real distribution of inclinations can be approximated by assuming they conform to a prolate or oblate spheroid distribution. By adjusting the ratio of horizontal to vertical axis of the spheroid, canopy element angle distributions of any canopy from erectophile to planophile can be simulated. Defining \(\small x\) as the ratio of average projected elements on horizontal surfaces, such that \(\small x=0\) for a vertical distribution, \(\small x=\infty\) for a horizontal distribution and \(\small x=1\) for a spherical distribution, \(\small K\) is approximated as

\[\begin{equation} \tag{2.7} K\approx\frac{\sqrt{x^2+\tan^2Z}}{x+1.774(x+1.182)^-0.733} \end{equation}\]

Typical values of \(\small x\) for different habitat types are shown in Table 2.1. In the special case where \(\small x=1\), the equation can be simplified to \(\small 1/(2\cos Z)\) and when \(\small x=0\), \(\small K=2/(\pi\tan(0.5\pi-Z))\). When \(\small x\to\infty\), \(\small K\to1\) and therefore does not depend on the sun’s position. \(\small K\) is computed using the function Kcanopy. In the example below, variation in \(\small K\) with \(\small Z\) is shown for various values of \(\small x\).

Z<-c(0:90)
K1<- Kcanopy(Z, 0)
K2<- Kcanopy(Z, 0.5)
K3<- Kcanopy(Z, 1)
K4<- Kcanopy(Z, 3)
K5<- Kcanopy(Z, Inf)
par(mar=c(6,6,2,2))
plot(K1 ~ Z, type = "l", xlab = "Solar zenith angle (degrees)", ylab = "Extinction coefficient (K)", ylim = c(0, 2))
par(new=T)
plot(K2 ~ Z, type = "l", xlab = "", ylab = "", ylim = c(0, 2))
par(new=T)
plot(K3 ~ Z, type = "l", xlab = "", ylab = "", ylim = c(0, 2))
par(new=T)
plot(K4 ~ Z, type = "l", xlab = "", ylab = "", ylim = c(0, 2))
par(new=T)
plot(K5 ~ Z, type = "l", xlab = "", ylab = "", ylim = c(0, 2))
text(5, 0.16, "x = 0.0")
text(5, 0.36, "x = 0.5")
text(5, 0.56, "x = 1.0")
text(5, 0.88, "x = 3.0")
text(8, 1.05, "horizontal")

The same principles can be extended to diffuse radiation, but here we must account for the fact that canopy elements are not black and entirely opaque and thus scatter radiation both upwards and downwards. We must also account for the upward scattering of radiation by the ground surface. Making the realistic assumption that the upward and downward scattering of radiation are diffuse (and in consequence that there is no upward direct radiation), equation 2b can be used to calculate direct radiative transfer irrespective of leaf reflectance and transmittance. However, a rather more convoluted approach must be used for diffuse radiation.

Diffuse radiatiave transfer: two-stream models

To account fully for the scattering of radiation one must consider the inclinations and position of individual canopy elements. Several such models that are rigorous and realistic in the way that they handle interception, scattering and absorption of radiation have been designed and tested (see e.g. Goudriaan 1977 and Gastellu-Etchegorry 2004). These models use a finite-element approach to describe the process of radiative transfer, whereby the rays of radiation and their interaction with individual canopy elements is traced numerically. Although these techniques allow for the inclusion of intricate details of the optical and structural characteristics of canopies they tend to be computationally expensive and cumbersome to use.

In many circumstances it is helpful to assume that a canopy is horizontally homogeneous at least over extents relevant for quantifying radiative transfer at any given location. It then becomes possible to solve the radiative transfer process using a two-stream approximation – i.e. by considering two discrete directions only: upward and downward. Various models of this type exist and here the similarities and differences among four broadly similar models are described. These are (1) the Ross (1981) model, (2) the Verhoef (1984) SAIL model; (3) the Pinty et al. (2006) model and (4) the model first proposed by Dickinson (1983), but given a full solution by Sellers (1985). This ‘Dickinson-Sellers’ model has been widely adopted in the land surface component of climate models. A more comprehensive description of these models and overview of their differences is given in Yuan et al. (2017), upon which the following draws heavily.

Some of elements of this and the following sections are quite mathematically involved. The salient points for readers who are more interested in the application of the modelling approaches are threefold. Firstly, downward diffuse radiation is attenuated as it passes through the the canopy but some of it is scattered in various directions. The scattering effects depend on properties of the canopy, but can be handled by quantifying the contribution to a downward and upward stream of radiation. Secondly, upward radiation can generally be considered as entirely diffuse, but both the downward diffuse and direct radiation contribute to the upward stream. Thirdly, the mathematical heavy lifting is done in function twostream and examples of how to use this function are provided.

With a bit of rearrangement, the diffuse radiation component of the four models can be written using two sets of differential equations in a form following that of equation (2.5a)

\[\begin{equation} \tag{2.8a} -\frac{dR_i^{d\uparrow}}{dP}=-(a+\gamma)R_i^{d\downarrow}+sR_i^{b\downarrow} \end{equation}\] \[\begin{equation} \tag{2.8b} \frac{dR_i^{d\downarrow}}{dP}=-(a+\gamma)R_i^{d\downarrow}+\gamma R_i^{d\uparrow}+s^,R_i^{b\downarrow} \end{equation}\]

Here \(\small R_i^{d\uparrow}\) is upward diffuse radiation (i.e. that reflected upward by canopy elements and the ground surface, \(\small R_i^{d\downarrow}\) is downward diffuse radiation (analogous to the processes described for direct radiation above) and \(\small R_i^{b\downarrow}\) is downward direct radiation quantified as in equation (2.6b).

Each of the equation’s components describes a different physical process. In equations (2.8a) and (2.8b), the left hand term describes the overall attenuation of upward or downward radiation through the canopy. The first right hand term in (2.8a) describes the fraction of upward diffuse radiation that is re-scattered in upward direction and the equivalent in (2.8b) describes the fraction of downward diffuse radiation that is re-scattered downward. This is the dominant process driving radiation attenuation and is broadly analogous to the processes described for direct radiation earlier. \(\small a\) is the absorption coefficient for incoming diffuse radiation per unit leaf area and \(\small \gamma\) is the backward scattering coefficient. In essence, when canopy elements are darker and there is more back-scattering, more of the radiation is intercepted. The second term on the right hand side in (2.8a) describes the fraction of the downward diffuse radiation that is converted into an upward diffuse flux by back-scattering and the equivalent in (2.8b) describes the fraction of the upward diffuse radiation that is converted into an downward diffuse flux by back-scattering. Both terms naturally increase with \(\small \gamma\). The final term on the right hand side of (2.8a) refers to the contribution to the upward diffuse flux by the scattering of direct radiation penetrating into the canopy to \(\small i\). The equivalent in (2.8b) is the contribution from scattering of direct radiation to the downward diffuse flux. \(\small s\) and \(\small s^,\) are thus the backward and forward scattering coefficient for direct radiation respectively.

In all four models, \(\small a\) is given by

\[\begin{equation} a=(1-\omega)K^* \end{equation}\]

where \(\small K^*\) is the extinction coefficient for diffuse radiation (equivalent to \(\small K\), but for diffuse radiation) and \(\small \omega\) is the single scattering albedo of individual canopy elements given by \(\small \omega=\alpha_P+\tau_P\) where \(\small \alpha_P\) is canopy element reflectance and \(\small \tau_P\) is canopy element transmittance (canopy elements such as leafs are partially translucent).

The different models calculate \(\small K^*\) in different ways. Under strictly isotropic conditions, i.e. those where flux densities are equal in all directions, \(\small K^*=1\) regardless of the angles of inclination of canopy foliage. The diffuse radiation within the canopy is not exactly isotropic once it has interacted with canopy and been further scattered. In the Ross and SAIL models, diffuse radiation is assumed strictly isotropic (i.e. \(\small K^*=1\)). In the Dickinson-Sellers and Pinty model, diffuse radiation is assumed to have a directional bias once it has interacted with canopy and been further scattered as given by \(\small K^*=[\int^0_1\mu/G(\mu)]^{-1}\) where \(\small \mu=\cos Z\) and \(\small G(\mu)\) is the mean projected area of leaves in the direction \(\small Z\). However, a consequence of calculating \(\small K^*\) in this way is that the calculations of \(\small s\) and \(\small s^,\) become inconsistent with \(\small K^*\) and, following Yuan et al. (2017), \(\small K^*\) is assumed in this book to equal 1.

The models also differ slightly in the way that \(\small \gamma\) is calculated, but can all be identically expressed as

\[\begin{equation} \gamma=0.5(\omega+J\delta)K^* \end{equation}\]

Here, \(\small \delta\) is simply reflectance less transmittance and is given by \(\small \delta=\alpha_P-\tau_P\) and \(\small J\) is an integral function of the inclination distribution of canopy elements. The Ross model makes a simplification that assumes the special case where \(\small \alpha_P=\tau_P=\omega/2\), which means that \(\small J=0\). The Dickinson-Sellers model, instead of considering the inclination density distribution of canopy elements explicitly, derives its effects from the mean inclination angle. \(\small J\) is thus given by \(\small J=\cos \hat{\Theta}_P\) where \(\small \hat{\Theta}_P\) is the mean inclination angle of canopy elements in the zenith direction. Campbell (1990) gives a method for deriving \(\small \hat{\Theta}_P\) (in Radians) from the \(\small x\) coefficient described previously as \(\small \hat{\Theta}_P=9.65(3+x)^{1.65}\). An alternative derivation is given in Wang & Jarvis (1988).

The Pinty and SAIL models give the most explicit formulation for backward scattering and in both models \(\small J\) is given by

\[\begin{equation} J=\int^{\pi/2}_0f(\Theta_P)\cos^2\Theta_Pd\Theta_P \end{equation}\]

where \(\small f(\Theta_P)\) is the inclination density distribution function of canopy elements and \(\small \Theta_P\) is the inclination angle of canopy elements in the zenith direction. Campbell (1990) gives a general analytical expression for \(\small f(\Theta_P)\) as

\[\begin{equation} f(\Theta_P)=\frac{2x^3\sin\Theta_P}{\Lambda_P(\cos^2\Theta_P+x\sin^2\Theta_P)^2} \end{equation}\]

where \(\small \Lambda_P\) is the normalised ellipse area given by \(\small \Lambda_P=x+1.774(x+1.182)^{-0.733}\) (cf equation 2.7). Then, by performing the integration

\[\begin{equation} J=\frac{x^2\arctan\lgroup\frac{\sqrt{1-x^2}}{x}\rgroup}{\Lambda_P(1=x^2)^\frac{3}{2}}-\frac{x^3}{\Lambda_P\lgroup(1-x^2)+x^2\rgroup} \end{equation}\]

when \(\small x<1\)

\[\begin{equation} J=\frac{1}{3} \end{equation}\]

when \(\small x=1\) and

\[\begin{equation} J=y^2\sqrt{1-y^2}\ln(\frac{\sqrt{1-y^2}-y^2+1}{\sqrt{1-y^2}+y^2-1})+2y^2-2 \end{equation}\]

when \(\small x>1\).

In the example below, the function scattercoef is used to calculate \(\small \gamma\) as a function of \(\small \alpha_P\) and a visual comparison of the three methods made. As \(\small \alpha_P\rightarrow0\), \(\small \gamma\rightarrow 0\) but increases when leaf reflectance increases, mostly rapidly when canopy elements are more horizontally orientated. The Dickenson-Sellers and the identical Pinty and SAIL methods give very similar results, but large differences are observed with the Ross method for planophiles since this method is assumed not to depend on leaf angle distributions. Given the rather tedious integration involved in deriving \(\small J\) using the Pinty and SAIL methods, the far simpler Dickenson-Sellers method is probably preferential for most applications.

leafr <- c(0:800) / 1000
leaftr <- 0.2 * leafr
gma1 <- scattercoef(leafr, leaftr, 0.2, method = "Pinty")
gma2 <- scattercoef(leafr, leaftr, 0.2, method = "Dickenson-Sellers")
gma3 <- scattercoef(leafr, leaftr, 0.2, method = "Ross")
plot(gma1 ~ leafr, type = "l", xlim = c(0,0.8), ylim = c(0,0.8), lwd = 2, 
     xlab = "Leaf reflectance", ylab = "Scatter coefficient")
par(new=T)
plot(gma2 ~ leafr, type = "l", xlim = c(0,0.8), ylim = c(0,0.8), lwd = 2,
     xlab = "", ylab = "", lty = 2)
par(new=T)
plot(gma3 ~ leafr, type = "l", xlim = c(0,0.8), ylim = c(0,0.8), lwd = 2,
     xlab = "", ylab = "", lty = 3)
text(0.4, 0.8, "x = 0.2")
gma1 <- scattercoef(leafr, leaftr, 5, method = "Pinty")
gma2 <- scattercoef(leafr, leaftr, 5, method = "Dickenson-Sellers")
gma3 <- scattercoef(leafr, leaftr, 5, method = "Ross")
plot(gma1 ~ leafr, type = "l", xlim = c(0,0.8), ylim = c(0,0.8), lwd = 2, 
     xlab = "Leaf reflectance", ylab = "Scatter coefficient")
par(new=T)
plot(gma2 ~ leafr, type = "l", xlim = c(0,0.8), ylim = c(0,0.8), lwd = 2,
     xlab = "", ylab = "", lty = 2)
par(new=T)
plot(gma3 ~ leafr, type = "l", xlim = c(0,0.8), ylim = c(0,0.8), lwd = 2,
     xlab = "", ylab = "", lty = 3)
text(0.4, 0.8, "x = 5")

In most of the models, the calculation of the backward scattering coefficient for direct radiation is similar to that of \(\small \gamma\). In the Ross model \(\small s=0.5\omega K\) and in the Pinty and SAIL models

\[\begin{equation} s=0.5(\omega+\frac{J\delta}{K})K \end{equation}\]

The largest difference is observed with the Dickinson-Seller model in which single scattering albedo is used to derive \(\small s\). However, deriving the single scattering albedo is not straightforward and solution offered in Sellers (1985) underestimates \(\small s\) in comparison to other methods particularly when leaf transmittance is low (Yuan et al. 2017). In this book, the Pinty and SAIL method is therefore adopted. In all four models, the forward scattering coefficient is given by

\[\begin{equation} s^,=\omega K-s \end{equation}\]

Solving the differential equations for a two-stream model is not entirely straightforward, but is achieved by setting boundary conditions at the bottom and the top of the canopy. Each of the four models uses a rather different approach and full explanation is given in Yuan et al (2017). Here the solution given by Sellers (1985) is presented, but with equations rearranged to match the terms described above. Using the Sellers method, the boundary conditions are set such that contribution of direct radiation to the diffuse flux is zero at the top of the canopy and determined by soil reflectance (\(\small \alpha_G\)) and the one-sided area of canopy elements per unit ground area of the entire canopy (\(\small P_{AI}\)). Consequently, the upward diffuse flux at the ground surface (\(\small R_G^{d\uparrow}\)) is given by \(\small \alpha_G\lgroup R_G^{d\downarrow}+\exp(-KP_{AI})\rgroup\), where \(\small R_G^{d\downarrow}\) is the downward diffuse fluxes at the ground surface. Assuming \(\small K^*=1\), the upward and downward diffuse radiation fluxes are then given by

\[\begin{equation} \tag{2.9a} R_i^{d\uparrow}=R_0^{d\downarrow}[p_1\exp(-hP_i)+p_2\exp(hP_i)]+R_b^{d\uparrow} \end{equation}\]

\[\begin{equation} \tag{2.9b} R_i^{d\downarrow}=R_0^{d\downarrow}[p_3\exp(-hP_i)+p_4\exp(hP_i)]+R_b^{d\downarrow} \end{equation}\]

where Where \(\small R_b^{d\uparrow}\) and \(\small R_b^{d\downarrow}\) are the contribution of direct radiation to the diffuse upward and downward flux respectively given by

\[\begin{equation} \tag{2.9c} R_b^{d\uparrow}=R_0^{b\downarrow}[\frac{p_5}{\sigma_p}\exp(-KP_i)+p_6\exp(-hP_i)+p_7\exp(hP_i)] \end{equation}\]

\[\begin{equation} \tag{2.9d} R_b^{d\downarrow}=R_0^{b\downarrow}[\frac{p_8}{\sigma_p}\exp(-KP_i)+p_9\exp(-hP_i)+p_{10}\exp(hP_i)] \end{equation}\]

The derivation of the terms \(p_{1..10}\) are as follows:

The function twostream implements the two stream models described above, giving options for different methods of calculation. Using this function, the example below shows how albedo can be calculated and how the flux density of upward and downward radiation changes through the canopy.

albs <- twostreamrad(Rb0 = 100, Rd0 = 400, Z = 45, PAI = 0, PAIt = 3, leafr = 0.5, leaftr = 0.2, x = 1, groundr = 0.15, out = "albedo")
albs
PAI<-c(3000:0)/1000
z<-c(0:3000) / 200
rads <- twostreamrad(Rb0 = 100, Rd0 = 400, Z = 45, PAI = PAI, PAIt = 3, leafr = 0.5, leaftr = 0.2, x = 1, groundr = 0.15, out = "radiation")
plot(z~rads$Rddown, type = "l", xlim = c(0, 420), lwd = 2, ylab = "Height (m)",
     xlab = expression(paste("Radiation (", W/m^2, ")")))
par(new = T)
plot(z~rads$Rbdown, type = "l", xlim = c(0, 420), lwd = 2, lty = 2, ylab = "", 
     xlab = "")
par(new = T)
plot(z~rads$Rup, type = "l", xlim = c(0, 420), lwd = 2, lty = 3, ylab = "", 
     xlab = "")
text(300, 10.6, "Downward diffuse")
text(38, 14.5, "Downward direct")
text(123, 11, "Upward")

At the top of the canopy \(\small P_i=0\), and hence the exponential terms in equation 2.9 become 1. Since albedo is defined as the ratio of upward to downward radiation, the terms \(\small p_{1..10}\) also relate ground and canopy element reflectance and transmittance to the ‘white sky’ and ‘black sky’ albedos of the Vegetated surface as follows

\[\begin{equation} \tag{2.10a} \alpha_d=p_1+p_2 \end{equation}\] \[\begin{equation} \tag{2.10b} \alpha_b=\frac{p_5}{\sigma_p}+p_6+p_7 \end{equation}\]

Technically the albedos relate to the combined effect of the canopy and the underlying ground surface, but as we have alluded to in the previous chapter and explain in more detail below, this is convenient.

Both the ‘white sky’ and ‘black sky’ albedos can thus be related to the reflectance of the canopy elements and the ground surface. While ‘white sky’ albedo varies somewhat through time because cloud cover and atmospheric water vapour filter radiation in specific wavelengths, these effects are fairly minor except seasonally, when the decidious nature of many canopies comes into play. By contrast, ‘black sky’ and hence also ‘blue sky’ albedo varies through time according to the solar zenith angle.

In a practical sense it would be convenient to apply the calculation in reverse to derive leaf reflectance. This is because albedo can often be derived spatially across regions from remotely-sensed imagery, whereas the measurement of canopy element and ground reflectance across entire regions would require a laborious number of field measurements Unfortunately, there is no unique solution for the reverse calculation. However, by specifying a fixed leaf transmittance or by defining it in terms of a ratio of leaf transmittance to reflectance, it is possible to derive solution numerically by iteration. In one of the lengthier examples at the end of the chapter we show how how this can be done using the function leafrfromalb, which applies the uniroot function in R to do so. It is moderately computationally efficient.

Heterogeneous canopies and clumping

Thus far we he have considered canopies in which elements such as leafs are small and randomly distributed such that they behave like a turbid medium. In reality, canopy elements are non-randomly distributed and radiation can penetrate through larger gaps in he canopy unimpeded. The easiest way to handle this is to define a gap fraction (\(\small F_G\)) representing the fraction of the canopy that is unobscured by vegetation when the sun is at its zenith. For the simple case when the sun is at its zenith, the transmission of direct beam radiation is described as

\[R_i^b=R_0^b\exp(-K\tilde{P_i})+F_G\] where \(\small\tilde{P_i}\) is the cumulative canopy element area now concentrated into ‘non-gaps’ given by \(\small\tilde{P_i}=P_i/(1-F_G)\). The probability of the sunbeam’s path through the canopy being unobscured diminishes as the optical path length increases, but the extent to which it does depends on the shape of the vegetation clumps. Where the clumps are spherical or equal in the vertical and horizontal direction, the optical path length (\(\small K_c\)) relative to that when the sun is at its zenith is simply given by \(\small 1/\cos Z\). However, it is perhaps more reasonable to assume that the shape of the clumps is dictated by the inclination angles of leaves such that canopies with planophile leafs will typically have planophile clumps of vegetation. Under these circumstances \(\small K_c = K(Z) / K(0)\) where \(\small K(Z)\) and \(\small K(0)\) are the values of \(\small K\) derived using equation 2.7, when the solar zenith angles are \(\small Z\) and \(\small 0\) respectively. The transmission of beam direct beam radiation for any given zenith angle is then described as follows

\[\begin{equation} \tag{2.4c} R_i^b=R_0^b\exp(-K\tilde{P_i})+F_G^{K_c} \end{equation}\]

In the example code below, the transmission fraction as a function of \(\small K\) is shown for three values of \(\small F_G\)

Z <- c(0:80)
K <- 1 / (2 * cos(Z * pi / 180))
PAI <- 3
tr <- exp(-K * PAI / (1 - 0.5)) + 0.5^(K / 0.5)
plot(tr ~ K, type = "l", lwd = 2, xlim = c(0.5, 3), ylim = c(0, 0.6), xlab = "K", ylab = "Radiation transmission")
text(0.92, 0.4, expression(paste(F[G]," = 0.5")))
par(new = TRUE)
tr <- exp(-K * PAI / (1 - 0.3)) + 0.3^(K / 0.5)
plot(tr ~ K, type = "l", lwd = 2, xlim = c(0.5, 3), ylim = c(0, 0.6), xlab = "", ylab = "")
text(1.27,0.1,expression(paste(F[G]," = 0.3")))
par(new = TRUE)
tr <- exp(-K * PAI / (1 - 0)) + 0^(K / 0.5)
plot(tr ~ K, type = "l", lwd = 2, xlim = c(0.5, 3), ylim = c(0, 0.6), xlab = "", ylab = "")
text(0.75,0.05,expression(paste(F[G]," = 0")))

For diffuse radiation, \(\small K^*=1\) irrespective of the shape of the path length, but since the gap fraction is defined in terms of radiation transmission when the sun is at its zenith, the following modifications can be made to equations 2.9a-d:

\[\begin{equation} \tag{2.9e} R_i^{d\uparrow}=R_0^{d\downarrow}[(1-\tau_d)(p_1\exp(-h\tilde{P_i})+p_2\exp(h\tilde{P_i}))+\alpha_G\tau_d]+R_b^{d\uparrow} \end{equation}\]

\[\begin{equation} \tag{2.9f} R_i^{d\downarrow}=R_0^{d\downarrow}[(1-\tau_d)(p_3\exp(-h\tilde{P_i})+p_4\exp(h\tilde{P_i}))+\tau_d]+R_b^{d\downarrow} \end{equation}\]

\[\begin{equation} \tag{2.9g} R_b^{d\uparrow}=R_0^{b\downarrow}[(1-\tau_b)(\frac{p_5}{\sigma_p}\exp(-K\tilde{P_i})+p_6\exp(-h\tilde{P_i})+p_7\exp(h\tilde{P_i}))+\alpha_G\tau_b] \end{equation}\]

\[\begin{equation} \tag{2.9h} R_b^{d\downarrow}=R_0^{b\downarrow}[(1=\tau_b)(\frac{p_8}{\sigma_p}\exp(-K\tilde{P_i})+p_9\exp(-h\tilde{P_i})+p_{10}\exp(h\tilde{P_i}))+\tau_b] \end{equation}\]

where \(\small \tau_d\) and \(\small \tau_b\) are the transmission of diffuse and direct beam radiation through larger gaps in the canopy given by:

\[\tau_d=F_G^{\frac{1}{K(0)}}\] and

\[\tau_b=F_G^{K_c}\]

and the equations for calculating albedo are modified to

\[\begin{equation} \tag{2.10c} \alpha_d=(1-\tau_d)(p_1+p_2)+\tau_d \end{equation}\] \[\begin{equation} \tag{2.10d} \alpha_b=(1-\tau_b)(\frac{p_5}{\sigma_p}+p_6+p_7)+\tau_b \end{equation}\]

In one of the examples below, we show how the coefficients \(\small x\) and \(\small F_G\) can be derived from a hemispherical photo of a canopy.

Canopies over inclined ground surfaces

For direct radiation, a correction must also be applied when the canopy overlies an inclined ground surface. While it can be assumed, for a given habitat type, that the inclination of canopy elements is unaffected by ground surface inclination, the optical path length through the canopy will be shorter for ground surfaces facing in the direction of the sun. A full equation for the transmission of direct radiation through the canopy is thus given by

\[\begin{equation} \tag{2.4d} R_i^b=R_0^b[(1-F_G)\exp(-\tilde{K}\tilde{P_i})+F_G] \end{equation}\]

where \(\small \tilde K=K\cos Z/s_c\) and \(\small s_c\) is given by equation 2.3b.

In the two-stream model, the effects of an inclined gorund surface can be handled by adjusting the ground reflectance when calculating the contribution of direct beam radiation to the upward diffuse beam such that \(\small \tilde{\alpha_G}=\alpha_G\cos Z/s_c\) where \(\small \tilde{\alpha_G}\) is the adjusted ground reflectance.

Longwave and emitted radiation

Before turning to describing the transfer of longwave radiation within vegetated canopies it is worth commenting briefly on four variables that are commonly available as gridded climate products. The first is the upward radiation flux (\(\small R_0^{l\uparrow}\)), which is related to surface temperature (\(\small T_s\)) as \(\small R_0^{l\uparrow}=\varepsilon_s\sigma T_s^4\), where \(\small \varepsilon_s\) is the emissivity of the surface. The second is the downward longwave radiation flux (\(\small R_0^{l\downarrow}\)). This is simply the radiation emitted downward from the sky from clouds and other particles in the atmosphere. The third is the effective sky temperature (\(\small T_{sky}\)). This is the temperature the sky would need to be to emit downward radiation assuming it to behave like a black body such that \(\small R_0^{l\downarrow}=\sigma T_{sky}^4\). The third is effective sky emissivity (\(\small \varepsilon_{sky}\)), which describes the relationship between upward and downward radiation as \(\small R_0^{l\downarrow}=\varepsilon_{sky}R_0^{l\uparrow}\).

Leafs and other elements of the canopy typically have very low reflectance to longwave radiation so there is limited scattering of longwave radiation. In consequence, any longwave radiation within the canopy can be thought of as originating from three sources. Firstly from the ground surface, where incident ongwave radiation at any given height below canopy is determined by ground temperature and the tranmission of radiation from the ground surface. The second is from all the individual canopy elements, which in turn is determined by the temperature of those surfaces and tranmission from those surfaces. The third is downward from the sky, which is determined by the effective sky temperature and tranmission from the sky. It is convenient here also to consider the flux density per unit two-sided area - i.e. that incident on both the upward-facing and downward-facing surface of the canopy element, but scaling these values such that the total is given in terms of two-sided as opposed to one-sided flux densities. In essence, this means dividing the one-sided flux densities by 2.

We first provide a more complex and explicit way in which longwave radiation within a canopy can be calculated. Next we provide a far simpler method in which the radiation is derived from mean canopy temperature.

Dividing the canopy into \(\small n\) layers and defining \(\small \Delta P_i\) as the total one sided plant area within a layer \(\small i\), the total (two-sided) flux density of longwave radiation (\(\small R_j^l\)) incident on a canopy element \(\small j\) below canopy is given by

\[\begin{equation} \tag{2.11a} R_j^l=0.5(\tau_{0,j}\varepsilon_G\sigma T_G^4 + \tau_{j,h}\sigma T_{sky}^4+\sum_{0}^{i=n}\Delta P_i \tau_{i,j}\varepsilon_i\sigma T_i^4) \end{equation}\]

where \(\small \tau_{0,j}\) is transmission from the ground surface to \(\small j\) given by \(\small \exp(-P_{0,j})\) where \(\small P_{0,j}\) is the total one sided plant area between \(\small 0\) and \(\small j\), \(\small \tau_{j,h}\) is transmission from the top of the canopy at height \(\small h\) to \(\small j\) given by \(\small \exp(-P_{j,h})\) where \(\small P_{j,h}\) is the total one sided plant area between \(\small j\) and \(\small h\) and \(\small \tau_{i,j}\) is transmission between \(\small i\) and \(\small j\) given by \(\small \exp(-P_{i,j})\) where \(\small P_{i,j}\) is the total one sided plant area between \(\small i\) and \(\small j\). \(\small \varepsilon_G\) and \(\small \varepsilon_i\) are the emissivitities of the ground surface and canopy elements respectively and \(\small T_G\), \(\small T_{sky}\) and \(\small T_i\) are absolute temperatures of the ground surface, sky and canopy elements respectively.

The issue, however, is that in order to calculate the longwave radiation incident on a canopy element one must know the vertical temperature profile of canopy elements, but in order to calculate this, one must know the longwave radiation absorbed by that element. While the problem can be tackled iteratively (see Chapter 5), it is also convenient to a simplification. If one assumes that the range of canopy temperatures is in the order of a few degrees and are sufficiently absolute zero, then the mean radiation emmitted by canopy elements \(\small \overline{R_{em}}\) is approximately equal to \(\small \varepsilon_i\sigma\overline{T_i}^4\) where \(\small \overline{T_i}\) is the mean temperature of the canopy. If one assumes further that ground temperature is not unduly different from the temperature of canopy elements then

\[\begin{equation} \tag{2.11b} R_j^l\approx 0.5\tau_{j,h}\sigma T_{sky}^4+(1-0.5\tau_{j,h})\overline{\varepsilon_i}\sigma\ T_c^4 \end{equation}\]

where \(\small \overline{\varepsilon}\) is the combined mean emissivity of canopy elements and the ground surface and \(\small T_c\) is the temperature of the heat exchange surface of the canopy calculated using the method detailed in Chapter 1. In effect, equation 2.5 assumes longwave radiation incident on a surface can be thought of as originating from two sources: that downward from the sky (determined by both effective sky temperature and the transmisison of radiation downward through the canopy) and from the ground and surrounding canopy (determined by the mean temperature of canopy elements and the ground and by the portion of the view sphere surrounding the canopy element that is not sky). Fig 2.9 compares the two approaches for effective sky emissivities of 0.7 and 0.95 corresponding to clear and cloudy days using the example code presented in the last section of this chapter. In Chapter 5 we explore these issues in greater detail, demonstrating that though \(\small T_c\) is approximately equal to the mean canopy and ground temperature it is, in fact, a useful but somewhat arbitrary construct and its relationship with ground and canopy temperatures is time-variant.

Absorption of radiation by the ground and canopy combined

In Chapter 1 we describe how the concept of a heat exchange surface for a canopy is useful when calculating air temperature above canopy. When calcuating the energy budget of this surface, the combined effects of both the ground and canopy must be considered, and the absorbed radiation component of this energy budget is, in effect, the radiation absorbed by both the canopy and the ground. Since the ground is opaque, the absorbed radiation is given simply by

\[R_{abs}^{CG}=R_0^{d\downarrow}+R_0^{b\downarrow}-R_0^{d\uparrow}+\varepsilon_{CG}R_0^l\]

where \(\small R_{abs}^{CG}\) is radiation absorbed by the ground and canopy, \(\small R_0^{d\uparrow\) is the upward shortwave radiation flux at the top of the canopy, \(\small\varepsilon_{GC}\) is combined ground and canopy emissivity and \(\small R_0^l\) is downward longwave radiation from the sky given by \(\small\sigma T_{sky}^4\). Since the ‘blue-sky’ albedo surface is ratio of the upward to downward (direct and diffuse) flux, it follows that

\[\begin{equation} \tag{2.12} R_{abs}^{CG}=(1-\alpha)(R_0^{d \downarrow}+R_0^{b\downarrow})+\varepsilon_{CG}R_0^l \end{equation}\]

where \(\alpha\) is the ‘blue-sky’ albedo given by equations 2.2 and 2.10.

Absorption of radiation by the ground

In chapter 5, we show why it is necessary to calculate radiation absorbed by the ground surface alone (\(\small R_{abs}^G\)) when calculating the ground heat flux. This is given simply by setting \(\small P_i\) to \(\small P_{AI}\) in equations 2.6, 2.9 and 2.11, such that

\[\begin{equation} \tag{2.13} R_{abs}^{G}=R_{P_{AI}}^{d\downarrow}+R_{P_{AI}}^{b\downarrow}-R_{P_{AI}}^{d\uparrow}+\varepsilon_GR_{P_{AI}}^l \end{equation}\]

where \(\small R_{P_{AI}}^{d\downarrow}\) and \(\small R_{P_{AI}}^{b\downarrow}\) are the downward diffuse and direct radiation flux incident on the ground surface below plant area \(\small P_{AI}\), \(\small R_{P_{AI}}^{d\uparrow}\) is the upward radiation flux at the ground surface, \(\small \varepsilon_G\) is the emissivity of the ground surface and \(\small R_{P_{AI}}^l\) is longwave radiation at the ground surface given equation 2.11 with \(\small \tau_{j,h}\) set to \(\small \tau_{j,0}\) - namely transmission of downward longwave radiation from the top of the canopy to the ground.

Absorption of radiation by leafs

In chapter 1, we introduce the concept of a humid operative temperature, for which it is necessary to compute the radiation absorbed by leafs. The flux density of radiation absorbed by leafs per two-sided leaf area is simply given by

\[\begin{equation} \tag{2.14} R_{abs}=0.5a(R^{d\downarrow}+s_c^LR^{b\downarrow}+R^{d\uparrow})+\varepsilon_LR^l \end{equation}\]

where here \(\small R^{d\downarrow}\), \(\small R^{b\downarrow}\), \(\small R^{d\uparrow}\) and \(\small R^l\) are downward diffuse and direct (beam) radiation, upward diffuse and total longwave radiation incident on the leaf, \(\small a\) is the absorption coefficient for the leaf, \(\small \varepsilon_L\) is leaf emissivity (also its absoption coefficient for longwave radiation) and \(\small s_c^L\) is the fraction of direct beam radiation intercepted by the leafs relative to that which would be intercepted by a flat surface of equivalent area, which for an individual leaf is given by equation 2.3b. For collection of leafs with a known probability density distribution of leaf inclinations this is given by \(\small s_c^L=K\cos Z\). The various radiation fluxes can be determined using the equations presented previously. Leaf emissivities are typically around 0.97.

Absorption of radiation by animals

In chapter 1, we introduce the concept of a humid operative temperature, for which it is necessary to compute the radiation absorbed by an animal. This is not an entirely straightforward task as the absorption of direct radiation depends in part on the animal’s shape, orientation and absorbtance, and the latter may be non-uniform with respect to its orientation. Nevertheless, it is possible to do so by making a number of simplifying assumptions, the primary one being that the animals shape can be approximated by a simpler geometric shape and that absorbtance differs only in a downward and upward direction. With this in mind, a general equation that describes the radiation absorbed by an animal (\(R_{abs}^A\)) above canopy is as follows:

\[\begin{equation} \tag{2.15} R_{abs}^A=a_\downarrow(0.5R_0^{d\downarrow}+s_c^AR_0^{b\downarrow})+0.5a_\uparrow R_0^{d\uparrow}+\varepsilon_A R^l \end{equation}\]

where \(\small a_\downarrow\) and \(\small a_\uparrow\) are absorption coefficients for downward and upward radiation and \(\small \varepsilon_A\) is the emissivity of the animal. In the case of an opaque animal these are given by 1 - its upperside and underside reflectance. Radiation fluxes are as those for a leaf, but when the animal is above vegetation with a known ‘blue-sky’ albedo, simplifies to \(\small R_0^{d\uparrow}=\alpha (R_0^{d\downarrow}+R_0^{b\downarrow})\) and \(\small R^l=R_0^l\)

The complex part is computing \(s_c^A\) - the solar coefficient for the animal. As mentioned briefly previously, there are three ways to conceptualise this coefficient. MOst formally it is the average flux density of direct (beam) radiation absorbed by the animal relative to that which would be absorbed by a surface of equivalent area perpendicular to the solar beam. It can also be conceptualised in terms of the ratio of its silhouette area (\(\small A_p\)) to total area (\(\small A\)), where the silhouette area is also the area of a shadow that would be cast on flat surface perpendicular to the solar beam. While it should in theory be possible to measure this area under a range of solar zenith and azimuth angles, it a further complication is the need to also quantify the total surface area of an animal. It is therefore often more convenient to approximate the animal by some simpler geometric shape, and here we present solutions for four such shapes, namely a sphere, a tri-axial ellipsoid, a elliptic cylinder and and an elliptic cylinder with spheroid-shaped ends. In all cases, \(h\) is the height of the animal, \(\small w\) its width and \(\small l\) its length. For mathematical brevity in subsequent equations we also define the semi-axis of a spheroid and cylinder as \(\small a=0.5l\), \(\small b=0.5w\) and \(\small c=0.5h\). Here only, \(\small a\) is the not an absorption coefficient and \(\small c\) is not the speed of light (cf earlier equations and term definitions and the end of the chapter)!

In the case of the ellipsoid, the vertical semi-axis of the ellipsoid is assumed to be perpendicular to the horizontal - i.e. the ellipsoid is not tilted with respect to horizontal. In the case of the cylindrical shapes, the length of the animal is the height of the cylinder - i.e. it assumed to ‘worm-shaped’, with the ends of the cylinder vertically orientated.

Some of the text below is quite mathematically involved. Readers more interested in application than underpinning theory may simply wish to skip this section and use the function silhoutte which handles all the maths. HOwever, in the case of a sphere (i.e. \(\small h=l=w\)), the solution is quite simple and

\[\begin{equation} \tag{2.15} \frac{A_p}{A}=\frac{\pi h^2}{4\pi h^2}=0.25 \end{equation}\]

In the case of a tri-axial ellipsoid (one in which the height, length and width of the spheroid all differ), from Vickers (1996)

\[\begin{equation} \tag{2.16a} A_p=\pi\sqrt{(l_c^2b^2c^2+m_c^2c^2a^2+n_c^2a^2b^2)} \end{equation}\]

where \(\small l_c\), \(\small m_c\) and \(\small n_c\) are the directional cosines given by \(\small l_c=\cos\theta\cos\phi\), \(\small m_c=\sin\theta\sin\phi\) and \(\small n_c=\cos\phi\), where \(\small \theta\) is the direction of the solar beam relative to longitudinal axis and \(\small \phi\) is the direction of the solar beam relative to the vertical axis. I.e.

\[\theta=\Lambda_a-\Lambda_s\] and \[\phi=\cos^{-1}(\cos Z\cos(-\phi_a)+\sin Z\sin(-\phi_a)\cos\theta)\] where \(\small \Lambda_a\) direction in which the animal is facing and \(\small \phi_a\) its tilt along the longitudinal axis (positive when the fornt is higher than the rear).

Rather surprisingly, there is no general solution for calculating the surface area of a spheroid, but a good approximation, accurate to within 1.061%, is given by Knud Thomsen’s Formula

\[\begin{equation} \tag{2.16b} A\approx 4\pi(\frac{(a^pb^p+a^pc^p+b^pc^p}{3})^{1/p} \end{equation}\]

where \(\small p=1.6075\). Astonishingly, this formula was only ‘discovered’ in 2004, and was simply remarked upon in an email to the author of a mathematical blog as opposed to published!

Vickers (1996) also provides a mathematical solution for the silhouette area of any cylinder-like shape that need not be perfectly elliptic

\[\begin{equation} \tag{2.17a} A_p=A_b\cos\phi+lW(\theta)\sin\phi \end{equation}\]

where \(\small A_b\) is the basal area of the cylinder-like shape and \(\small W(\theta)\) is its width-function defined as the width of the base in the direction \(\theta\). Here the cylinder is orientated such that \(\small \phi\) is the angle of the solar beam relative to the longitudinal axis of the cylinder and \(\small \theta\) is the angle of the solar beam relative to the horizontal axis across the width. I.e.

\[\theta=\tan^{-1}(\frac{cos(Z-\phi_a)}{\sin(Z-\phi_a)\sin(\Lambda_a-\Lambda_s)})\] and \[\phi=\cos^{-1}(\sin(Z-\phi_a)\cos(\Lambda_a-\Lambda_s))\]

In the case of an elliptic cylinder (one in which the cross-sectional height and width are unequal) \(\small A_b=\pi bc\) and \(\small W(\theta)=2\sqrt{b^2\sin\theta+c^2\cos\theta}\)

In the case of a circular cylinder (one in which cross-sectional height and width are equal) this simplifies to \(\small A_b=\pi (w/2)^2\) and \(\small W(\theta)=w\)

As with a tri-axial spheroid, there is no general solution for calculating the surface area of an elliptic cylinder, but a close approximation is given Ramanujan’s formula as

\[\begin{equation} \tag{2.17b} A=2\pi bc+Pl \end{equation}\]

where \(\small P\) is the perimeter of the cylinder approximated by \[P\approx\pi(b+c)(1+\frac{3\lambda_{bc}^2}{10+\sqrt{4-3\lambda_{bc}^2}})\] where \(\small \lambda_{bc}=(b-c)/(b+c)\)

The formulae can also be combined or adapted to determine the solar coefficient of combinations of shapes and are all included in the silhouette function accompanying this book. The particular case where a elliptic cylinder shape has hemi-ellipsoidal instead of flat ends is included with this function. We demonstrated the application of this function in the exercises section below, in which we show how the orientation of an animal with respect to the solar azimuth affects the ration of the silhouette to total surface area.

In the special case of a sphere, neither the solar zenith angle nor the orientation of the animal makes any difference to the silhouette area and in the case of an ellipsoid in which the length and width are equal, only the solar zenith affect the silhouette area. However, for more complex shapes the silhouette area of an animal will change with respect to its orientation and an issue therefore, is that in many situations this is unknown. Broadly, speaking there are three ways in which this issue can be tackled. The first is simply to model the two extreme cases which invariable occur when the animal side-on to or facing the sun. If cross-section area when viewed from above is smaller than the cross-sectional area when viewed side-on, the silhouette area will be greatest when the sun is on the horizon. Conversely, if cross-section area when viewed from above is greater than the cross-sectional area when viewed side-on, the silhouette area will be greatest when the sun is directly overhead.

One might go a step further and assume that the animal is flexible with regards to its orientation and would therefore chose to minimuse or maximise its radiation absorption depending on whether it is too hot or too cold. A third solution, albeit a rather implausible one with respect to the ‘tilt’ of the animal is to assume that it orientates itself entirely at random with respect to the solar zenith and azimuth. In this particular case, from Couchy’s theorem, the expected solhoutte area of any convex body is simply a quarter of the surface area irrespective of its shape and size.

The above formulae are useful for approximating how much radiation is absorbed by convex shaped animals, and can be used in various combinations to do so. However, they are less useful for conical or non-concave shaped animals. For example, with butterflies, internal reflection of radiation from the wings onto the thorax is of some importance. Typically, a finite-element approach is needed (though in some circumstances it may be possible to derive approximate functions that are derived empirically or to create look-up tables). Further treatment of this issue is provided in Rhodes et al (2022).

Examples

1. Deriving leaf reflectance form remotely-sensed imagery.
2. Simulation of the effects of canopy gaps
3. The effects of animal orientation on radiation absorption
4. Why upright animals absorb less radiation than recumbent animals

tme <- as.POSIXlt(c(0:720) * 60, origin = "2022-06-21 06:00:00", tz = "UTC")
jd <- julday(tme$year + 1900, tme$mon + 1, tme$mday)
lt <- tme$hour + tme$min / 60
Z <- solarzen(lt, 49.97, -5.22, jd)
azi <- solarazi(lt, 49.97, -5.22, jd)
# Ellipsoid-shaped animal
adir <- 0  # initial animal direction
sc <- silhouette(Z, azi, adir,  atilt = 0, height = 0.02, width = 0.03, len = 0.07)$solarcoef
par(mar=c(6,6,2,2))
plot(sc ~ as.POSIXct(tme), type = "l", lwd = 2, xlab = "Time", ylab = "Solar coefficient",
     ylim = c(0, 0.5), col = rgb(0,0,0,0.02), cex.axis = 2, cex.lab = 2)
par(new=T)
for (adir in 1:359) {
  sc <- silhouette(Z, azi, adir,  atilt = 0, height = 0.02, width = 0.03, len = 0.07)$solarcoef
  plot(sc ~ as.POSIXct(tme), type = "l", lwd = 2, xlab = "", ylab = "",
     ylim = c(0, 0.5), col = rgb(0,0,0,0.02), cex.axis = 2, cex.lab = 2)
  par(new=T)
}
# Cylinder-shaped animal
par(mar=c(6,6,2,2))
adir <- 0 # initial animal direction
sc <- silhouette(Z, azi, adir,  atilt = 0, height = 0.02, width = 0.03, len = 0.07, shape = "cylinder")$solarcoef
par(mar=c(6,6,2,2))
plot(sc ~ as.POSIXct(tme), type = "l", lwd = 2, xlab = "Time", ylab = "Solar coefficient",
     ylim = c(0, 0.5), col = rgb(0,0,0,0.02), cex.axis = 2, cex.lab = 2)
par(new=T)
for (adir in 1:359) {
  sc <- silhouette(Z, azi, adir,  atilt = 0, height = 0.02, width = 0.03, len = 0.07, shape = "cylinder")$solarcoef
  plot(sc ~ as.POSIXct(tme), type = "l", lwd = 2, xlab = "", ylab = "",
     ylim = c(0, 0.5), col = rgb(0,0,0,0.02), cex.axis = 2, cex.lab = 2)
  par(new=T)
}