In spite of extensive experimental and theoretical studies of land surface albedo, there is a need of generalization and systematic quantification of its dependence on solar zenith angle and type of surface, and of atmospheric effects in these relationships. This work presents a systematic numerical study of surface albedo for 12 land cover types, described by realistic anisotropic reflectance, under widely varying atmospheric conditions. Our research is based on the rigorous solution of the radiative transfer equation with the spherical harmonics method. It uses a new original analytical formalism for treatment of the anisotropic boundary conditions, which avoids numeric integration. This formalism is described in the first part of the paper along with the major steps in solution of the radiative transfer equation. The second part of the paper presents the results of our study, which can be summarized as follows: for a large number of vegetation and soil surfaces, the range of relative variation of surface albedo with atmospheric opacity does not exceed 10–15% at θ0 ≤ 50° and 20–30% at θ0 ≤ 70°. Among atmospheric parameters, optical thickness is the main factor affecting albedo. The effects of the scattering function and single‐scattering albedo are small and can be neglected. For small solar zenith angles, surface albedo increases with atmospheric opacity, while for big angles θ0, it decreases. The watershed lies at θ0 ≈ 52°–57° when the albedo is almost insensitive to the atmospheric opacity. Albedo measurements at these angles can be used directly for the surface characterization. At large solar zenith angles (θ0 > 70°–80°), albedo drops off with the further increase of angle θ0.