Modeling
Light Pollution From Population Data and Implications for National Park Service
Lands
Steve Albers*
NOAA-FSL
325 Broadway
Boulder, CO 80305-3328
Dan Duriscoe
NPS Night Sky Team
Forest Ecologist
Sequoia and Kings Canyon National Parks
Three Rivers, CA 93271
dan_duriscoe@nps.gov
INTRODUCTION
There are many factors that affect night-time sky brightness, both natural and man-made. It is useful to think of what the main light sources are and how this light is scattered. The natural sources come from stars, the Milky Way, airglow, and sometimes moonlight. Manmade sources include streetlights and other outdoor lights, concentrated largely in towns and cities. Light is scattered by air molecules, natural and anthropogenic particulates, and haze (an enlargement of these particulates related to atmospheric moisture). The result of all these factors is what we see at night in terms of the sky brightness. To help clarify the further discussion, some simplifications will be helpful. We will assume no moonlight and relatively low levels of particulates and haze, namely that we are looking among the best nights from a given location. We also neglect things like surface albedo, that affects how much light is directed upward from city lights. The main remaining factor is city lights - related approximately to population, and natural airglow (a continuous auroralike glow) that actually varies during the course of the sunspot cycle. The darkest sites have a brighter glow than being in outer space for two main reasons, the scattering of starlight by the atmosphere, and airglow.
Brief Review of previous modeling efforts
A number of people have modeled light pollution in various ways. As an example, Roy Garstang (1986) has done detailed calculations for a number of observatory sites, creating maps showing how the skyglow varies at different altitudes and azimuths from each site.
William Burton (2000) is analyzing DMSP (Defense Meteorological Satellite Program, run by the U.S. Air Force) satellite data to estimate sky glow in the close vicinity of urban areas. This has the advantage of considering actual satellite data at high resolution both spatially and in terms of intensity. However, limited consideration is given to atmospheric scattering, especially over large distances.
DMSP satellite data has been linked with a robust scattering model in Europe by Cinzano et. al. (2000,2001). When properly calibrated, this provides greater spatial information about the sources of light pollution than population data alone.
Description of the model
Assumptions
and data source
The present effort is unique in that it produces an areal map of zenithal sky brightness over the entire United States. It works both within and at large distances from urban centers. This also employs assumptions about scattering embodied in Walker's Law with additional consideration of the earth's curvature. The model is based on the location and population (1990 census) of significant U.S. cities and towns (over 50 population).
Mechanics
of the model
The model creates a map of expected sky glow at the zenith. For each location in the map the light pollution contribution from each city is assumed to be related linearly to the population and the inverse 2.5 power of the distance. This is similar to the relation used in Walker's law (Walker 1977), except that we are estimating light pollution at the zenith instead of 45 degrees high in the azimuth of the brightest city. The relation used here for each city i is in equation 1, where
Ii = 11,300,000pr-2.5 (1)
Ii is sky glow in nanoLamberts, p is the city population, and r is the distance to the city in meters. This is corrected for earth's curvature at large distances (this necessity was pointed out by Roy Garstang, 1989). The correction is done by calculating the fraction f of the air molecules and other scatterers over the observer that lie above the earth's shadow that is formed from light traveling in a straight line from the city. The overall scale height s for these scatterers (defined as the altitude increase required to see a drop off by a factor of e) is currently set to 4000m. This is less than the "clear air" value of 8000m accounting for a typical amount of aerosols. The scattering from this mixture is more strongly concentrated at low altitudes than that from air molecules alone.
f = e-h/s (2)
The height h is the amount of the air column above the observing location that is not in a direct line of sight to the light polluting city due to the curvature of the earth. Adding this correction term f into equation 1 yields a further modified form of Walker's law as shown in equation 3.
Ii = 11300000pr-2.5f (3)
Finally the light pollution from each city is summed to get the total artifical skyglow I at a given location on the light pollution map as in equation 4.
I = summation Ii (4)
A number of other ideas and equations (such as conversions between sky brightness and limiting magnitude) were also adapted from publications by Roy Garstang (1986,1989).
The assumed radius of each city is a function of city population, ranging from 2.5 km to 24 km. Walker's law applies if we are outside the city radius. Inside the city radius, the sky glow increases linearly toward the center by another factor of 2.5.
A value for the natural sky glow is added onto the light pollution contribution. The natural sky glow is assumed to be equal to 60 nanoLamberts (V = 21.9 mag / sq sec) at solar minimum. The last step in arriving at a pixel value is scaling the brightness with respect to the logarithm of the sky glow. The brightest city has pixel values of (255,255,255), and the darkest country site has pixel values of (42,52,67). The calibration bar is intended to linearly represent the sky brightness in terms of magnitudes per square arcsec. The model result is shown in Figures 1a and 1b.

Fig 1a. Modeled skyglow over the continental United States.
Schaaf Sky
Quality Scale
Fred Schaaf has provided much helpful discussion that has led to substantial improvements in this image. As part of the calibration process of the image, we are comparing the expected amount of light pollution for various locations with observations of limiting magnitude and sky quality according to the Schaaf Scale (Schaaf 1994).
Table 1. The Schaaf Scale. Within the context of this light pollution model, the following conversion can be done between Schaaf Scale and Zenithal Limiting Magnitude (ZLM). This has been modified slightly from the original scale so that the verbal sky descriptions in the Schaaf Scale are more consistent with modeled limiting magnitudes.
|
Schaaf Class |
Zenithal Limiting Magnitude |
|
1 |
<4.75 |
|
2 |
4.75-5.25 |
|
3 |
5.25-5.75 |
|
4 |
5.75-6.20 |
|
5 |
6.20-6.55 |
|
6 |
6.55-6.76 |
|
7 |
6.76-6.81 |

Fig 1b. Same as Figure 1a except that it is rescaled to represent the seven level
Schaaf Scale.
Verification of the model with observations
Field data has been received from locations throughout the United States, from the International Meteor Organization and other interested observers. Most of this data is in the form of zenith limiting magnitude, and can be compared to values predicted by the model for the same location. We also are doing comparisons of some sky brightness observations, mainly from Garstang (1986, 1989). Figure 2 is a scatterplot of observations received near sunspot maximum compared with predicted limiting magnitude for the site. The graph shows that the predicted values are in relatively good agreement with observed, especially between magnitudes 5 and 6. For very dark skies, observers typically do not see stars as faint as the model predicts, perhaps relating to differing visual acuity and weather conditions among observers. For brighter skies, observers consistently see stars fainter than the model predicts, implying the model may need some adjustment within urban areas.

Figure 2. Plot of predicted vs. observed limiting magnitude for observations near sunspot maximum. The solid line represents observed=predicted, the dashed fine lines show one magnitude lower and higher than predicted, and the bold dashed line is the best fit linear regression of the data.
Future Work
A natural extension of this work would be to incorporate DMSP satellite data, perhaps with the population data, to gain a better idea of where the light sources are in the United States. This has been done over Europe and may be extended over other parts of the world by Cinzano et. al. (2001). The use of population data from other census decades could provide a time series of light pollution in the United States over many years.
NATIONAL PARK Land Area Analysis
The model may be used to evaluate the effects of light pollution on areas administered by the National Park Service for the purpose of protecting night sky visibility. When the image generated by the model is imported into a Geographic Information System such as ArcInfo, and the park boundaries superimposed, a simple intersection of the two themes yields data on the relative proportion of each park that falls within each of the Schaaf scale classes. Figures 3 through 6 show selected regions of the United States (in a negative image for clarity) with National Park Service areas superimposed upon the light pollution model. The state boundaries are also added, and the maps are produced in Lambert's Conformal Conic projection.

Figure 3. Light pollution and national parks in the southwest.

Figure 4. Light pollution and national parks in the northwest.

Figure 5. Light pollution and national parks in the northeast.

Figure 6. Light pollution and national parks in the southeast.

Figure 7. Legend for figures 3-6.
Visual examination of these maps reveals that large areas of the west were predicted to still possess Schaaf class 7 skies, or "pristine", while such sites occurred very rarely east of the Mississippi River as of 1990. Large areas of class 7 skies are seen in Nevada, Montana, North Dakota, eastern Oregon, southeastern Utah/northern Arizona, west Texas, and Wyoming. These regions were far enough from large urban centers that the influence of light pollution was minimal. Several large and well known national parks fall within these areas, including Glacier, Yellowstone, Canyonlands, Grand Canyon, and Death Valley.
The GIS analysis produced a table of park areas showing what percentage of the area within each park fell within each of the Schaaf classes. Also, a mean Schaaf class was computed, and the total acreage of each park was calculated (Table 2). Not all park areas are shown, many were edited out for brevity; the authors apologize in advance if your favorite park was left out. The table is ordered first from darkest to brightest mean Schaaf class, then alphabetically by park name for parks with identical means. Note that the acreage may be off by as much as 2-3% of the true acreage because of the scale of the park boundary data used for the analysis.
Table 2. Analysis of sky quality at selected National Park Service areas.
|
Park Name |
Total Acres |
Mean Schaaf Class |
Percentage of Land Area Within Each Schaaf Class |
||||||
|
1 |
2 |
3 |
4 |
5 |
6 |
7 |
|||
|
Badlands NP |
241,284 |
7.00 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
100.0 |
|
Big Bend NP |
827,169 |
7.00 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
100.0 |
|
Canyonlands NP |
331,342 |
7.00 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
100.0 |
|
Capitol Reef NP |
241,505 |
7.00 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
100.0 |
|
Carlsbad Caverns NP |
46,921 |
7.00 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
100.0 |
|
Chaco Culture NHP |
34,504 |
7.00 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
0.5 |
99.5 |
|
Chiricahua NMON |
12,225 |
7.00 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
100.0 |
|
Crater Lake NP |
180,631 |
7.00 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
100.0 |
|
Craters of the Moon NMON |
750,312 |
7.00 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
100.0 |
|
Devils Tower NMON |
1,341 |
7.00 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
100.0 |
|
Death Valley NP |
3,370,969 |
7.00 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
100.0 |
|
Dinosaur NMON |
208,650 |
7.00 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
0.1 |
99.9 |
|
Dry Tortugas NP |
72,382 |
7.00 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
100.0 |
|
Gila Cliff Dwellings NMON |
526 |
7.00 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
100.0 |
|
Glacier NP |
1,026,615 |
7.00 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
100.0 |
|
Great Basin NP |
76,349 |
7.00 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
100.0 |
|
Great Sand Dunes NP |
38,202 |
7.00 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
100.0 |
|
Guadalupe Mountains NP |
88,254 |
7.00 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
100.0 |
|
Hovenweep NMON |
797 |
7.00 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
100.0 |
|
Isle Royale NP |
143,269 |
7.00 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
100.0 |
|
Lava Beds NMON |
46,004 |
7.00 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
100.0 |
|
Lassen Volcanic NP |
106,239 |
7.00 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
100.0 |
|
Natural Bridges NMON |
7,324 |
7.00 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
100.0 |
|
Navajo NMON |
597 |
7.00 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
100.0 |
|
North Cascades NP |
510,531 |
7.00 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
100.0 |
|
Organ Pipe Cactus NMON |
331,119 |
7.00 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
100.0 |
|
Petrified Forest NP |
94,430 |
7.00 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
100.0 |
|
Rainbow Bridge NMON |
161 |
7.00 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
100.0 |
|
Theodore Roosevelt NP |
71,048 |
7.00 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
100.0 |
|
Voyageurs NP |
208,263 |
7.00 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
100.0 |
|
Wupatki NMON |
36,164 |
7.00 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
100.0 |
|
Grand Teton NP |
308,640 |
6.99 |
0.0 |
0.0 |
0.0 |
0.0 |
0.0 |
0.5 |
99.5 |
|
Yellowstone NP |
2,197,269 |
6.99 |
0.0 |
0.0 |
0.0 |
0.0 |
0.2 |
0.3 |
99.5 |
|
Grand Canyon NP |
1,197,475 |
6.98 |
0.0 |
0.0 |
0.0 |
0.0 |
0.5 |
1.2 |
98.3 |