2.3. Surface Albedos

2.3.1. Canopy Radiative Transfer

Radiative transfer within vegetative 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 } =\alpha _{\Lambda }^{leaf} w_{leaf} +\alpha _{\Lambda }^{stem} w_{stem}\]

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

(2.3.12)\[\tau _{\Lambda } =\tau _{\Lambda }^{leaf} w_{leaf} +\tau _{\Lambda }^{stem} w_{stem} .\]

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.

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 } e^{-K\left(L+S\right)} +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)e^{-K\left(L+S\right)}\]
(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 trees and shrubs are from Dorman and Sellers (1989). Leaf and stem optical properties (VIS and NIR reflectance and transmittance) were derived for grasslands and crops 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.07

0.35

0.16

0.39

0.05

0.10

0.001

0.001

NET Boreal

0.01

0.07

0.35

0.16

0.39

0.05

0.10

0.001

0.001

NDT Boreal

0.01

0.07

0.35

0.16

0.39

0.05

0.10

0.001

0.001

BET Tropical

0.10

0.10

0.45

0.16

0.39

0.05

0.25

0.001

0.001

BET temperate

0.10

0.10

0.45

0.16

0.39

0.05

0.25

0.001

0.001

BDT tropical

0.01

0.10

0.45

0.16

0.39

0.05

0.25

0.001

0.001

BDT temperate

0.25

0.10

0.45

0.16

0.39

0.05

0.25

0.001

0.001

BDT boreal

0.25

0.10

0.45

0.16

0.39

0.05

0.25

0.001

0.001

BES temperate

0.01

0.07

0.35

0.16

0.39

0.05

0.10

0.001

0.001

BDS temperate

0.25

0.10

0.45

0.16

0.39

0.05

0.25

0.001

0.001

BDS boreal

0.25

0.10

0.45

0.16

0.39

0.05

0.25

0.001

0.001

C3 arctic grass

-0.30

0.11

0.35

0.31

0.53

0.05

0.34

0.120

0.250

C3 grass

-0.30

0.11

0.35

0.31

0.53

0.05

0.34

0.120

0.250

C4 grass

-0.30

0.11

0.35

0.31

0.53

0.05

0.34

0.120

0.250

C3 Crop

-0.30

0.11

0.35

0.31

0.53

0.05

0.34

0.120

0.250

Temp Corn

-0.50

0.11

0.35

0.31

0.53

0.05

0.34

0.120

0.250

Spring Wheat

-0.50

0.11

0.35

0.31

0.53

0.05

0.34

0.120

0.250

Temp Soybean

-0.50

0.11

0.35

0.31

0.53

0.05

0.34

0.120

0.250

Cotton

-0.50

0.11

0.35

0.31

0.53

0.05

0.34

0.120

0.250

Rice

-0.50

0.11

0.35

0.31

0.53

0.05

0.34

0.120

0.250

Sugarcane

-0.50

0.11

0.35

0.31

0.53

0.05

0.34

0.120

0.250

Tropical Corn

-0.50

0.11

0.35

0.31

0.53

0.05

0.34

0.120

0.250

Tropical Soybean

-0.50

0.11

0.35

0.31

0.53

0.05

0.34

0.120

0.250

Miscanthus

-0.50

0.11

0.35

0.31

0.53

0.05

0.34

0.120

0.250

Switchgrass

-0.50

0.11

0.35

0.31

0.53

0.05

0.34

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 are from Paterson (1994)

\[\alpha _{soi,\, vis}^{\mu } =\alpha _{soi,\, vis} =0.6\]
\[\alpha _{soi,\, nir}^{\mu } =\alpha _{soi,\, nir} =0.4.\]

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

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

Frozen lake albedos are from NCAR LSM (Bonan 1996)

\[\alpha _{soi,\, vis}^{\mu } =\alpha _{soi,\, vis} =0.60\]
\[\alpha _{soi,\, nir}^{\mu } =\alpha _{soi,\, nir} =0.40.\]

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

(2.3.61)\[\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).

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), which incorporates a two-stream radiative transfer solution from Toon et al. (1989). Albedo and the vertical 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), and ice effective grain size (\(r_{e}\)), which 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 in CLM is also described somewhat by Flanner and Zender (2005) and Flanner et al. (2007).

The 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). The bulk optical properties are weighted functions of each constituent k, computed for each snow layer and spectral band as

(2.3.62)\[\tau =\sum _{1}^{k}\tau _{k}\]
(2.3.63)\[\omega =\frac{\sum _{1}^{k}\omega _{k} \tau _{k} }{\sum _{1}^{k}\tau _{k} }\]
(2.3.64)\[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\) (m2 kg-1) are computed offline with Mie Theory, e.g., 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.65)\[\tau _{k} =\psi _{k} w_{k}\]

The 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. Solar fluxes are computed in five spectral bands, listed in Table 2.3.4. Because snow albedo varies strongly across the solar spectrum, it was determined that four bands were needed to accurately represent the near-infrared (NIR) characteristics of snow, whereas only one band was needed for the visible spectrum. Boundaries of the NIR bands were selected to capture broad radiative features and maximize accuracy and computational efficiency. We partition NIR (0.7-5.0 \(\mu\) m) surface downwelling flux from CLM according to the weights listed in Table 2.3.4, which are unique for diffuse and direct incident flux. These fixed weights were determined with offline hyperspectral radiative transfer calculations for an atmosphere typical of mid-latitude winter (Flanner et al. (2007)). The tri-diagonal solution includes intermediate terms that allow for easy interchange of two-stream techniques. We apply the Eddington solution for the visible band (following Wiscombe and Warren 1980) and the hemispheric mean solution ((Toon et al. (1989)) for NIR bands. These choices were made because the Eddington scheme works well for highly scattering media, but can produce negative albedo for absorptive NIR bands with diffuse incident flux. Delta scalings are applied to \(\tau\), \(\omega\), and \(g\) (Wiscombe and Warren 1980) in all spectral bands, producing effective values (denoted with \(*\)) that are applied in the two-stream solution

(2.3.66)\[\tau ^{*} =\left(1-\omega g^{2} \right)\tau\]
(2.3.67)\[\omega ^{*} =\frac{\left(1-g^{2} \right)\omega }{1-g^{2} \omega }\]
(2.3.68)\[g^{*} =\frac{g}{1+g}\]
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.586

Band 3: 1.0-1.2\(\mu\)m (near-IR)

0.181

0.202

Band 4: 1.2-1.5\(\mu\)m (near-IR)

0.121

0.109

Band 5: 1.5-5.0\(\mu\)m (near-IR)

0.204

0.103

Under direct-beam conditions, singularities in the radiative approximation are occasionally approached in spectral bands 4 and 5 that produce unrealistic conditions (negative energy absorption in a layer, negative albedo, or total absorbed flux greater than incident flux). When any of these three conditions occur, the Eddington approximation is attempted instead, and if both approximations fail, the cosine of the solar zenith angle is adjusted by 0.02 (conserving incident flux) and a warning message is produced. This situation occurs in only about 1 in 10 6 computations of snow albedo. After looping over the five spectral bands, absorption fluxes and albedo are averaged back into the bulk NIR band used by the rest of CLM.

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 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.69)\[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 spectral bands are derived offline and stored in a namelist-defined lookup table for online retrieval (see CLM5.0 User’s Guide). Mie properties are first computed at fine spectral resolution (470 bands), and are then weighted into the five bands applied by CLM according to incident solar flux, \(I^{\downarrow } (\lambda )\). For example, the broadband mass-extinction cross section (\(\bar{\psi }\)) over wavelength interval \(\lambda _{1}\) to \(\lambda _{2}\) is

(2.3.70)\[\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.71)\[\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 470-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 are described in Flanner et al. (2007). Single-scatter albedo, mass extinction cross-section, and asymmetry parameter values for all snowpack species, 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. Weighting into the five CLM spectral bands was determined only with incident solar flux, as in equation (2.3.69).

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.516

0.434

0.346

0.276

0.139

Hydrophobic black carbon

0.288

0.187

0.123

0.089

0.040

Hydrophilic organic carbon

0.997

0.994

0.990

0.987

0.951

Hydrophobic organic carbon

0.963

0.921

0.860

0.814

0.744

Dust 1

0.979

0.994

0.993

0.993

0.953

Dust 2

0.944

0.984

0.989

0.992

0.983

Dust 3

0.904

0.965

0.969

0.973

0.978

Dust 4

0.850

0.940

0.948

0.953

0.955

Ice (\(r _{e}\) = 30 \(\mu\) m)

0.9999

0.9999

0.9992

0.9938

0.9413

Ice (\(r _{e}\) = 1500 \(\mu\) m)

0.9998

0.9960

0.9680

0.8730

0.5500
Table 2.3.6 Mass extinction values (m2 kg-1) used for snowpack impurities and ice

Species

Band 1

Band 2

Band 3

Band 4

Band 5

Hydrophilic black carbon

25369

12520

7739

5744

3527

Hydrophobic black carbon

11398

5923

4040

3262

2224

Hydrophilic organic carbon

37774

22112

14719

10940

5441

Hydrophobic organic carbon

3289

1486

872

606

248

Dust 1

2687

2420

1628

1138

466

Dust 2

841

987

1184

1267

993

Dust 3

388

419

400

397

503

Dust 4

197

203

208

205

229

Ice (\(r _{e}\) = 30 \(\mu\) m)

55.7

56.1

56.3

56.6

57.3

Ice (\(r _{e}\) = 1500 \(\mu\) m)

1.09

1.09

1.09

1.09

1.1
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.52

0.34

0.24

0.19

0.10

Hydrophobic black carbon

0.35

0.21

0.15

0.11

0.06

Hydrophilic organic carbon

0.77

0.75

0.72

0.70

0.64

Hydrophobic organic carbon

0.62

0.57

0.54

0.51

0.44

Dust 1

0.69

0.72

0.67

0.61

0.44

Dust 2

0.70

0.65

0.70

0.72

0.70

Dust 3

0.79

0.75

0.68

0.63

0.67

Dust 4

0.83

0.79

0.77

0.76

0.73

Ice (\(r _{e}\) = 30\(\mu\)m)

0.88

0.88

0.88

0.88

0.90

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.72)\[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 m2 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.73)\[\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.74)\[\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.75)\[\frac{dr_{e} }{dt} =\frac{10^{18} C_{1} f_{liq} ^{3} }{4\pi r_{e} ^{2} }\]

The constant C1 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.76)\[\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.77)\[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.78)\[\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.79)\[\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.80)\[\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.81)\[\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.82)\[\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.83)\[e=\sqrt{\left(e^{\cos } \right)^{2} +\left(e^{\sin } \right)^{2} }\]

where

(2.3.84)\[\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)). The longitude of the perihelion relative to the moving vernal equinox \(\tilde{\omega }\) (degrees) is

(2.3.85)\[\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.86)\[\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))). The longitude of the perihelion \(\Pi\) (radians) depends on the sine and cosine series expansions for the eccentricity \(e\)as follows:

(2.3.87)\[\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