An analytical model characterizing the atmospheric radiance field over a non-Lambertian surface divides the radiation field into three components: unscattered radiance, single-scattering radiance, and multiple-scattering radiance. The first two components are calculated exactly. A δ-four-stream approximation is extended to calculate the azimuth-independent multiple-scattering radiance over a non-Lambertian surface, which is modeled by a statistical bidirectional reflectance distribution function (BRDF). Accuracy is assessed with respect to the exact results computed from a Gauss-Seidel iterative algorithm. Experiments comparing the results obtained with Lambertian and non-Lambertian surfaces show that incorporating the BRDF into the four-stream approximation significantly improves the accuracy in calculating radiance as well as radiative flux.