2.3. Surface Albedos¶
2.3.1. Canopy Radiative Transfer¶
Radiative transfer within vegetative canopies is calculated from the twostream approximation of Dickinson (1983) and Sellers (1985) as described by Bonan (1996)
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 nearinfrared (\(\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
where \(\phi _{1} ={\rm 0.5}0.633\chi _{L} 0.33\chi _{L}^{2}\) and \(\phi _{2} =0.877\left(12\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
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 snowcovered fraction \(f_{can,\, sno}\) (Chapter 2.7). The optical parameters are
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 twostream equations. If there is no snow on the canopy, this reduces to
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}\) )
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}\))
The upscatter for diffuse radiation is
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
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
where the single scattering albedo is
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
The downward diffuse fluxes per unit incident direct beam and diffuse radiation, respectively, are
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
These fluxes are partitioned to the sunlit and shaded canopy using an analytical solution to the twostream 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
and for shaded leaves is
with
For diffuse radiation, the absorbed radiation for sunlit leaves is
and for shaded leaves is
with
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)]:
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).
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 
C_{3} arctic grass 
0.30 
0.11 
0.35 
0.31 
0.53 
0.05 
0.34 
0.120 
0.250 
C_{3} grass 
0.30 
0.11 
0.35 
0.31 
0.53 
0.05 
0.34 
0.120 
0.250 
C_{4} grass 
0.30 
0.11 
0.35 
0.31 
0.53 
0.05 
0.34 
0.120 
0.250 
C_{3} 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 
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
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)
Unfrozen lake albedos depend on the cosine of the solar zenith angle \(\mu\)
Frozen lake albedos are from NCAR LSM (Bonan 1996)
As in NCAR LSM (Bonan 1996), soil albedos vary with color class
where \(\Delta\) depends on the volumetric water content of the first soil layer \(\theta _{1}\) (section 2.7.3) as \(\Delta =0.110.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 allsky surface albedo as described in Strahler et al. (1999) and Schaaf et al. (2002). The CLM twostream 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 allsky albedo in the twostream radiation model was selected as the best fit for the month. The fitted monthly soil colors were averaged over all snowfree months to specify a representative soil color for the grid cell. In cases where there was no snowfree surface albedo for the year, the soil color derived from snowaffected albedo was used to give a representative soil color that included the effects of the minimum permanent snow cover.
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 twostream 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 atmosphericdeposited 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 twostream solution requires the following bulk optical properties for each snow layer and spectral band: extinction optical depth (\(\tau\)), singlescatter 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 freshlyfallen 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
For each constituent (ice, two black carbon species, two organic carbon species, and four dust species), \(\omega\), g, and the mass extinction crosssection \(\psi\) (m^{2} 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 crosssection and layer mass, \(w _{k}\) (kgm^{1}) as
The twostream solution (Toon et al. (1989)) applies a tridiagonal 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 nearinfrared (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.75.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 midlatitude winter (Flanner et al. (2007)). The tridiagonal solution includes intermediate terms that allow for easy interchange of twostream 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 twostream solution
Spectral band 
Directbeam weight 
Diffuse weight 

Band 1: 0.30.7\(\mu\)m (visible) 
(1.0) 
(1.0) 
Band 2: 0.71.0\(\mu\)m (nearIR) 
0.494 
0.586 
Band 3: 1.01.2\(\mu\)m (nearIR) 
0.181 
0.202 
Band 4: 1.21.5\(\mu\)m (nearIR) 
0.121 
0.109 
Band 5: 1.55.0\(\mu\)m (nearIR) 
0.204 
0.103 
Under directbeam 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 directbeam 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
This weighting is performed for directbeam and diffuse, visible and NIR fluxes. After the groundincident 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 groundincident 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 namelistdefined 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 massextinction cross section (\(\bar{\psi }\)) over wavelength interval \(\lambda _{1}\) to \(\lambda _{2}\) is
Broadband singlescatter albedo (\(\bar{\omega }\)) is additionally weighted by the diffuse albedo for a semiinfinite snowpack (\(\alpha _{sno}\))
Inclusion of this additional albedo weight was found to improve accuracy of the fiveband albedo solutions (relative to 470band solutions) because of the strong dependence of opticallythick snowpack albedo on ice grain singlescatter 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. Singlescatter albedos for the endmembers of this size range are listed in Table 2.3.5.
Optical properties for black carbon are described in Flanner et al. (2007). Singlescatter albedo, mass extinction crosssection, 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).
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 
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 
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 areatovolume 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 areaweighted 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 waterinduced metamorphism (\(dr_{e,wet}\)), refreezing of liquid water, and addition of freshlyfallen snow. The mass of each snow layer is partitioned into fractions of snow carrying over from the previous time step (\(f_{old}\)), freshlyfallen snow (\(f_{new}\)), and refrozen liquid water (\(f_{rfz}\)), such that snow \(r_{e}\) is updated each time step t as
Here, the effective radius of freshlyfallen snow (\(r_{e,0}\)) is based on a simple linear temperaturerelationship. Below 30 degrees Celsius, a minimum value is enforced of 54.5 \(\mu\) m (corresponding to a specific surface area of 60 m^{2} 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 interparticle 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
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 midlayer temperatures (\(T_{n}\)) and snow layer thicknesses (\(dz_{n}\)), as
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_{n1}\) = \(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
The constant C_{1} 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 massweighted 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 301500\(\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
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
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
where \(\varepsilon\) is Earth’s obliquity and \(\lambda\) is the true longitude of the Earth.
The obliquity of the Earth \(\varepsilon\) (degrees) is
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)
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
where \(d_{ve} =80.5\) is the calendar day at vernal equinox (March 21 at noon), and
where \(\beta =\sqrt{1e^{2} }\). Earth’s eccentricity \(e\) is
where
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
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
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:
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)).
Parameter 


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