Model Equations

Radiation and albedo

The albedo of the vegetated surface is, in effect, the proportion of shortwave radiation that is reflected - i.e. total downward less than absorbed and dependent both on the fraction absorbed by the canopy and by underlying ground surface. Whereas the reflectance of a surface - the ratio of the radiant flux reflected from a material to the incident radiant flux is an innate property of the surface, the albedo of a vegetated surface depends both the ptroperties of the vegetation and the underlying ground surface and also changes in relation to the transmission of radiation through the vegetation and hence also on the angle of the reflecting surfaces relative to the solar beam. It changes therefore in relation to solar angles and also depending on what fraction of the radiation received by the surface is direct radiation relative to diffuse radiation.

From Sellers (1985) and Yuan et al (2017), but with adaptation allowing for radiation to pass through larger gaps in the canopy, albedos are derived using a two-stream radiative transfer model. The two-stream model is also used to derive the radiative flux within the canopy, later used to compute below canopy microclimate. We therefore describe the model in full.

The underpinning assumption of a two-stream model is that solar radiation arrives at the top of the canopy in the form of either diffuse or direct radiation. Because the leaves are partially reflective and translucent, a portion of the direct radiation at any given point within the canopy is scattered in either an upward or downward direction, but because at that point it has interacted with multiple surfaces within the canopy it is essentially diffuse.The a portion of the diffuse radiation is likewise scattered both upward and downward. Consequently there are several radiation streams that need to be derived: (i) Downward direct, (ii) downward diffuse, (iii) upward diffuse, (iv) the contribution of downward direct to downward diffuse and (v), the contribution of downward direct to upward diffuse. There is is no upward direct, as any radiation upward is assumed to be scattered isotropically and therefore diffuse. The albedo is essentially determined from the upward and downward streams at the top of the canopy, but it is also possible to derive upward and downward streams within the canopy.

Downward direct radiation

It is convenient to envisage the canopy being divided into multiple layers, each with a series of nodes \(\small i\) situation at the top of each layer and numbered sequentially downward from the top of the canopy. Let us combine the effects of leafs and other canopy elements such as branches and tree trunks 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. From e.g. Campbell (1985), the attenuation of radiation at node \(\small i\) is then described as \(\small R_i^{\downarrow b}=R_0^{\downarrow b}\exp(-KP_i)\) where \(\small R_i^{\downarrow b}\) is the downward flux of direct radiation that would be intercepted by a surface perpendicular to the solar beam at node \(\small i\), \(\small \downarrow R_0^{\downarrow b}\) the is the downward flux of direct radiation intercepted perpendicular to the solar beam at the top of the canopy, and \(\small K\) is an radiation extinction coefficient, which depends on the inclination angle of canopy elements relative to the solar beam (vertically orientated leaves will block out more sunlight when the sun is low above the horizon).

However, 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

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

where \(\small Z\) is the solar zenith angle (see below). 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.

This assumes canopy elements are small and randomly distributed. In reality, canopy elements are non-randomly distributed and radiation can penetrate through larger gaps in he canopy unimpeded. We thus 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. 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. It is assumed 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.

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

\[R_i^{\downarrow b}=R_0^{\downarrow b}[{\tau_i}^b+(1-{\tau_i}^b)\exp(-\tilde{K}\tilde{P_i})]\] where \(\small {\tau_i}^b=F_G^{m_zK_C}\) is the fraction of direct radiation transmitted through larger canopy gaps, \(\small K_C=\tilde{K}/K(0)\) and \(\small \tilde{K}=K(Z)\cos Z/s_c\), where \(\small K(Z)\) and \(\small K(0)\) are the values of \(\small K\) derived using the equation above, when the solar zenith angles are \(\small Z\) and \(\small 0\) respectively, \(\small m_z=(h-z_i)/h\) adjust the gap transmission for height above ground \(\small z_i\) relative to total canopy height \(\small h\) and \(\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)\) and \(\small s_c\) is the fraction of direct beam radiation intercepted by an inclined surface given by

\[s_c=\cos Z \cos s+\sin Z\sin s\cos(\Lambda_s-\Lambda)\] where \(\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). Slope and aspect are given as user inputs.

The position of the sun

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

\[\cos Z=\sin\delta\sin \phi_y+\cos\delta\cos\phi_y\cos h_s\] where \(\small h_s\) is the hour angle in solar time (\(\small t_s\)) given by

\[h_s=(\pi⁄12)(t_s-12)\] and \(\small \delta\) is the current declination angle of the sun calculated from the Astronomical Julian day (\(\small j\)) as

\[\delta=\frac{\pi 23.5}{180} \cos\lgroup 2\pi\frac{j-159.5}{365.25}\rgroup\] 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

\[t_s=t_l+\frac{4(\phi_x-\phi_m)-\Delta_t}{60}\] 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)

\[\Delta_t=-7.659\sin(M+9.683)\sin(2M-3.5932)\] where \(\small M\) is the mean anomaly given by \(\small 6.24004077+0.01720197(j-2451545)\).

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

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

Diffuse radiation

The diffuse radiation streams can be written using two sets of differential equations:

\[-\frac{dR_i^{d\uparrow}}{dP}=-(a+\gamma)R_i^{d\downarrow}+sR_i^{b\downarrow}\] \[\frac{dR_i^{d\downarrow}}{dP}=-(a+\gamma)R_i^{d\downarrow}+\gamma R_i^{d\uparrow}+s^,R_i^{b\downarrow}\] 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.

Each of the equation’s components describes a different physical process. In both equations, the left hand term describes the overall attenuation of upward or downward radiation through the canopy. The first right hand term in the top equation describes the fraction of upward diffuse radiation that is re-scattered in upward direction and the equivalent in the bottom equation describes the fraction of downward diffuse radiation that is re-scattered downward. \(\small a\) is the absorption coefficient for incoming diffuse radiation per unit leaf area and \(\small \gamma\) is the backward scattering coefficient. The second term on the right hand side in the top equation describes the fraction of the downward diffuse radiation that is converted into an upward diffuse flux by back-scattering and the equivalent in the bottom equation describes the fraction of the upward diffuse radiation that is converted into an downward diffuse flux by back-scattering.The final term on the right hand side of the top equation 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 the bottom equation 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.

\(\small a\) is given by \(\small a=(1-\omega)\) where \(\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, both provided as user inputs.

\(\small \gamma\) is given by \(\small \gamma=0.5(\omega+J\delta)\) where \(\small \delta\) is given by \(\small \delta=\alpha_P-\tau_P\) and \(\small J\) is an integral function of the inclination distribution of canopy elements approximated by \(\small J=\cos^2 \hat{\Theta}_P\) where \(\small \hat{\Theta}_P\) is the mean inclination angle of canopy elements in the zenith direction. From Campbell (1990) \(\small \hat{\Theta}_P=9.65(3+x)^{1.65}\) and from Verhoef (1984) and Pinty et al (2006) \(\small s=0.5(\omega+\frac{J\delta}{K})K\). The forward scattering coefficient for incident direct radiation, \(\small s^,\) is given by \(\small s^,=\omega K-s\).

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. A full explanation of their derivation is given in Sellers (1985) and Yuan et al (2017). Only the final equations are presented here. These are as follows:

\[R_i^{d\uparrow}=R_0^{d\downarrow}[(1-\tau_i^d)(p_1\exp(-h\tilde{P_i})+p_2\exp(h\tilde{P_i}))+\alpha_G\tau_i^d]+R_b^{d\uparrow}\] \[R_i^{d\downarrow}=R_0^{d\downarrow}[(1-\tau_i^d)(p_3\exp(-h\tilde{P_i})+p_4\exp(h\tilde{P_i}))+\tau_i^d]+R_b^{d\downarrow}\] where \(\small R_i^{d\uparrow}\) is upward diffuse radiation and \(\small R_i^{d\downarrow}\) is downward diffuse at \(\small i\), \(R_0^{d\downarrow}\) is downward diffuse at the top of the canopy, \(\small \alpha_G\) is ground reflectance. Here \(\small \tau_d={F_G}^{2m_z}\) is transmission of diffuse radiation through larger gaps in the canopy and \(\small R_i^{db\uparrow}\) and \(\small R_i^{db\downarrow}\) are the contribution of direct radiation to the diffuse upward and downward flux respectively given by

\[R_i^{db\uparrow}=R_0^{b\downarrow}[(1-\tau_i^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_i^b]\] and

\[R_i^{db\downarrow}=R_0^{b\downarrow}[(1-\tau_i^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})]\] The equations for calculating albedos are given by

\[\alpha_d=(1-\tau_0^d)(p_1+p_2)+\tau_0^b\] and

\[\alpha_b=(1-\tau_0^d)(\frac{p_5}{\sigma_p}+p_6+p_7)+\tau_0^d\] where \(\small \alpha_d\) is white-sky albedo (i.e. Bi-Hemispherical Reflectance: that reflecting diffuse radiation), \(\small \alpha_b\) is black-sky albedo (i.e. Directional Hemispherical Reflectance - that reflecting direct radiation) and \(\small \tau_0^d\) and \(\small \tau_0^b\) are the values of \(\small \tau_i^d\) and \(\small \tau_i^b\) with \(\small m_z=1\). The total albedo is given by

\[\alpha=f_d\alpha_d+(1-f_d)\alpha_b\] where \(\small f_d\) is the diffuse radiation fraction given by \(\small R_0^{d\downarrow}/R_{sw}^\downarrow\).

The derivation of the terms \(p_{1..10}\) is given below.

Here \(\small \alpha^*_G\) is the bluesky albedo of the gorund surface, effectively the the ground reflectance \(\small \alpha_G\) adjusted for inclination given by \(\small \alpha_G\cos Z/s_c\)

Longwave and emitted radiation

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, the longwave radiation streams within the canopy can be thought of as originating from three sources. Firstly from the ground surface, where incident longwave 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.

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 downward flux of radiation \(\small R_j^{\downarrow lw}\) at \(\small j\) is given by

\[R_j^{\downarrow lw}=\tau_{j,h}^dR_{lw}^\downarrow+\sum_{i=j}^{i=n}\Delta \tilde{P}_i \tau_{i,j}\varepsilon_i\sigma T_i^4)\] and the upward flux \(\small R_j^{\uparrow lw}\) is given by

\[R_j^{\uparrow lw}=\tau_{0,j}^d\varepsilon_G\sigma T_G^4+\sum_{i=0}^{i=j}\Delta \tilde{P}_i \tau_{i,j}^d\varepsilon_i\sigma T_i^4)\] where \(\small \tau_{i,j}^d\) is the transmission of longwave radiation between \(\small i\) and \(\small j\) given by

\[\tau_{i,j}^d=F_G^{2m_{i,j}}+(1-F_G^{2m_{i,j}})\exp(-\tilde{P}_{i,j})\] where \(\small m_{i,j}=|z_i-z_j|/h\) adjust the grap transmission for distance from \(\small i\) to \(\small j\) and \(\small \tilde{P}_{i,j}=|\tilde{P}_i-\tilde{P}_j|\) is the is the total one sided plant area between \(\small i\) and \(\small j\). Here \(\small \varepsilon_G\) and \(\small \varepsilon_i\) are the emissivities of the ground surface and canopy elements respectively, \(\small \sigma\) the Stefan–Boltzmann constant and \(\small T_G\) and \(\small T_i\) are absolute temperatures of the ground surface and canopy elements respectively in Kelvin. Here plant areas \(\small \tilde{P}\) refer to the plant areas concentrated into non-gaps.

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. Resultantly the model is run iteratively.

Emitted radiation is given by

\[R_{lw}^\uparrow=\varepsilon_C\sigma T_C^4\]

The linearising the emitted radiation, used to solve for temperature in the Penman-Monteith equation assumes

\[\varepsilon_C\sigma T_C^4\approx\varepsilon_C\sigma (T_A^4+4T_A^3(T_C-T_A) =\varepsilon_C\sigma T_A^4+\frac{\tilde{\rho}c_p}{r_R}(T_L-T_A)\] where the radiative resistance \(\small r_R\) is given by

\[r_R=\frac{\tilde{\rho}c_p}{4\varepsilon_C\sigma T_A^3}\] ### Radiation absorption by the ground and by leaves The two-stream model also used to calculate radiation absorbed by the ground surface (\(\small R_{abs}(G)\)), given by:

\[R_{abs}(G)=(1-\alpha_G)(R_{i=0}^{d\downarrow}+s_cR_{i=0}^{b\downarrow})+\varepsilon_GR_{j=0}^{\downarrow lw}\]

where \(\small \varepsilon_G\) is the emissivity of the ground surface. Radiation absorption by canopy elements at height \(\small i\), averaged over the upward and downward-facing surfaces, is given by

\[R_{abs}(L_i)=(1-\alpha_L)(0.5(R_i^{d\downarrow}+R_i^{d\uparrow})+0.5s_L(R_i^{\downarrow b}+R_i^{\uparrow b})+0.5\varepsilon_L(R_{j=i}^{\downarrow lw}+R_{j=i}^{\uparrow lw})\] where \(\small \varepsilon_L\) is the emissivity of canopy elements and \(\small s_L = K/\cos Z\)

Diabatic influencing factors and coefficients

From Harman & Finnigen (2008) the diabatic influencing factors and coefficients feature when the canopy surface is strongly heated or cooled and are given by

\[\psi_M(z)=\psi_M(\frac{z_M}{L_O})-\psi_M(\frac{z-d}{L_O})\] \[\psi_H(z)=\psi_H(\frac{z_H}{L_O})-\psi_H(\frac{z-d}{L_O})\]

where from Businger et al (1971)

\[\psi_M(\zeta)=\ln[(\frac{(1+x)}{2})^2(\frac{1+x^2}{2})]-2\arctan x+\frac{\pi}{2}\] with \(\small x=(1-15\zeta)^{1/4}\) under unstable conditions (\(\small H>0\)) and

\[\psi_M(\zeta)=-4.7\zeta\] under stable conditions (\(\small H<0\)). The equivalent formulae for heat are given by

\[\psi_H(\zeta)=\ln[(\frac{1+y}{2})^{2}]\] with \(\small y=\) With \(\small y=1⁄\phi_H =\sqrt{1-9\zeta}\) under unstable conditions and \(\small \psi_H(\zeta)=-4.7\zeta/0.74\) under stable conditions. Here \(\small L_0\) is the Obukhov length given by

\[L_O=-\frac{\tilde{\rho}c_pu_*^3\overline{T}}{\kappa gH}\] where \(\small g\) is the gravitational constant (~9.81) and \(\small \overline{T}\) is the average temperature of the height profile (in Kelvin), taken as the average of the canopy and the air.

From Yasuda (1988)

\[\sqrt{\phi_H}=\frac{1}{(1-16\zeta)^{1/4}}\] under unstable conditions and

\[\phi_H=1+\frac{6\zeta}{1+\zeta}\] under stable conditions.

To avoid unrealistic reversals of the temperature profiles caused by an extreme diabatic influence, the diabatic coeffcients are capped such that

\[|\psi_M|<0.9\ln\frac{z_R-d}{z_M},\space\space\space |\psi_H|<0.9\ln\frac{z_R-d}{z_H}\] Additionally, to avoid the diabatic influencing factor for heat resulting in a negative roughness length for momentum, it is capped to

\[\psi_H<0.9\frac{\kappa}{\beta}\]

Soil thermal properties

From Campbell (1986), thermal conductivity (\(\small k_s\)) is given by

\[k_s=c_1+c_2S_M-(c_1-c_4)\exp(-c_3S_M^4)\] where

\[c_1=\frac{0.57+1.73Q_v+0.93M_v}{1-0.74Q_v-0.49M_v}-2.8(Q_v+M_v)(1-(Q_v-M_v))\] \[c_2=1.06\rho_DS_M,\space\space c_3=1+2.6C_M^{-1/2}, \space\space c_4=0.03+0.7(M_v+Q_V)^2\] where \(\small \rho_D\) is the bulk density of soil (\(\small Mg·m^{-3}\)), \(\small Q_v\) and \(\small M_v\) the volumetric quartz and mineral content of the soil respectively, \(\small C_M\) the mass fraction of clay and \(\small S_M\) the volumetric water fraction of the soil.

The volumetric density (\(\small \rho_s\)) and specific heat (\(\small c_s\)) of the soil are then given as follows

\[\rho_s=\rho_D(1-S_M)+1000S_M,\space\space c_s=\frac{2400\rho_D}{2.64}+4180S_M\] Typical values of \(\small \rho_D\), \(\small Q_v\), \(\small M_V\) and \(\small C_M\) for given soil types are shown in Table 1

Temperatures below ground

In the soil, heat storage is significant, and both the annual and diurnal temperature cycles are attenuated and occur later in the day or year. Assuming that the annual and diurnal cycles are approximately sinusoidal, that the soil layer is infinitely deep and with uniform thermal properties, from Campbell & Norman (2012) and de Vries et al (1967), the temperature at depth \(\small z\) and time \(\small t\) is given by

\[T_z(t)=\overline{T_G}+A_D(0)\exp(\frac{z}{D_D})\sin((\omega_D(t-t_0(D))+\frac{z}{D_D})+A_A(0)\exp(\frac{z}{D_A})\sin((\omega_A(t-t_0(A))+\frac{z}{D_A})\] where \(\small \overline{T_G}\) is mean annual temperature of the ground surface (K), \(\small z\) is depth below ground (m) and negative, \(\small A_D\) and \(\small A_A\) are the amplitude in the diurnal and annual temperature cycles and \(\small D_D\) and \(\small D_A\) are the diurnal and annual damping depths. It can be shown empirically, that the above relationship is approximated well by the following relationship:

\[T_z(t)=\delta(\frac{1}{n}\sum_{i=t-n+1}^{t}T_G(i))+\overline{T_G}(1-\delta)\] where \(\small \frac{1}{n}\sum_{i=t-n+1}^{t}T_G(i)\) is the rolling mean of the ground surface temperature over the last \(\small n\) hours. The number of hours over which the rolling mean is calculated is contingent on depth and is given by

\[n=\frac{24z}{\pi D_D}\] and the weighting parameter \(\small \delta\) is given by

\[\delta=0.00069\frac{z}{D_D}+0.87142\] Comparisons of this approach with those produced by the full multi-layer soil model described in Campbell (1986) and included in the NicheMapR package are shown in Fig. 1.

Soil water content

This package implements a modified version of the Mahrt and Pan (1984) two-layer model of soil hydrology. The volumetric soil moisture fraction at each daily time-increment \(\small j\) in a top shallow soil layer is computed as

\[\theta_j=\theta_{j-1}-s_1 P-s_2 R_{netp}\] where \(\small \theta_{j-1}\) is soil moisture in the previous time increment, \(\small P\) is daily precipitation, \(\small R_{netp}\) is the daily mean of hourly positive net radiation. Soil moistures are capped so that \(\small θ_j\) cannot exceed the volumetric water fraction at saturation (\(\small S_{mx}\)) or drop below the residual volumetric water fraction for a given soil type (\(\small S_{mn}\)).The assumption is thus that daily evaportransporation is linearly related to \(\small R_{netp}\) and thus that the Bowen ratio remains relative constant during the day when net radiation budget is positive.

At each daily time increment, water is then exchanged with an underlying deeper soil layer of 10 times the volume as follows:

\[\theta_j=\theta_{j-1}+s_3 k(\theta_j^d-\theta_j )\] \[\theta_j^d=\theta_{j-1}^d+0.1s_3 k(\theta_j-\theta_j^d )\] where \(\small \theta_j^d\) is the soil moisture in underlying deeper soil layer at time \(\small j\) and \(\small k\) is the hydraulic conductivity given by

\[k=k_{sat}(\frac{\overline{\theta_j}}{\theta_S})^{s_4}\] where \(\small k_{sat}\) is saturated hydraulic conductivity, \(\small \overline{\theta_j}\) is the mean volumetric water content of the two soil layers given by \(\small (\theta_j+\theta_j^d)/2\). The terms \(\small s_{1..4}\) are coefficients estimated by iteration for given soil types using the outputs the more complex multi-layer hydrological model given in Campbell (1986) and included in the NicheMapR package (Kearney & Porter 2017). Examples of the model output in comparison to a more complex multi-layer model in NicheMapR are shown in Fig. 1.

Typical values of \(\small S_{mx}\), \(\small S_{mn}\) and \(\small k_{sat}\) and estimated values for \(\small s_{1..4}\) for given soil types are shown in Table 1

Conceptual overview

In contrast to above canopy, it is now widely recognised that the transfer of heat and vapour below canopy cannot be predicted using a MOST approach (e.g. Bonan, et al. 2021). To replace MOST in this context, Raupach (1989a,b) developed an analytic Lagrangian theory, which predicts the concentration heat and vapour emanating from a spatially extensive source downwind of the point of interest. In this ‘localized near-field’ theory (LNZT), the mean concentrations of heat or vapour are expressed as the sum of a diffusive far-field contribution that obeys MOST, and a non-diffusive near-field contribution, which is determined from local sources by assuming the turbulence to be locally homogeneous.

To recap formulae presented above:

\[\tilde{\rho}c_pT_z=C_f(z)+C_n(z)\] \[C_f(z)=\tilde{\rho}c_pT_h-C_n(h)+\sum_{z}^{h}\frac{H(z_i)}{K_H(z_t)}dz\] \[H(z_i)=\tilde{\rho}c_p/r_L(T_{c_i}-T(z_i))+F_G\] where \(\small C_f(z)\) and \(\small C_n(z)\) are the far-field and near-field contributions, \(T_h\) is air temperature at the top of the canopy (\(\small C_f(h)=\tilde{\rho}c_pT_h\) scales the source concentration to temperature), \(\small K_H\) is thermal diffusivity, \(\small r_K\) resistance to heat loss by canopy elements such as leaves with temperature \(\small T_{c_i}\) and \(\small F_G\) is the ground heat flux.

The formulae can seem quite confusing, but essentially capture four processes. First, near the top of the canopy temperatures will closely match temperatures at the top of the canopy because \(\small \sum_{z}^{h}\frac{H(z_i)}{K_H(z_t)}dz\to0\) as \(\small z_i\to h\) (i.e. there is nothing to sum over). Second, as one descends through the canopy, temperatures will typically increase if the sensible heat flux from canopy elements is positive as would typically be the case when the net radiation budget is strongly positive. This is because \(\small \sum_{z}^{h}\frac{H(z_i)}{K_H(z_t)}dz\) increases as \(\small z\to0\). This seems a bit counter-intuitive as one typically imagines increased shading to result in cooler temperatures, but the key point is that heat arriving at height \(\small z_i\) is emanating from canopy elements at all heights. However, near the top of the canopy, temperatures are more closely coupled to the air above it. The cooling effect of canopies is the result of two processes: lower down in the canopy, negative sensible heat flux from individual canopy elements that are shaded is possible even if the total flux from the canopy is positive. More importantly though, the ground heat flux comes into play: if the ground surface is cooler than the leaves at the top of the canopy there is a strong cooling effect simply because \(\small F_g\to\infty\) as \(\small z\to0\) (see below). The last effect arises due the near-field contribution given by

\[C_n(z)=\int_{0}^{\infty}\frac{S(z_i)}{\sigma_\omega(z)}k_n[\frac{z-z_i}{\sigma_\omega(z_i)T_L}+\frac{z+z_i}{\sigma_\omega(z_i)T_L}]dz_i\] where \(\small S(z_i)\) is the source concentration given by \(\small f_d(z_i)H(z_i)\) where \(\small f_d(z_i)\) is foliage density for layer \(\small z_i\) and \(\small k_n\) is a kernal function. The important take homes from this formula is canopy elements heated by the sun will warm the air immediately surrounding the elements with the effect influenced strongly by the density of canopy foliage, and diminishing very sharply with distance (as dictated by the shape of the kernal function).

The Langrangian timescale problem

The term Langrangian essential means that the heat or vapour in a canopy can be considered as being emitted from (or absorbed by) a large number of approximately point sources (or sinks),the individual leaves. The spread of the plume from each leaf, is determined by the motion of heat or vapour ‘particles’ emanating from the leaf in question. Of key importance is the effect of the persistence of the motions within the canopy on the dispersion properties of these marked particles, ultimately determined by two parameters: the particle position variance (\(\small \sigma_\omega^2\)) and the Lagrangian time scale \(\small T_L\), both of which appear in the formulae above, but neither of which are actually known.

However, plausable vertical profiles are propsoed as:

\[T_L=\frac{a_2h}{u_*}\] \[\sigma_\omega(z_i)=u_*[0.5(a_1+a_0)+0.5(a_1-a_0)\cos(\pi(1-\frac{z_i}{h}))]\] with \(\small a_1=1.25\) (the value of \(\small \sigma_\omega(z_i)/u_*\) at the top of the canopy at height \(\small h\)) and \(\small a_0=0.25\). The parameter \(\small a_2\), though generally taken to be vertically constant within the canopy is effectively unknown. However, although it is questionable whether conventional MOST and LNZT theory can be reconciled in this way, following several authtos (e.g. Ogee et al 2003) a solution to the problem is found by making thermal diffusivity equivalent to its above canopy formulation at height \(\small z=h\). Since above canopy

\[K_H=\kappa u_*(z-d)/\phi_H\] and below canopy

\[K_H=\sigma_\omega(z_i)^2T_L\] Setting \(\small a_1=1.25\) and \(\small a_0=0.25\) It follows that

\[\sigma_\omega(z_i)=\frac{a_2h}{u_*}[u_*(0.75+0.5\cos(0))]^2\] Thus

\[a_2=\frac{\kappa(1-\frac{d}{h})}{1.25^2}\phi_H\]

Ground heat flux

In micropoint, ground surface temperature \(\small T_G\) is inferred prior to running the below canopy model, using the Big-leaf model. Thus, from Fick’s law:

\[F_G=\frac{\tilde{\rho}c_p}{\overline{K_H^{-1}(z_i)}}(T_G-T_{z})\] where \(\small T_z\) is air temperature at height \(\small z\) and \(\small \overline{K_H^{-1}(z_i)}\) is the mean inverse thermal diffusivity between heights \(\small z_i=0\) and \(\small z_i=z\).

Snowpack energy balance

The energy balance of both the ground and canopy snow models are given by

\[R_{abs}-R_{lw}^\uparrow-H_S-L_S-G_s-F_M=0\] where \(\small R_{abs}\) is radiation absorbed by the snowpack, \(\small R_{lw}^\uparrow\) is radiation emitted by the snowpack, \(\small H_S\) is the sensible heat flux, \(\small L_S\) is the latent heat flux due to sublimation, \(\small G_S\) is the rate of heat storage by the snowpack and \(\small F_M\) is the energy flux removed from the pack during melt.

Radiation budget

Radiation absorbed by the canopy is given by

\[R_{abs}=(1-\alpha_s)R_{sw}^\downarrow+0.97 R_{lw}^\downarrow\] where \(\small a_s\) is snow surface albedo calculated as

\[\alpha_s=-9.874\ln(d_s+78.3434)\]

where \(\small d_s\) is the number of decimal days since he last snowfall. other terms are as defined for the canopy big-leaf model.

To calculate the temperature of the snow-covered ground surface, radiation is adjusted to accommodate the presence of the canopy. The shortwave component of absorption is calculated using the two-stream model described above, but with the following adjustments made. First, snow is assumed to absorb radiation isotropically and in consequence the absorption of direct radiation is assumed independent of ground slope and aspect. Thus, the coefficient \(\small x\) is set to one irrespective of leaf inclination angles. Second, both ground (\(\small \alpha_G\)) and leaf (\(\small \alpha_P\) reflectance are set to snow \(\small \alpha_S\) snow albedo. Third, canopy element transmittance (\(\small \tau_P\)) is adjusted to accommodate the additional presence of snow, and is adjusted such that transmittance both through vegetation and snow is given by

\[\tau=\tau_P\tau_S,\space\space \tau_s=\exp(-k_{S}z_{i_s})\] where from Warren et al (2006) \(\small k_s=0.1\) in an optical extinction coefficient for radiation transmitted through snow and \(\small z_{i_s}\) is the average snow water equivalent thickness of snow within each canopy layer given by

\[z_{i_s}=\frac{h_{sc}\rho_s}{(h-h_{sg})\rho_w}\] where \(\small h_sG\) is the depth of ground snow (m) and \[\small h_sC\] is the dpeth of canopy snow (m).

Lastly, the plant area index is adjusted to accommodate both the presence of both snow lying on the ground and distributed through the canopy as follows

\[P_{AI}^S=\frac{h-g_Fh_s}{h}P_{AI}+\frac{(1-g_F)h_s}{h-g_Fh_s}\] The absorption of longwave radiation is given by

\[R_0^{\downarrow lw}=\tau_c0.97R_{lw}^\downarrow+(1-\tau_c)0.97^2\sigma{T_s}^4\] where \(\small \tau_c\) is calculated as for the plant canopy, but with \(\small P_{AI}\) and \(\small h\) adjusted for now.

transmission of longwave radiation through the canopy given by

\[R_0^{\downarrow lw}=\tau_0^d0.97R_{lw}^\downarrow+(1-\tau_0^d){0.97}^2\sigma T_s^4\] where \(\small \tau_0^d\) is calculated as for the canopuy without snow, but wiht \(\small P_{AI}\) and \(\small h\) adjusted for snow.

Where \(\small h_s>h\) \(\small R_{}\) is calculated as for the bigleaf snow model - i.e. canopy effects are ignored

Radiation emitted by the snowpacks are given by

\[R_{lw}^\uparrow=0.97\sigma {T_{S_C}}^4\] and

\[R_{lw}^\uparrow=0.97\sigma {T_{S_G}}^4\]

where \(\small T_{S_C}\) and \(\small T_{S_G}\) are snow surface temperatures for the ground and canopy derived by solving the energy balance equation using the Penman-Monteith equation.

Sensible heat fluxes

The sensible heat fluxes are given by

\[ H_S=\frac{\tilde{\rho}c_p}{r_{Ha}}(T_{S}-T_A)\] where, as for the canopy \(\small T_s=T_{S_c}\) and for the ground \(\small T_S=T_{S_G}\) and \(\small r_{Ha}\) is the bulk surface aerodynamic resistance of the snow pack given by

\[r_{Ha}=\frac{\ln\frac{z_R-d}{z_H}+\psi_H(z_R)}{\kappa u_*}\] where \(\small u_*\) is the friction velocity of wind given by \(\small u_*=κu_{z_R}/[\ln\frac{z_R-d}{z_M}+\psi_M(z_R)]\). The key distinctions are that the zero plane displacement height and roughness lengths accommodate the presence of snow and are thus given by

\[d=h\frac{1-(1-\exp(-\sqrt{7.5P_{AI}^S}))}{\sqrt{7.5P_{AI}^S}},\space\space h_s<h\] \[d=h_s,\space\space,h_s\ge h\] where \(\small h\) is vegetation height (m), \(\small h_s\) is snow depth given by \(\small h_s=S\rho_w/\rho_s\) where \(\small S\) is snow water equivalent (m), \(\small \rho_w=1000 \space kg·m^{-3}\) is the density of water and \(\small \rho_s\) is the density of snow. Here and \(\small P_{AI}^S\) is the snow-covered plant area index above the snow calculated as described below. Following Sturm et al. (2010) and Kearney (2020), snow density is computed as

\[\rho_s=(\rho_{max}-\rho_0)(1-\exp(-k_1h_s-k_2a_s)+\rho_0\]

where \(\small \rho_{max}\) is the maximum allowable density, \(\small \rho_0\) is the initial density, \(\small a_s\) is snow age (in decimal days) and \(\small k_1\) and \(\small k_2\) are fitting parameters, derived form the snow environment (Table 2).

The roughness length for momentum is computed as

\[z_M=(h-d)\exp[\frac{-\kappa}{\beta}-\psi_H(z_R)],\space\space h>h_s\] \[z_M=0.0005,\space\space h_s\ge h\] where \(\small Be=\sqrt{0.003+0.1P_{AI}^S}\)

Latent heat fluxes

Latent heat fluxes are given by

\[L_S=\frac{\lambda\tilde{\rho}}{r_{Ha}p_A}(e_S-e_A)\] where \(\small e_S\) is the effective vapour pressure of either the canopy or ground surfaces calculated from snow temperatures using Tetens equation. \(\small lambda\) is the latent of sublimation of the snow pack temperature is below zero.

Rate of heat storage by snowpack

As for the diurnal ground heat flux cycle, the annual rate of heat storage by the canopy and ground snow packs are given by

\[G_s=\frac{\sqrt{2}A_{T_S}(0)k_s\sin(\omega(t-t_0+\pi/4)}{D_S}\] where \(\small A_{T_S}\) is the amplitude of the diurnal snow surface temperature cycles, \(\small k_s\) is snow thermal conductivity (\(\small W·m^{-1}·K^{-1}\)), \(\small \omega\) is the angular frequency given by \(\small \omega_D=2\pi/(24 \times 3600)\) and \(\small D_S\) is the damping depth given by \(\small D_S=\sqrt{2\kappa/\omega}\) where \(\small \kappa\) is thermal diffusivity given by \(\small \kappa=k/\rho_sc_s\) where \(\small \rho_s\) and \(\small c_s\) are the volumetric density and specific heat of the snow pack. Snow thermal conductivity is estimated from Djachkova’s formula, following Anderson (2006) and Kearney (2020)

\[k_s=0.0442\exp(5.181\rho_s)\] The Energy Balance equation, excluding \(\small F_M\) is then solved for temperature using the Penman-Monteith equation. If snow surface temperature exceeds zero, the snow pack is then assumed to melt until temperatures attain zero, so snow pack temperature is set to zero. The energy flux removed from the pack during melt is thus given by:

\[F_M=\frac{(T_s-273.15)c_s\rho_wS}{\Delta_t}\]

where \(\small c_s\approx 2100 J·kg^{-1}·K^{-1}\) is the specific heat of ice and \(\small Delta_t=3600\space s\) is the time-step of the model. Since \(\small \rho_w\) is also approximately constant, the equation simplifies to

\[F_M=583.3(T_s-273.15)S\] The relevance of explicit calculation of \(\small F_M\) is the calculation of the mass balance.

Snowpack mass balance

The mass balance of the canopy and ground snow models are given by

\[\Delta S_C=I_C-(1-f_G)M_R-((1-f_G)\frac{L_s}{\rho_w\lambda_s}+\frac{F_M}{\rho_w\lambda_f})\Delta_t\] \[\Delta S_G=(P_s-I_C)-f_GM_R-(f_G\frac{L_s}{\rho_w\lambda_s}+\frac{F_M}{\rho_w\lambda_f})\Delta_t\] where \(\small \Delta S_C\) and \(\small \Delta S_G\) are the change in canopy and ground snow balance (m snow-water equivelent), \(\small P_S\) is precipitation (m) falling as snow (assumed to do so when air temperature is below 0 degrees C), \(\small I_C\) is canopy interception of snow, \(\small \lambda_s\) and \(\small \lambda_f\) are the latent heat of sublimation and fusion respectively here measured in \(\small J·kg^{-1}\) where the conversion is given by \(\small \lambda(J·kg^{-1})=\lambda(J·mol^{-1})/M_w\) where \(\small M_w=18.015×10^{-3} kg·mol\) is the molar mass of water. Here \(\small M_R\) is rain melt calculated following Anderson (2006) and Kearney (2020) as

\[R_M=0.0125(T_A-273.15)P,\space\space T_s>273.15\] \[R_M=0,\space\space T_A\le273.15\] \(\small f_G\approx\exp(-P_{AI}^S)\) partitions sublimation and rain melt between the canopy and the underlying round layer (the portioning of melt is handled explicitely as it is contingent on mas balance of the snowpacks)

Canopy interception

Following Hedstrom & Pomeroy (1998) snow interception by the canopy is given by

\[I_c=0.678(L^*-L_{t-1}(1-\exp(-\frac{C_f(0)}{L^*} P))\] where \(\small L^*\) is the maximum canopy snow load (mm snow water equivelent), \(\small L_{t-1}\) the snow load in the previous time-step and \(\small C_f(0)\) is effective canopy cover perpendicular to the direction of snow fall given by \(\small C_p(0)=1-\exp(-KP_{AI}^A)\) where \(\small K=1/2\cos(Z_s)\) where \(\small Z_s\) is the zenith angle of the direction of snow fall given by \(\small \arctan(\overline{u_c}/0.8)\) where 0.8 m/s is the terminal velocity of snow fall and \(\small \overline{u_c}\) is mean canopy wind speed derived by integrating the wind speed profile below canopy such that

\[\overline{u_c}=\frac{8u_h\beta^2}{P_{AI}^A}(1-\exp(-\frac{P_{AI}^A}{8\beta^2}))\] where \(\small u_h\) the wind speed at the top of the canopy is given by \(\small u_h=u_*/\beta\).

Microclimate model with snow

The snow microclimate model forms an extension of the non-snow model and its formulation is broadly similar, but on hours with snow present, the following adjustments are made.

1. As for radiation absorption canopy elements, ground and leaf reflectance are replaced by snow albedo, the coefficient \(\small x\) is set to 1 and foliage density, canopy height and leaf transmittance adjusted to accommodate the presence of snow.

2. When calculating the latent heat flux from canopy elements in the small leaf model, the surface is assumed to be free evaporating or sublimating and such that resistance to heat and vapour loss are identical.

3. When calculating the ground heat flux in the small leaf model, ground snow surface temperature rather than ground surface temperature is used and \(\small z_i\) is adjusted to \(\small z_i-g_Fh_s\) where \(\small g_Fh_s\) is the depth of the ground surface layer

4. In applying the small leaf model, heat fluxes are summed only for those canopy layers that are above the height of ground snow.

5. If the height for which model outputs are requiired lies below the snow pack but above the ground surface, temperatures are calculated as if there were below ground but wiht the thermal properties of the gorund replaced by the thermal properties for snow.Thus

\[T_z(t)=\delta(\frac{1}{n}\sum_{i=t-n+1}^{t}T_S(i))+\overline{T_S}(1-\delta)\]

where \(\small n\) is given by

\[n=\frac{24(g_fh_s-z)}{\pi D_S}\] and \(\small \delta\) is given by

\[\delta=0.00069\frac{z}{D_S}+0.87142\]