2.3. Surface Albedos

2.3.1. Canopy Radiative Transfer

Radiative transfer within vegetated canopies is calculated from the two-stream approximation of Dickinson (1983) and Sellers (1985) as described by Bonan (1996)

(2.3.1)\[-\bar{\mu }\frac{dI\, \uparrow }{d\left(L+S\right)} +\left[1-\left(1-\beta \right)\omega \right]I\, \uparrow -\omega \beta I\, \downarrow =\omega \bar{\mu }K\beta _{0} e^{-K\left(L+S\right)}\]

(2.3.2)\[\bar{\mu }\frac{dI\, \downarrow }{d\left(L+S\right)} +\left[1-\left(1-\beta \right)\omega \right]I\, \downarrow -\omega \beta I\, \uparrow =\omega \bar{\mu }K\left(1-\beta _{0} \right)e^{-K\left(L+S\right)}\]

where \(I\, \uparrow\) and \(I\, \downarrow\) are the upward and downward diffuse radiative fluxes per unit incident flux, \(K={G\left(\mu \right)\mathord{\left/ {\vphantom {G\left(\mu \right) \mu }} \right.} \mu }\) is the optical depth of direct beam per unit leaf and stem area, \(\mu\) is the cosine of the zenith angle of the incident beam, \(G\left(\mu \right)\) is the relative projected area of leaf and stem elements in the direction \(\cos ^{-1} \mu\), \(\bar{\mu }\) is the average inverse diffuse optical depth per unit leaf and stem area, \(\omega\) is a scattering coefficient, \(\beta\) and \(\beta _{0}\) are upscatter parameters for diffuse and direct beam radiation, respectively, \(L\) is the exposed leaf area index, and \(S\) is the exposed stem area index (section 2.2.1.4). Given the direct beam albedo \(\alpha _{g,\, \Lambda }^{\mu }\) and diffuse albedo \(\alpha _{g,\, \Lambda }\) of the ground (section 2.3.2), these equations are solved to calculate the fluxes, per unit incident flux, absorbed by the vegetation, reflected by the vegetation, and transmitted through the vegetation for direct and diffuse radiation and for visible (\(<\) 0.7\(\mu {\rm m}\)) and near-infrared (\(\geq\) 0.7\(\mu {\rm m}\)) wavebands. The absorbed radiation is partitioned to sunlit and shaded fractions of the canopy. The optical parameters \(G\left(\mu \right)\), \(\bar{\mu }\), \(\omega\), \(\beta\), and \(\beta _{0}\) are calculated based on work in Sellers (1985) as follows.

The relative projected area of leaves and stems in the direction \(\cos ^{-1} \mu\) is

(2.3.3)\[G\left(\mu \right)=\phi _{1} +\phi _{2} \mu\]

where \(\phi _{1} ={\rm 0.5}-0.633\chi _{L} -0.33\chi _{L}^{2}\) and \(\phi _{2} =0.877\left(1-2\phi _{1} \right)\) for \(-0.4\le \chi _{L} \le 0.6\). \(\chi _{L}\) is the departure of leaf angles from a random distribution and equals +1 for horizontal leaves, 0 for random leaves, and –1 for vertical leaves.

The average inverse diffuse optical depth per unit leaf and stem area is

(2.3.4)\[\bar{\mu }=\int _{0}^{1}\frac{\mu '}{G\left(\mu '\right)} d\mu '=\frac{1}{\phi _{2} } \left[1-\frac{\phi _{1} }{\phi _{2} } \ln \left(\frac{\phi _{1} +\phi _{2} }{\phi _{1} } \right)\right]\]

where \(\mu '\) is the direction of the scattered flux.

The optical parameters \(\omega\), \(\beta\), and \(\beta _{0}\), which vary with wavelength (\(\Lambda\) ), are weighted combinations of values for vegetation and snow, using the canopy snow-covered fraction \(f_{can,\, sno}\) (Chapter 2.7). The optical parameters are

(2.3.5)\[\omega _{\Lambda } =\omega _{\Lambda }^{veg} \left(1-f_{can,\, sno} \right)+\omega _{\Lambda }^{sno} f_{can,\, sno}\]

(2.3.6)\[\omega _{\Lambda } \beta _{\Lambda } =\omega _{\Lambda }^{veg} \beta _{\Lambda }^{veg} \left(1-f_{can,\, sno} \right)+\omega _{\Lambda }^{sno} \beta _{\Lambda }^{sno} f_{can,\, sno}\]

(2.3.7)\[\omega _{\Lambda } \beta _{0,\, \Lambda } =\omega _{\Lambda }^{veg} \beta _{0,\, \Lambda }^{veg} \left(1-f_{can,\, sno} \right)+\omega _{\Lambda }^{sno} \beta _{0,\, \Lambda }^{sno} f_{can,\, sno}\]

The snow and vegetation weights are applied to the products \(\omega _{\Lambda } \beta _{\Lambda }\) and \(\omega _{\Lambda } \beta _{0,\, \Lambda }\) because these products are used in the two-stream equations. If there is no snow on the canopy, this reduces to

(2.3.8)\[\omega _{\Lambda } =\omega _{\Lambda }^{veg}\]

(2.3.9)\[\omega _{\Lambda } \beta _{\Lambda } =\omega _{\Lambda }^{veg} \beta _{\Lambda }^{veg}\]

(2.3.10)\[\omega _{\Lambda } \beta _{0,\, \Lambda } =\omega _{\Lambda }^{veg} \beta _{0,\, \Lambda }^{veg} .\]

For vegetation, \(\omega _{\Lambda }^{veg} =\alpha _{\Lambda } +\tau _{\Lambda }\). \(\alpha _{\Lambda }\) is a weighted combination of the leaf and stem reflectances (\(\alpha _{\Lambda }^{leaf},\alpha _{\Lambda }^{stem}\) )

(2.3.11)\[\alpha _{\Lambda } =\max(\alpha _{\Lambda }^{leaf} w_{leaf} +\alpha _{\Lambda }^{stem} w_{stem},1e-6)\]

where \(w_{leaf} =\max({L\mathord{\left/ {\vphantom {L \left(L+S\right)}} \right.} \left(L+S\right)},1e-6)\) and \(w_{stem} =\max({S\mathord{\left/ {\vphantom {S \left(L+S\right)}} \right.} \left(L+S\right)},1e-6)\). \(\tau _{\Lambda }\) is a weighted combination of the leaf and stem transmittances (\(\tau _{\Lambda }^{leaf}, \tau _{\Lambda }^{stem}\))

(2.3.12)\[\tau _{\Lambda } =\max(\tau _{\Lambda }^{leaf} w_{leaf} +\tau _{\Lambda }^{stem} w_{stem},1e-6) .\]

The upscatter for diffuse radiation is

(2.3.13)\[\omega _{\Lambda }^{veg} \beta _{\Lambda }^{veg} =\frac{1}{2} \left[\alpha _{\Lambda } +\tau _{\Lambda } +\left(\alpha _{\Lambda } -\tau _{\Lambda } \right)\cos ^{2} \bar{\theta }\right]\]

where \(\bar{\theta }\) is the mean leaf inclination angle relative to the horizontal plane (i.e., the angle between leaf normal and local vertical) (Sellers (1985)). Here, \(\cos \bar{\theta }\) is approximated by

(2.3.14)\[\cos \bar{\theta }=\frac{1+\chi _{L} }{2}\]

Using this approximation, for vertical leaves (\(\chi _{L} =-1\), \(\bar{\theta }=90^{{\rm o}}\) ), \(\omega _{\Lambda }^{veg} \beta _{\Lambda }^{veg} =0.5\left(\alpha _{\Lambda } +\tau _{\Lambda } \right)\), and for horizontal leaves (\(\chi _{L} =1\), \(\bar{\theta }=0^{{\rm o}}\) ), \(\omega _{\Lambda }^{veg} \beta _{\Lambda }^{veg} =\alpha _{\Lambda }\), which agree with both Dickinson (1983) and Sellers (1985). For random (spherically distributed) leaves (\(\chi _{L} =0\), \(\bar{\theta }=60^{{\rm o}}\) ), the approximation yields \(\omega _{\Lambda }^{veg} \beta _{\Lambda }^{veg} ={5\mathord{\left/ {\vphantom {5 8}} \right.} 8} \alpha _{\Lambda } +{3\mathord{\left/ {\vphantom {3 8}} \right.} 8} \tau _{\Lambda }\) whereas the approximate solution of Dickinson (1983) is \(\omega _{\Lambda }^{veg} \beta _{\Lambda }^{veg} ={2\mathord{\left/ {\vphantom {2 3}} \right.} 3} \alpha _{\Lambda } +{1\mathord{\left/ {\vphantom {1 3}} \right.} 3} \tau _{\Lambda }\). This discrepancy arises from the fact that a spherical leaf angle distribution has a true mean leaf inclination \(\bar{\theta }\approx 57\) (Campbell and Norman 1998) in equation (2.3.13), while \(\bar{\theta }=60\) in equation (2.3.14). The upscatter for direct beam radiation is

(2.3.15)\[\omega _{\Lambda }^{veg} \beta _{0,\, \Lambda }^{veg} =\frac{1+\bar{\mu }K}{\bar{\mu }K} a_{s} \left(\mu \right)_{\Lambda }\]

where the single scattering albedo is

(2.3.16)\[\begin{split}\begin{array}{rcl} {a_{s} \left(\mu \right)_{\Lambda } } & {=} & {\frac{\omega _{\Lambda }^{veg} }{2} \int _{0}^{1}\frac{\mu 'G\left(\mu \right)}{\mu G\left(\mu '\right)+\mu 'G\left(\mu \right)} d\mu '} \\ {} & {=} & {\frac{\omega _{\Lambda }^{veg} }{2} \frac{G\left(\mu \right)}{\max (\mu \phi _{2} +G\left(\mu \right),1e-6)} \left[1-\frac{\mu \phi _{1} }{\max (\mu \phi _{2} +G\left(\mu \right),1e-6)} \ln \left(\frac{\mu \phi _{1} +\max (\mu \phi _{2} +G\left(\mu \right),1e-6)}{\mu \phi _{1} } \right)\right].} \end{array}\end{split}\]

Note here the restriction on \(\mu \phi _{2} +G\left(\mu \right)\). We have seen cases where small values can cause unrealistic single scattering albedo associated with the log calculation, thereby eventually causing a negative soil albedo.

Note that the parameters \(h_{1}\) –\(h_{10}\), \(\sigma\), \(h\), \(s_{1}\), and \(s_{2}\) in the following equations [(2.3.17) - (2.3.30)] are defined in (2.3.31) - (2.3.57).

The upward diffuse fluxes per unit incident direct beam and diffuse flux (i.e., the surface albedos) are

(2.3.17)\[I\, \uparrow _{\Lambda }^{\mu } =\frac{h_{1} }{\sigma } +h_{2} +h_{3}\]

(2.3.18)\[I\, \uparrow _{\Lambda } =h_{7} +h_{8} .\]

The downward diffuse fluxes per unit incident direct beam and diffuse radiation, respectively, are

(2.3.19)\[I\, \downarrow _{\Lambda }^{\mu } =\frac{h_{4} }{\sigma } s_{2} +h_{5} s_{1} +\frac{h_{6} }{s_{1} }\]

(2.3.20)\[I\, \downarrow _{\Lambda } =h_{9} s_{1} +\frac{h_{10} }{s_{1} } .\]

With reference to Figure 2.4.1, the direct beam flux transmitted through the canopy, per unit incident flux, is \(e^{-K\left(L+S\right)}\), and the direct beam and diffuse fluxes absorbed by the vegetation, per unit incident flux, are

(2.3.21)\[\vec{I}_{\Lambda }^{\mu } =1-I\, \uparrow _{\Lambda }^{\mu } -\left(1-\alpha _{g,\, \Lambda } \right)I\, \downarrow _{\Lambda }^{\mu } -\left(1-\alpha _{g,\, \Lambda }^{\mu }\right) s_{2}\]

(2.3.22)\[\vec{I}_{\Lambda } =1-I\, \uparrow _{\Lambda } -\left(1-\alpha _{g,\, \Lambda } \right)I\, \downarrow _{\Lambda } .\]

These fluxes are partitioned to the sunlit and shaded canopy using an analytical solution to the two-stream approximation for sunlit and shaded leaves (Dai et al. 2004), as described by Bonan et al. (2011). The absorption of direct beam radiation by sunlit leaves is

(2.3.23)\[\vec{I}_{sun,\Lambda }^{\mu } =\left(1-\omega _{\Lambda } \right)\left[1-s_{2} +\frac{1}{\bar{\mu }} \left(a_{1} +a_{2} \right)\right]\]

and for shaded leaves is

(2.3.24)\[\vec{I}_{sha,\Lambda }^{\mu } =\vec{I}_{\Lambda }^{\mu } -\vec{I}_{sun,\Lambda }^{\mu }\]

with

(2.3.25)\[a_{1} =\frac{h_{1} }{\sigma } \left[\frac{1-s_{2}^{2} }{2K} \right]+h_{2} \left[\frac{1-s_{2} s_{1} }{K+h} \right]+h_{3} \left[\frac{1-{s_{2} \mathord{\left/ {\vphantom {s_{2} s_{1} }} \right.} s_{1} } }{K-h} \right]\]

(2.3.26)\[a_{2} =\frac{h_{4} }{\sigma } \left[\frac{1-s_{2}^{2} }{2K} \right]+h_{5} \left[\frac{1-s_{2} s_{1} }{K+h} \right]+h_{6} \left[\frac{1-{s_{2} \mathord{\left/ {\vphantom {s_{2} s_{1} }} \right.} s_{1} } }{K-h} \right].\]

For diffuse radiation, the absorbed radiation for sunlit leaves is

(2.3.27)\[\vec{I}_{sun,\Lambda }^{} =\left[\frac{1-\omega _{\Lambda } }{\bar{\mu }} \right]\left(a_{1} +a_{2} \right)\]

and for shaded leaves is

(2.3.28)\[\vec{I}_{sha,\Lambda }^{} =\vec{I}_{\Lambda }^{} -\vec{I}_{sun,\Lambda }^{}\]

with

(2.3.29)\[a_{1} =h_{7} \left[\frac{1-s_{2} s_{1} }{K+h} \right]+h_{8} \left[\frac{1-{s_{2} \mathord{\left/ {\vphantom {s_{2} s_{1} }} \right.} s_{1} } }{K-h} \right]\]

(2.3.30)\[a_{2} =h_{9} \left[\frac{1-s_{2} s_{1} }{K+h} \right]+h_{10} \left[\frac{1-{s_{2} \mathord{\left/ {\vphantom {s_{2} s_{1} }} \right.} s_{1} } }{K-h} \right].\]

The parameters \(h_{1}\) –\(h_{10}\), \(\sigma\), \(h\), \(s_{1}\), and \(s_{2}\) are from Sellers (1985) [note the error in \(h_{4}\) in Sellers (1985)]:

(2.3.31)\[b=1-\omega _{\Lambda } +\omega _{\Lambda } \beta _{\Lambda }\]

(2.3.32)\[c=\omega _{\Lambda } \beta _{\Lambda }\]

(2.3.33)\[d=\omega _{\Lambda } \bar{\mu }K\beta _{0,\, \Lambda }\]

(2.3.34)\[f=\omega _{\Lambda } \bar{\mu }K\left(1-\beta _{0,\, \Lambda } \right)\]

(2.3.35)\[h=\frac{\sqrt{b^{2} -c^{2} } }{\bar{\mu }}\]

(2.3.36)\[\sigma =\left(\bar{\mu }K\right)^{2} +c^{2} -b^{2}\]

(2.3.37)\[u_{1} =b-{c\mathord{\left/ {\vphantom {c \alpha _{g,\, \Lambda }^{\mu } }} \right.} \alpha _{g,\, \Lambda }^{\mu } } {\rm \; or\; }u_{1} =b-{c\mathord{\left/ {\vphantom {c \alpha _{g,\, \Lambda } }} \right.} \alpha _{g,\, \Lambda } }\]

(2.3.38)\[u_{2} =b-c\alpha _{g,\, \Lambda }^{\mu } {\rm \; or\; }u_{2} =b-c\alpha _{g,\, \Lambda }\]

(2.3.39)\[u_{3} =f+c\alpha _{g,\, \Lambda }^{\mu } {\rm \; or\; }u_{3} =f+c\alpha _{g,\, \Lambda }\]

(2.3.40)\[s_{1} =\exp \left\{-\min \left[h\left(L+S\right),40\right]\right\}\]

(2.3.41)\[s_{2} =\exp \left\{-\min \left[K\left(L+S\right),40\right]\right\}\]

(2.3.42)\[p_{1} =b+\bar{\mu }h\]

(2.3.43)\[p_{2} =b-\bar{\mu }h\]

(2.3.44)\[p_{3} =b+\bar{\mu }K\]

(2.3.45)\[p_{4} =b-\bar{\mu }K\]

(2.3.46)\[d_{1} =\frac{p_{1} \left(u_{1} -\bar{\mu }h\right)}{s_{1} } -p_{2} \left(u_{1} +\bar{\mu }h\right)s_{1}\]

(2.3.47)\[d_{2} =\frac{u_{2} +\bar{\mu }h}{s_{1} } -\left(u_{2} -\bar{\mu }h\right)s_{1}\]

(2.3.48)\[h_{1} =-dp_{4} -cf\]

(2.3.49)\[h_{2} =\frac{1}{d_{1} } \left[\left(d-\frac{h_{1} }{\sigma } p_{3} \right)\frac{\left(u_{1} -\bar{\mu }h\right)}{s_{1} } -p_{2} \left(d-c-\frac{h_{1} }{\sigma } \left(u_{1} +\bar{\mu }K\right)\right)s_{2} \right]\]

(2.3.50)\[h_{3} =\frac{-1}{d_{1} } \left[\left(d-\frac{h_{1} }{\sigma } p_{3} \right)\left(u_{1} +\bar{\mu }h\right)s_{1} -p_{1} \left(d-c-\frac{h_{1} }{\sigma } \left(u_{1} +\bar{\mu }K\right)\right)s_{2} \right]\]

(2.3.51)\[h_{4} =-fp_{3} -cd\]

(2.3.52)\[h_{5} =\frac{-1}{d_{2} } \left[\left(\frac{h_{4} \left(u_{2} +\bar{\mu }h\right)}{\sigma s_{1} } \right)+\left(u_{3} -\frac{h_{4} }{\sigma } \left(u_{2} -\bar{\mu }K\right)\right)s_{2} \right]\]

(2.3.53)\[h_{6} =\frac{1}{d_{2} } \left[\frac{h_{4} }{\sigma } \left(u_{2} -\bar{\mu }h\right)s_{1} +\left(u_{3} -\frac{h_{4} }{\sigma } \left(u_{2} -\bar{\mu }K\right)\right)s_{2} \right]\]

(2.3.54)\[h_{7} =\frac{c\left(u_{1} -\bar{\mu }h\right)}{d_{1} s_{1} }\]

(2.3.55)\[h_{8} =\frac{-c\left(u_{1} +\bar{\mu }h\right)s_{1} }{d_{1} }\]

(2.3.56)\[h_{9} =\frac{u_{2} +\bar{\mu }h}{d_{2} s_{1} }\]

(2.3.57)\[h_{10} =\frac{-s_{1} \left(u_{2} -\bar{\mu }h\right)}{d_{2} } .\]

Plant functional type optical properties (Table 2.3.1) for \(\chi _{L}\), \(\alpha _{vis}^{leaf}\), \(\alpha _{nir}^{leaf}\), \(\alpha _{vis}^{stem}\), \(\alpha _{nir}^{stem}\), \(\tau _{vis}^{leaf}\), and \(\tau _{nir}^{leaf}\) are from Majasalmi and Bright (2019). One of these parameter values was tuned based on a Perturbed Parameter Ensemble. Specifically, \(\chi _{L}\) for BET Tropical was changed from 0.32 to 0.45. \(\tau _{vis}^{stem}\) and \(\tau _{nir}^{stem}\) for trees and shrubs are from Dorman and Sellers (1989). \(\tau _{vis}^{stem}\) and \(\tau _{nir}^{stem}\) for grassland and crops were derived from full optical range spectra of measured optical properties (Asner et al. 1998). Optical properties for intercepted snow (Table 2.3.2) are from Sellers et al. (1986).

Table 2.3.1 Plant functional type optical properties

Plant Functional Type	\(\chi _{L}\)	\(\alpha _{vis}^{leaf}\)	\(\alpha _{nir}^{leaf}\)	\(\alpha _{vis}^{stem}\)	\(\alpha _{nir}^{stem}\)	\(\tau _{vis}^{leaf}\)	\(\tau _{nir}^{leaf}\)	\(\tau _{vis}^{stem}\)	\(\tau _{nir}^{stem}\)
NET Temperate	0.01	0.09	0.41	0.12	0.36	0.04	0.32	0.001	0.001
NET Boreal	0.01	0.09	0.41	0.12	0.36	0.04	0.32	0.001	0.001
NDT Boreal	0.01	0.08	0.39	0.12	0.36	0.06	0.42	0.001	0.001
BET Tropical	0.45	0.11	0.46	0.21	0.49	0.06	0.33	0.001	0.001
BET temperate	0.32	0.11	0.46	0.21	0.49	0.06	0.33	0.001	0.001
BDT tropical	0.20	0.08	0.41	0.21	0.49	0.06	0.43	0.001	0.001
BDT temperate	0.59	0.08	0.41	0.21	0.49	0.06	0.43	0.001	0.001
BDT boreal	0.59	0.08	0.41	0.21	0.49	0.06	0.43	0.001	0.001
BES temperate	0.32	0.11	0.46	0.21	0.49	0.06	0.33	0.001	0.001
BDS temperate	0.59	0.08	0.41	0.21	0.49	0.06	0.43	0.001	0.001
BDS boreal	0.59	0.08	0.41	0.21	0.49	0.06	0.43	0.001	0.001
C3 arctic grass	-0.23	0.05	0.28	0.31	0.53	0.05	0.40	0.120	0.250
C3 grass	-0.23	0.05	0.28	0.31	0.53	0.05	0.40	0.120	0.250
C4 grass	-0.23	0.05	0.28	0.31	0.53	0.05	0.40	0.120	0.250
C3 Crop	0.30	0.08	0.42	0.31	0.53	0.05	0.40	0.120	0.250
Temp Corn	0.30	0.08	0.42	0.31	0.53	0.05	0.40	0.120	0.250
Spring Wheat	0.30	0.08	0.42	0.31	0.53	0.05	0.40	0.120	0.250
Temp Soybean	0.30	0.08	0.42	0.31	0.53	0.05	0.40	0.120	0.250
Cotton	0.30	0.08	0.42	0.31	0.53	0.05	0.40	0.120	0.250
Rice	0.30	0.08	0.42	0.31	0.53	0.05	0.40	0.120	0.250
Sugarcane	0.30	0.08	0.42	0.31	0.53	0.05	0.40	0.120	0.250
Tropical Corn	0.30	0.08	0.42	0.31	0.53	0.05	0.40	0.120	0.250
Tropical Soybean	0.30	0.08	0.42	0.31	0.53	0.05	0.40	0.120	0.250
Miscanthus	0.30	0.08	0.42	0.31	0.53	0.05	0.40	0.120	0.250
Switchgrass	0.30	0.08	0.42	0.31	0.53	0.05	0.40	0.120	0.250

Table 2.3.2 Intercepted snow optical properties

Parameter	vis	nir
\(\omega ^{sno}\)	0.8	0.4
\(\beta ^{sno}\)	0.5	0.5
\(\beta _{0}^{sno}\)	0.5	0.5

2.3.2. Ground Albedos

The overall direct beam \(\alpha _{g,\, \Lambda }^{\mu }\) and diffuse \(\alpha _{g,\, \Lambda }\) ground albedos are weighted combinations of "soil" and snow albedos

(2.3.58)\[\alpha _{g,\, \Lambda }^{\mu } =\alpha _{soi,\, \Lambda }^{\mu } \left(1-f_{sno} \right)+\alpha _{sno,\, \Lambda }^{\mu } f_{sno}\]

(2.3.59)\[\alpha _{g,\, \Lambda } =\alpha _{soi,\, \Lambda } \left(1-f_{sno} \right)+\alpha _{sno,\, \Lambda } f_{sno}\]

where \(f_{sno}\) is the fraction of the ground covered with snow (section 2.8.1).

\(\alpha _{soi,\, \Lambda }^{\mu }\) and \(\alpha _{soi,\, \Lambda }\) vary with glacier, lake, and soil surfaces. Glacier albedos were originally from Paterson (1994) but were adjusted to the following values agreed upon by the Land Ice Working Group leadership team

\[\alpha _{soi,\, vis}^{\mu } =\alpha _{soi,\, vis} =0.5\]

\[\alpha _{soi,\, nir}^{\mu } =\alpha _{soi,\, nir} =0.3.\]

Note that these are read in from the clm_inparm namelist group (variable albice).

Unfrozen lake direct albedos depend on the cosine of the solar zenith angle \(\mu\)

(2.3.60)\[\alpha _{soi,\, \Lambda }^{\mu } =0.05\left(\max (0.001,\mu +0.15\right))^{-1} .\]

For diffuse radiation the expression in equation (2.3.60) is integrated over the full sky to yield \(\alpha _{soi,\, \Lambda } =0.10\).

For frozen lakes without resolved snow layers (\(snl=0\)), the albedo at cold temperatures \(a_0\) is 0.60 for visible and 0.40 for near infrared radiation. As the temperature at the ice surface, \({T}_{g}\), approaches freezing [ \({T}_{f}\) (K) (Table 2.2.7)], the albedo is relaxed towards 0.10 based on Mironov et al. (2010):

(2.3.61)\[\alpha _{soi,\, \Lambda }^{\mu }=a_{0} \left(1-x\right)+\max(0.10,\alpha _{soi,\, \Lambda }^{\mu })x,x=\exp \left(-95\frac{T_{f} -T_{g} }{T_{f} } \right)\]

(2.3.62)\[\alpha _{soi,\, \Lambda }=a_{0} \left(1-x\right)+0.10x,x=\exp \left(-95\frac{T_{f} -T_{g} }{T_{f} } \right)\]

where \(\alpha _{soi,\, \Lambda }^{\mu }\) is calculated from (2.3.60) for direct radiation.

For frozen lakes with resolved snow layers, the reflectance of the ice surface is fixed at \(a_0\), and the snow reflectance is calculated as over non-vegetated surfaces (Chapter 2.3). These two reflectances are combined to obtain the snow-fraction-weighted albedo as in over non-vegetated surfaces (Chapter 2.3).

As in NCAR LSM (Bonan 1996), soil albedos vary with color class

(2.3.63)\[\alpha _{soi,\, \Lambda }^{\mu } =\alpha _{soi,\, \Lambda } =\left(\alpha _{sat,\, \Lambda } +\Delta \right)\le \alpha _{dry,\, \Lambda }\]

where \(\Delta\) depends on the volumetric water content of the first soil layer \(\theta _{1}\) (section 2.7.3) as \(\Delta =0.11-0.40\theta _{1} >0\), and \(\alpha _{sat,\, \Lambda }\) and \(\alpha _{dry,\, \Lambda }\) are albedos for saturated and dry soil color classes (Table 2.3.3).

Albedos associated with CLM soil colors were updated for CLM4 using the procedure described below (Lawrence et al. (2011)) but have not been updated since then. CLM soil colors are prescribed so that they best reproduce observed MODIS local solar noon surface albedo values at the CLM grid cell following the methods of Lawrence and Chase (2007). The soil colors are fitted over the range of 20 soil classes shown in Table 2.3.3 and compared to the MODIS monthly local solar noon all-sky surface albedo as described in Strahler et al. (1999) and Schaaf et al. (2002). The CLM two-stream radiation model was used to calculate the model equivalent surface albedo using climatological monthly soil moisture along with the vegetation parameters of PFT fraction, LAI, and SAI. The soil color that produced the closest all-sky albedo in the two-stream radiation model was selected as the best fit for the month. The fitted monthly soil colors were averaged over all snow-free months to specify a representative soil color for the grid cell. In cases where there was no snow-free surface albedo for the year, the soil color derived from snow-affected albedo was used to give a representative soil color that included the effects of the minimum permanent snow cover.

Table 2.3.3 Dry and saturated soil albedos

	Dry		Saturated			Dry		Saturated
Color Class	vis	nir	vis	nir	Color Class	vis	nir	vis	nir
1	0.36	0.61	0.25	0.50	11	0.24	0.37	0.13	0.26
2	0.34	0.57	0.23	0.46	12	0.23	0.35	0.12	0.24
3	0.32	0.53	0.21	0.42	13	0.22	0.33	0.11	0.22
4	0.31	0.51	0.20	0.40	14	0.20	0.31	0.10	0.20
5	0.30	0.49	0.19	0.38	15	0.18	0.29	0.09	0.18
6	0.29	0.48	0.18	0.36	16	0.16	0.27	0.08	0.16
7	0.28	0.45	0.17	0.34	17	0.14	0.25	0.07	0.14
8	0.27	0.43	0.16	0.32	18	0.12	0.23	0.06	0.12
9	0.26	0.41	0.15	0.30	19	0.10	0.21	0.05	0.10
10	0.25	0.39	0.14	0.28	20	0.08	0.16	0.04	0.08

2.3.2.1. Snow Albedo

Snow albedo and solar absorption within each snow layer are simulated with the Snow, Ice, and Aerosol Radiative Model (SNICAR-ADv3; Flanner et al., 2021), which incorporates an adding-doubling radiative transfer solution (Dang et al., 2019; Briegleb and Light, 2007; Note that the original SNICAR uses the tri-diagonal matrix two-stream radiative transfer solution from Toon et al., 1989). Snow albedo and the vertical light absorption profile depend on solar zenith angle, albedo of the substrate underlying snow, mass concentrations of atmospheric-deposited aerosols (black carbon, mineral dust, and organic carbon), snow layer thickness, snow density, and ice effective grain size (\(r_{e}\)) that is simulated with a snow aging routine described in section 2.3.2.3. Representation of impurity mass concentrations within the snowpack is described in section 2.8.4. Implementation of SNICAR-ADv3 in CTSM/CLM is described by He et al. (2024), while the original SNICAR implementation is described by Flanner and Zender (2005), and Flanner et al. (2007).

The adding-doubling two-stream solution requires the following bulk optical properties for each snow layer and spectral band: extinction optical depth (\(\tau\)), single-scatter albedo (\(\omega\)), and scattering asymmetry parameter (g). The snow layers used for radiative calculations are identical to snow layers applied elsewhere in CLM, except for the case when snow mass is greater than zero but no snow layers exist. When this occurs, a single radiative layer is specified to have the column snow mass and an effective grain size of freshly-fallen snow (section 2.3.2.3). For aerosols externally mixed with snow grains, the bulk optical properties are weighted functions of each constituent k, computed for each snow layer and spectral band as

(2.3.64)\[\tau =\sum _{1}^{k}\tau _{k}\]

(2.3.65)\[\omega =\frac{\sum _{1}^{k}\omega _{k} \tau _{k} }{\sum _{1}^{k}\tau _{k} }\]

(2.3.66)\[g=\frac{\sum _{1}^{k}g_{k} \omega _{k} \tau _{k} }{\sum _{1}^{k}\omega _{k} \tau _{k} }\]

For each constituent (ice, two black carbon species, two organic carbon species, and four dust species), \(\omega\), g, and the mass extinction cross-section \(\psi\) (m² kg^-1) are pre-computed offline with Mie Theory by assuming lognormal size distribution for each constituent and applying the computational technique from Bohren and Huffman (1983). The extinction optical depth for each constituent depends on its mass extinction cross-section and layer mass, \(w _{k}\) (kgm^-1) as

(2.3.67)\[\tau _{k} =\psi _{k} w_{k}\]

Compared to the original SNICAR implementation, the SNICAR-ADv3 implementation includes two new features: nonspherical snow grain shape (He et al., 2017) and aerosol-snow internal mixing for black carbon (BC) and dust (He et al., 2017; He et al., 2019). The original SNICAR assumes spherical snow grains, which however may not be a realistic representation since nonspherical snow grains are ubiquitous in reality. Thus, in SNICAR-ADv3, four typical snow grain shapes representative of real‐world observations, including sphere, spheroid, hexagonal plate/column, and fractal snowflake are used. The He et al. (2017) parameterizations are used to quantify snow grain nonsphericity impacts on snow asymmetry factor (g) before the calculation of bulk optical properties (equation (2.3.66)). Snow extinction optical depth (\(\tau\)) and single-scattering albedo (\(\omega\)) are not modified. The default snow grain shape is set to hexagonal plate/column (Flanner et al., 2021). The snow grain shape can be controlled via namelist option (He et al., 2024). For BC internal mixing with snow grain, the original SNICAR assumes BC‐snow external mixing, with hydrophilic BC treated as coated BC. However, in reality BC can also be internally mixed with snow grains. The He et al. (2017) parameterizations are used to quantify BC-snow internal mixing effects on snow single-scattering albedo (\(\omega\)), with no changes in snow extinction optical depth (\(\tau\)) and asymmetry factor (g) due to neglibile impacts. The hydrophilic BC is no longer treated as coated BC because sulfate coating on the BC particle surface is dissolvable into water during wet deposition. The hydrophilic BC is hence treated as internally mixed with snow grains. For dust-snow internal mixing, the He et al. (2019) parameterizations are used to quantify its impact on snow single-scattering albedo (\(\omega\)), with no changes in snow extinction optical depth (\(\tau\)) and asymmetry factor (g). Currently, dust-snow internal mixing (if activated) applies to all dust size bins. Both BC-snow and dust-snow internal/external mixing can be controlled by namelist options (He et al., 2024). By default, BC and dust are externally mixed with snow grains, and it is recommended not to activate both internal mixing together which has not been fully tested.

The original SNICAR two-stream solution (Toon et al., 1989) applies a tri-diagonal matrix solution to produce upward and downward radiative fluxes at each layer interface, from which net radiation, layer absorption, and surface albedo are easily derived. The SNICAR-ADv3 instead uses a more accurate adding-doubling radiative transfer solution (Dang et al., 2019; Briegleb and Light, 2007). Solar fluxes are computed in either five spectral bands (default) listed in Table 2.3.4 or 480 hyperspectral bands (from 200 nm to 5000 nm with 10-nm spectral resolution; Flanner et al., 2021). Because snow albedo varies strongly across the solar spectrum, it was determined that at least four bands were needed to accurately represent the near-infrared (NIR) characteristics of snow, whereas only one band was needed for the visible spectrum, which is why the five bands are used by default to achieve a balance between computational time and accuracy. Boundaries of the four NIR bands were selected to capture broad radiative features and maximize accuracy and computational efficiency. The 480-band capability can be used via namelist option (He et al., 2024). We partition NIR (0.7-5.0 \(\mu\) m) surface downwelling flux from CLM according to the weights from six typical atmospheric conditions (Flanner et al., 2021), including Mid‐latitude winter, Mid‐latitude summer, Sub‐Arctic winter, Sub‐Arctic summer, Summit (Greenland), and High mountain, which are unique for diffuse and direct incident flux and can be selected via namelist options (He et al., 2024). By default, the mid‐latitude winter downward solar spectrum (Table 2.3.4) is used. These prescribed weights were determined with offline hyperspectral atmospheric radiative transfer calculations. The radiative transfer solution includes intermediate terms that allow for easy interchange of two-stream techniques. We apply the delta-Eddington solution to the layer bulk optical properties following Briegleb and Light (2007). Specifically, Delta scalings are applied to \(\tau\), \(\omega\), and \(g\) in all spectral bands, producing effective values (denoted with \(*\)) that are applied in the two-stream solution.

(2.3.68)\[\tau ^{*} =\left(1-\omega g^{2} \right)\tau\]

(2.3.69)\[\omega ^{*} =\frac{\left(1-g^{2} \right)\omega }{1-g^{2} \omega }\]

(2.3.70)\[g^{*} =\frac{g-g^{2}}{1-g^{2}}\]

Table 2.3.4 Spectral bands and weights used for snow radiative transfer

Spectral band	Direct-beam weight	Diffuse weight
Band 1: 0.3-0.7\(\mu\)m (visible)	(1.0)	(1.0)
Band 2: 0.7-1.0\(\mu\)m (near-IR)	0.494	0.634
Band 3: 1.0-1.2\(\mu\)m (near-IR)	0.180	0.186
Band 4: 1.2-1.5\(\mu\)m (near-IR)	0.123	0.094
Band 5: 1.5-5.0\(\mu\)m (near-IR)	0.203	0.086

Under direct-beam conditions, the two-stream approximations become poor for large solar zenith angle, which is mostly contributed by the errors of near-IR band calculations, especially for optically thick snowpacks. To improve the performance of two-stream algorithms, the Dang et al. (2019) parameterization that corrects the underestimated near-IR snow albedo at large solar zenith angles (>75deg) is used for NIR bands (He et al., 2024).

Soil albedo (or underlying substrate albedo), which is defined for visible and NIR bands, is a required boundary condition for the snow radiative transfer calculation. Currently, the bulk NIR soil albedo is applied to all four NIR snow bands. With ground albedo as a lower boundary condition, SNICAR-ADv3 simulates solar absorption in all snow layers as well as the underlying soil or ground. With a thin snowpack, penetrating solar radiation to the underlying soil can be quite large and heat cannot be released from the soil to the atmosphere in this situation. Thus, if the snowpack has total snow depth less than 0.1 m (\(z_{sno} < 0.1\)) and there are no explicit snow layers, the solar radiation is absorbed by the top soil layer. If there is a single snow layer, the solar radiation is absorbed in that layer. If there is more than a single snow layer, 75% of the solar radiation is absorbed in the top snow layer, and 25% is absorbed in the next lowest snow layer. This prevents unrealistic soil warming within a single timestep.

The radiative transfer calculation is performed twice for each column containing a mass of snow greater than \(1 \times 10^{-30}\) kgm^-2 (excluding lake and urban columns); once each for direct-beam and diffuse incident flux. Absorption in each layer \(i\) of pure snow is initially recorded as absorbed flux per unit incident flux on the ground (\(S_{sno,\, i}\) ), as albedos must be calculated for the next timestep with unknown incident flux. The snow absorption fluxes that are used for column temperature calculations are

(2.3.71)\[S_{g,\, i} =S_{sno,\, i} \left(1-\alpha _{sno} \right)\]

This weighting is performed for direct-beam and diffuse, visible and NIR fluxes. After the ground-incident fluxes (transmitted through the vegetation canopy) have been calculated for the current time step (sections 2.3.1 and 2.4.1), the layer absorption factors (\(S_{g,\, i}\)) are multiplied by the ground-incident fluxes to produce solar absorption (W m^-2) in each snow layer and the underlying ground.

2.3.2.2. Snowpack Optical Properties

Ice optical properties for the five or 480 spectral bands are derived offline based on Mie calculations and stored in a namelist-defined lookup table for online retrieval (He et al., 2024). The ice refractive index used in Mie calculation follows the updates in (Flanner et al., 2021), which is a compilation of the Picard et al. (2016) and Warren and Brandt (2008) datasets. Mie properties are first computed at fine spectral resolution (480 bands from 200 nm to 5000 nm with 10-nm spectral resolution), and are then weighted into the five bands according to incident solar flux, \(I^{\downarrow } (\lambda )\), under six different atmospheric conditions (section 2.3.2.1). For example, the broadband mass-extinction cross section (\(\bar{\psi }\)) over wavelength interval \(\lambda _{1}\) to \(\lambda _{2}\) is

(2.3.72)\[\bar{\psi }=\frac{\int _{\lambda _{1} }^{\lambda _{2} }\psi \left(\lambda \right) I^{\downarrow } \left(\lambda \right){\rm d}\lambda }{\int _{\lambda _{1} }^{\lambda _{2} }I^{\downarrow } \left(\lambda \right){\rm d}\lambda }\]

Broadband single-scatter albedo (\(\bar{\omega }\)) is additionally weighted by the diffuse albedo for a semi-infinite snowpack (\(\alpha _{sno}\))

(2.3.73)\[\bar{\omega }=\frac{\int _{\lambda _{1} }^{\lambda _{2} }\omega (\lambda )I^{\downarrow } ( \lambda )\alpha _{sno} (\lambda ){\rm d}\lambda }{\int _{\lambda _{1} }^{\lambda _{2} }I^{\downarrow } ( \lambda )\alpha _{sno} (\lambda ){\rm d}\lambda }\]

Inclusion of this additional albedo weight was found to improve accuracy of the five-band albedo solutions (relative to 480-band solutions) because of the strong dependence of optically-thick snowpack albedo on ice grain single-scatter albedo (Flanner et al., 2007). The lookup tables contain optical properties for lognormal distributions of ice particles over the range of effective radii: 30\(\mu\)m \(< r _{e} < \text{1500} \mu \text{m}\), at 1 \(\mu\) m resolution. Single-scatter albedos for the end-members of this size range are listed in Table 2.3.5.

Optical properties for black carbon, organic carbon, and mineral dust are described in Flanner et al. (2021). Three types of dust can be selected via namelist options (He et al., 2024), including Saharan dust (default), Colorado dust, and Greenland dust, due to their substantially different optical properties. Single-scatter albedo, mass extinction cross-section, and asymmetry parameter values for all snowpack species under diffuse radiation with the default atmospheric condition (mid-latitude winter), in the five spectral bands used, are listed in Table 2.3.5, Table 2.3.6, and Table 2.3.7. These properties were also derived with Mie Theory, using various published sources of indices of refraction and assumptions about particle size distribution (Flanner et al., 2021). Weighting into the five CLM spectral bands was determined using the 480-band values with incident solar flux under six different atmospheric conditions (section 2.3.2.1), as in equation (2.3.72).

Table 2.3.5 Single-scatter albedo values used for snowpack impurities and ice

Species	Band 1	Band 2	Band 3	Band 4	Band 5
Hydrophilic black carbon	0.366	0.302	0.252	0.217	0.152
Hydrophobic black carbon	0.366	0.302	0.252	0.217	0.152
Hydrophilic organic carbon	0.772	0.990	0.987	0.983	0.971
Hydrophobic organic carbon	0.772	0.990	0.987	0.983	0.971
Dust 1	0.945	0.991	0.992	0.992	0.983
Dust 2	0.870	0.976	0.989	0.992	0.991
Dust 3	0.802	0.948	0.965	0.974	0.984
Dust 4	0.730	0.913	0.943	0.954	0.965
Ice (\(r _{e}\) = 30 \(\mu\) m)	0.9999	0.9999	0.9993	0.9954	0.9510
Ice (\(r _{e}\) = 1500 \(\mu\) m)	0.9998	0.9963	0.9678	0.8735	0.5492

Table 2.3.6 Mass extinction values (m² kg^-1) used for snowpack impurities and ice

Species	Band 1	Band 2	Band 3	Band 4	Band 5
Hydrophilic black carbon	12389	7971	5744	4654	3155
Hydrophobic black carbon	12389	7971	5744	4654	3155
Hydrophilic organic carbon	4933	1390	628	375	143
Hydrophobic organic carbon	4933	1390	628	375	143
Dust 1	2543	2242	1469	1013	458
Dust 2	803	950	1144	1205	1000
Dust 3	369	399	378	380	488
Dust 4	188	193	197	197	212
Ice (\(r _{e}\) = 30 \(\mu\) m)	55.7	56.1	56.4	56.7	57.2
Ice (\(r _{e}\) = 1500 \(\mu\) m)	1.09	1.09	1.09	1.09	1.09

Table 2.3.7 Asymmetry scattering parameters used for snowpack impurities and ice.

Species	Band 1	Band 2	Band 3	Band 4	Band 5
Hydrophilic black carbon	0.44	0.34	0.28	0.24	0.18
Hydrophobic black carbon	0.44	0.34	0.28	0.24	0.18
Hydrophilic organic carbon	0.58	0.47	0.39	0.34	0.25
Hydrophobic organic carbon	0.58	0.47	0.39	0.34	0.25
Dust 1	0.71	0.72	0.67	0.62	0.48
Dust 2	0.73	0.66	0.71	0.74	0.73
Dust 3	0.82	0.76	0.68	0.63	0.67
Dust 4	0.87	0.80	0.78	0.76	0.73
Ice (\(r _{e}\) = 30\(\mu\)m)	0.88	0.88	0.88	0.88	0.89
Ice (\(r _{e}\) = 1500\(\mu\)m)	0.89	0.90	0.90	0.92	0.97

2.3.2.3. Snow Aging

Snow aging is represented as evolution of the ice effective grain size (\(r_{e}\)). Previous studies have shown that use of spheres which conserve the surface area-to-volume ratio (or specific surface area) of ice media composed of more complex shapes produces relatively small errors in simulated hemispheric fluxes (e.g., Grenfell and Warren 1999). Effective radius is the surface area-weighted mean radius of an ensemble of spherical particles and is directly related to specific surface area (SSA) as \(r_{e} ={3\mathord{\left/ {\vphantom {3 \left(\rho _{ice} SSA\right)}} \right.} \left(\rho _{ice} SSA\right)}\), where \(\rho_{ice}\) is the density of ice. Hence, \(r_{e}\) is a simple and practical metric for relating the snowpack microphysical state to dry snow radiative characteristics.

Wet snow processes can also drive rapid changes in albedo. The presence of liquid water induces rapid coarsening of the surrounding ice grains (e.g., Brun 1989), and liquid water tends to refreeze into large ice clumps that darken the bulk snowpack. The presence of small liquid drops, by itself, does not significantly darken snowpack, as ice and water have very similar indices of refraction throughout the solar spectrum. Pooled or ponded water, however, can significantly darken snowpack by greatly reducing the number of refraction events per unit mass. This influence is not currently accounted for.

The net change in effective grain size occurring each time step is represented in each snow layer as a summation of changes caused by dry snow metamorphism (\(dr_{e,dry}\)), liquid water-induced metamorphism (\(dr_{e,wet}\)), refreezing of liquid water, and addition of freshly-fallen snow. The mass of each snow layer is partitioned into fractions of snow carrying over from the previous time step (\(f_{old}\)), freshly-fallen snow (\(f_{new}\)), and refrozen liquid water (\(f_{rfz}\)), such that snow \(r_{e}\) is updated each time step t as

(2.3.74)\[r_{e} \left(t\right)=\left[r_{e} \left(t-1\right)+dr_{e,\, dry} +dr_{e,\, wet} \right]f_{old} +r_{e,\, 0} f_{new} +r_{e,\, rfz} f_{rfrz}\]

Here, the effective radius of freshly-fallen snow (\(r_{e,0}\)) is based on a simple linear temperature-relationship. Below -30 degrees Celsius, a minimum value is enforced of 54.5 \(\mu\) m (corresponding to a specific surface area of 60 m² kg^-1). Above 0 degrees Celsius, a maximum value is enforced of 204.5 \(\mu\) m. Between -30 and 0 a linear ramp is used.

The effective radius of refrozen liquid water (\(r_{e,rfz}\)) is set to 1000\(\mu\) m.

Dry snow aging is based on a microphysical model described by Flanner and Zender (2006). This model simulates diffusive vapor flux amongst collections of ice crystals with various size and inter-particle spacing. Specific surface area and effective radius are prognosed for any combination of snow temperature, temperature gradient, density, and initial size distribution. The combination of warm snow, large temperature gradient, and low density produces the most rapid snow aging, whereas aging proceeds slowly in cold snow, regardless of temperature gradient and density. Because this model is currently too computationally expensive for inclusion in climate models, we fit parametric curves to model output over a wide range of snow conditions and apply these parameters in CLM. The functional form of the parametric equation is

(2.3.75)\[\frac{dr_{e,\, dry} }{dt} =\left(\frac{dr_{e} }{dt} \right)_{0} \left(\frac{\eta }{\left(r_{e} -r_{e,\, 0} \right)+\eta } \right)^{{1\mathord{\left/ {\vphantom {1 \kappa }} \right.} \kappa } }\]

The parameters \({(\frac{dr_{e}}{dt}})_{0}\), \(\eta\), and \(\kappa\) are retrieved interactively from a lookup table with dimensions corresponding to snow temperature, temperature gradient, and density. The domain covered by this lookup table includes temperature ranging from 223 to 273 K, temperature gradient ranging from 0 to 300 K m^-1, and density ranging from 50 to 400 kg m^-3. Temperature gradient is calculated at the midpoint of each snow layer n, using mid-layer temperatures (\(T_{n}\)) and snow layer thicknesses (\(dz_{n}\)), as

(2.3.76)\[\left(\frac{dT}{dz} \right)_{n} =\frac{1}{dz_{n} } abs\left[\frac{T_{n-1} dz_{n} +T_{n} dz_{n-1} }{dz_{n} +dz_{n-1} } +\frac{T_{n+1} dz_{n} +T_{n} dz_{n+1} }{dz_{n} +dz_{n+1} } \right]\]

For the bottom snow layer (\(n=0\)), \(T_{n+1}\) is taken as the temperature of the top soil layer, and for the top snow layer it is assumed that \(T_{n-1}\) = \(T_{n}\).

The contribution of liquid water to enhanced metamorphism is based on parametric equations published by Brun (1989), who measured grain growth rates under different liquid water contents. This relationship, expressed in terms of \(r_{e} (\mu \text{m})\) and subtracting an offset due to dry aging, depends on the mass liquid water fraction \(f_{liq}\) as

(2.3.77)\[\frac{dr_{e} }{dt} =\frac{10^{18} C_{1} f_{liq} ^{3} }{4\pi r_{e} ^{2} }\]

The constant C₁ is 4.22\(\times\)10^-13, and: \(f_{liq} =w_{liq} /(w_{liq} +w_{ice} )\)(Chapter 2.8).

In cases where snow mass is greater than zero, but a snow layer has not yet been defined, \(r_{e}\) is set to \(r_{e,0}\). When snow layers are combined or divided, \(r_{e}\) is calculated as a mass-weighted mean of the two layers, following computations of other state variables (section 2.8.7). Finally, the allowable range of \(r_{e}\), corresponding to the range over which Mie optical properties have been defined, is 30-1500\(\mu\) m.

2.3.3. Solar Zenith Angle

The CLM uses the same formulation for solar zenith angle as the Community Atmosphere Model. The cosine of the solar zenith angle \(\mu\) is

(2.3.78)\[\mu =\sin \phi \sin \delta -\cos \phi \cos \delta \cos h\]

where \(h\) is the solar hour angle (radians) (24 hour periodicity), \(\delta\) is the solar declination angle (radians), and \(\phi\) is latitude (radians) (positive in Northern Hemisphere). The solar hour angle \(h\) (radians) is

(2.3.79)\[h=2\pi d+\theta\]

where \(d\) is calendar day (\(d=0.0\) at 0Z on January 1), and \(\theta\) is longitude (radians) (positive east of the Greenwich meridian).

The solar declination angle \(\delta\) is calculated as in Berger (1978a,b) and is valid for one million years past or hence, relative to 1950 A.D. The orbital parameters may be specified directly or the orbital parameters are calculated for the desired year. The required orbital parameters to be input by the user are the obliquity of the Earth \(\varepsilon\) (degrees, \(-90^{\circ } <\varepsilon <90^{\circ }\) ), Earth's eccentricity \(e\) (\(0.0<e<0.1\)), and the longitude of the perihelion relative to the moving vernal equinox \(\tilde{\omega }\) (\(0^{\circ } <\tilde{\omega }<360^{\circ }\) ) (unadjusted for the apparent orbit of the Sun around the Earth (Berger et al. 1993)). The solar declination \(\delta\) (radians) is

(2.3.80)\[\delta =\sin ^{-1} \left[\sin \left(\varepsilon \right)\sin \left(\lambda \right)\right]\]

where \(\varepsilon\) is Earth's obliquity and \(\lambda\) is the true longitude of the Earth.

The obliquity of the Earth \(\varepsilon\) (degrees) is

(2.3.81)\[\varepsilon =\varepsilon^{*}+\sum _{i=1}^{i=47}A_{i} \cos \left(f_{i} t+\delta _{i} \right)\]

where \(\varepsilon^{*}\) is a constant of integration (Table 2.3.8), \(A_{i}\), \(f_{i}\), and \(\delta _{i}\) are amplitude, mean rate, and phase terms in the cosine series expansion (Berger (1978a,b), and \(t=t_{0} -1950\) where \(t_{0}\) is the year. The series expansion terms are not shown here but can be found in the source code file shr_orb_mod.F90.

The true longitude of the Earth \(\lambda\) (radians) is counted counterclockwise from the vernal equinox (\(\lambda =0\) at the vernal equinox)

(2.3.82)\[\lambda =\lambda _{m} +\left(2e-\frac{1}{4} e^{3} \right)\sin \left(\lambda _{m} -\tilde{\omega }\right)+\frac{5}{4} e^{2} \sin 2\left(\lambda _{m} -\tilde{\omega }\right)+\frac{13}{12} e^{3} \sin 3\left(\lambda _{m} -\tilde{\omega }\right)\]

where \(\lambda _{m}\) is the mean longitude of the Earth at the vernal equinox, \(e\) is Earth's eccentricity, and \(\tilde{\omega }\) is the longitude of the perihelion relative to the moving vernal equinox. The mean longitude \(\lambda _{m}\) is

(2.3.83)\[\lambda _{m} =\lambda _{m0} +\frac{2\pi \left(d-d_{ve} \right)}{365}\]

where \(d_{ve} =80.5\) is the calendar day at vernal equinox (March 21 at noon), and

(2.3.84)\[\lambda _{m0} =2\left[\left(\frac{1}{2} e+\frac{1}{8} e^{3} \right)\left(1+\beta \right)\sin \tilde{\omega }-\frac{1}{4} e^{2} \left(\frac{1}{2} +\beta \right)\sin 2\tilde{\omega }+\frac{1}{8} e^{3} \left(\frac{1}{3} +\beta \right)\sin 3\tilde{\omega }\right]\]

where \(\beta =\sqrt{1-e^{2} }\). Earth's eccentricity \(e\) is

(2.3.85)\[e=\sqrt{\left(e^{\cos } \right)^{2} +\left(e^{\sin } \right)^{2} }\]

where

(2.3.86)\[\begin{split}\begin{array}{l} {e^{\cos } =\sum _{j=1}^{19}M_{j} \cos \left(g_{j} t+B_{j} \right) ,} \\ {e^{\sin } =\sum _{j=1}^{19}M_{j} \sin \left(g_{j} t+B_{j} \right) } \end{array}\end{split}\]

are the cosine and sine series expansions for \(e\), and \(M_{j}\), \(g_{j}\), and \(B_{j}\) are amplitude, mean rate, and phase terms in the series expansions (Berger (1978a,b)) (series expansion terms not shown but can be found in the source code file shr_orb_mod.F90). The longitude of the perihelion relative to the moving vernal equinox \(\tilde{\omega }\) (degrees) is

(2.3.87)\[\tilde{\omega }=\Pi \frac{180}{\pi } +\psi\]

where \(\Pi\) is the longitude of the perihelion measured from the reference vernal equinox (i.e., the vernal equinox at 1950 A.D.) and describes the absolute motion of the perihelion relative to the fixed stars, and \(\psi\) is the annual general precession in longitude and describes the absolute motion of the vernal equinox along Earth's orbit relative to the fixed stars. The general precession \(\psi\) (degrees) is

(2.3.88)\[\psi =\frac{\tilde{\psi }t}{3600} +\zeta +\sum _{i=1}^{78}F_{i} \sin \left(f_{i} ^{{'} } t+\delta _{i} ^{{'} } \right)\]

where \(\tilde{\psi }\) (arcseconds) and \(\zeta\) (degrees) are constants (Table 2.3.8), and \(F_{i}\), \(f_{i} ^{{'} }\), and \(\delta _{i} ^{{'} }\) are amplitude, mean rate, and phase terms in the sine series expansion (Berger (1978a,b)) (series expansion terms not shown but can be found in the source code file shr_orb_mod.F90). The longitude of the perihelion \(\Pi\) (radians) depends on the sine and cosine series expansions for the eccentricity \(e\) as follows:

(2.3.89)\[\begin{split}\Pi =\left\{\begin{array}{lr} 0 & \qquad {\rm for\; -1}\times {\rm 10}^{{\rm -8}} \le e^{\cos } \le 1\times 10^{-8} {\rm \; and\; }e^{\sin } = 0 \\ 1.5\pi & \qquad {\rm for\; -1}\times {\rm 10}^{{\rm -8}} \le e^{\cos } \le 1\times 10^{-8} {\rm \; and\; }e^{\sin } < 0 \\ 0.5\pi & \qquad {\rm for\; -1}\times {\rm 10}^{{\rm -8}} \le e^{\cos } \le 1\times 10^{-8} {\rm \; and\; }e^{\sin } > 0 \\ \tan ^{-1} \left[\frac{e^{\sin } }{e^{\cos } } \right]+\pi & \qquad {\rm for\; }e^{\cos } <{\rm -1}\times {\rm 10}^{{\rm -8}} \\ \tan ^{-1} \left[\frac{e^{\sin } }{e^{\cos } } \right]+2\pi & \qquad {\rm for\; }e^{\cos } >{\rm 1}\times {\rm 10}^{{\rm -8}} {\rm \; and\; }e^{\sin } <0 \\ \tan ^{-1} \left[\frac{e^{\sin } }{e^{\cos } } \right] & \qquad {\rm for\; }e^{\cos } >{\rm 1}\times {\rm 10}^{{\rm -8}} {\rm \; and\; }e^{\sin } \ge 0 \end{array}\right\}.\end{split}\]

The numerical solution for the longitude of the perihelion \(\tilde{\omega }\) is constrained to be between 0 and 360 degrees (measured from the autumn equinox). A constant 180 degrees is then added to \(\tilde{\omega }\) because the Sun is considered as revolving around the Earth (geocentric coordinate system) (Berger et al. 1993).

Table 2.3.8 Orbital parameters

Parameter
\(\varepsilon^{*}\)	23.320556
\(\tilde{\psi }\) (arcseconds)	50.439273
\(\zeta\) (degrees)	3.392506