Jay Cammermeyer - Paul Conrecode - Jeff Hansen - Pierre Kwan - Miranda Maupin

Nutrient loading to Lakes Washington and Sammamish is primarily a function of the surrounding non-point sources. Diversion of secondary effluent form both basins n the 1960’s resulted in rapid recovery of water quality (Edmonson, 1994). Lower chlorophyll a counts, reduced levels of total phosphorus, and greater Secchi disk transparency following recovery, demonstrated that external phosphorus loading largely determines lake water quality.

Recently, accelerated rates of urban development in the Lake Washington/Lake Sammamish watersheds are suspected of altering the total phosphorus load to the lakes. Changes in land use have increased the impervious surface area and altered stormwater runoff dynamics, which are probable causal agents of current lake quality decline. The spatially-based phosphorous loading model provides a tool to understand the relationship between land use change and phosphorous loading.

In addition to being a component of PRISM (Puget Sound Regional Synthesis Model), the model will support SWAMP (Sammamish Washington Analysis and Modeling Program). SWAMP has multiple objectives related to non-point pollution including evaluation of indirect water reuse, risks for drinking water supply, and risks for fisheries.

This report outlines the work completed and future steps to accomplish the following objectives.

1. Research

• Identify a method to model phosphorus transfer from land to lakes through a review of existing spatial models
• Determine land use specific phosphorus flux coefficients

1. Model Development

• Compile existing spatial data
• Create new spatial data layers
• Compute algorithms

1. Model Calibration

• Determine phosphorous concentration and loading data
• Adjust Transmission coefficients to fit data from reference watersheds

For large watersheds with heterogeneous land use, phosphorus loading models based on physical mechanisms driving phosphorus export have been largely unsuccessful in predicting stream phosphorus concentrations. Spatially distributed parameter models such as AGNPS and CREAMS are widely used to predict the effects on water quality due to nutrient and sediment loading. These models rely heavily on runoff equations that are based on cover type, antecedent runoff condition (ARC), and the US Soil Conservation Service’s hydrologic soil group classification. Total phosphorus loading or concentration is then calculated as a function of surface runoff. The spatial model (Saunders et al, 1996) applied to the San Antonio-Nueces Coastal Basin, Texas uses a runoff method to determine nutrient loadings. However, the relationship between runoff and total phosphorus export is highly dependent on the phosphorus source of the preferential adsorption of phosphorus to particular soil fractions, and the soil biota.

An alternative to a runoff-based model is an attenuation-based model. Attenuation based models account for the diminished export of phosphorus with increased distance from the surface water body. Land use determines what proportion of phosphorus will be transported overland a given distance; for impervious surfaces, the model used in this study assumes all phosphorus is transported to surface waters and (at the other extreme) a very low proportion of total phosphorus is transported across forested areas.

To account for the relative distance of phosphorus source from water body, several approaches have been suggested. Models from the 1970’s used a linear relationship between distance of source from shoreline and phosphorus export to the waterbody. Hakanson (1990) utilizes the drainage area zonation (DAZ) method to account for the decreased effective phosphorus loading with increasing distance from water body. Contour lines are drawn around the water body at increasing distances to delineate zones of reduced loading impacts. Separate contour maps can be developed for each land use type to account for differences in attenuation abilities. Soranno (1996) includes a transmission coefficient to characterize the attenuation effects of land use types. The transmission coefficient is also adjusted to reflect seasonal variations in precipitation.

The Lake Mendota attenuation-based model developed by Soranno serves as the template for this Lake Washington/Sammamish study. Further comparison between the Lake Mendota attenuation and San Antonio runoff model can be found in Table 6 in the Appendix.

It is important to distinguish between the phosphorus export coefficient, transmission coefficient, and phosphorus flux coefficient. The phosphorus export coefficient (c) is in units of kilograms per hectare per year, and is estimated by monitoring runoff from a particular land use and then dividing by the area drained (Dillon and Kirchner 1975, Rast and Lee 1978, Reckhow et al. 1980, Clesceri et al. 1986). The phosphorous export coefficient is part of the total phosphorous loading equation:

L is total phosphorus loading from land (in kilograms per year), m is number of land use types, c_i is the phosphorus export coefficient for land use i (in kilograms per hectare per year), and A_i is area of land use i (in hectares).

However, phosphorus loading as measured by the sum of the products of each land use's phosphorus export coefficient and each land use's area has considerable uncertainty. This uncertainty is due in part to the attenuation of phosphorus export from the watershed to the surface water bodies. That is, phosphorus export is not proportional to area for all land uses.

Soranno et al. (1996) account for attenuation over distance by including a transmission coefficient (T) that represents the proportion of phosphorus transported to the next pixel in the path of surface overland flow. The phosphorus loading equation can be modified:

T is the transmission coefficient (O<T<1) representing the proportion of phosphorus transported to the next pixel. The coefficient f is used (for phosphorus flux coefficient) instead of c (phosphorus export coefficient) because phosphorus flux is the phosphorus produced and transported to the next pixel. In contrast, phosphorus export coefficients assume 100% of phosphorus transported from land will reach surface water. The total number of pixels is represented by n, and m is the total number of land uses.

One step in the development of the spatial model entails assigning phosphorus flux coefficients to land use/land cover categories. These coefficients need to quantify the mass of total phosphorus fluxing off any hectare of a particular land use type and therefore are closely related to the traditional phosphorus export coefficients found in literature. The large number of particulate and dissolved forms of phosphorus that can be found in freshwater bodies and soil profiles makes attempts to predict loading rates to the terminal lake of a large (>100,000 ha) watershed quite difficult. Haygarth (1997) discusses the numerous phosphorus transport pathways that occur across many temporal scales, the nearly limitless adsorption capacities of many soils, and the variable effects of storms due to duration, frequency, and intensity. It has only been recently that researchers have attempted to identify the essential processes that drive the majority of nutrient transport from a landscape.

Previous studies have simply reported the nutrient concentration of outflow from small watershed studies, as measured in lysimeters or from stream grab samples. From these data, simple phosphorus export coefficients have been developed for the dominant land use type of the study watershed. The majority of studies have come from intensively farmed regions of the north-central or southeastern United States as well as from rural counties of the United Kingdom. The effects of crop types, tillage practices, mean slope, and fertilizer application rates on nitrogen and phosphorus export from small agricultural catchments have been well-documented (Reckhow & Beaulac, 1980). However there have been few studies that have examined phosphorus export coefficients in urban areas or the Pacific Northwest.

For the purposes of this study, the P-flux coefficients are based on studies conducted in small watersheds where the effects of attenuation are minimal. In these settings, it is assumed that phosphorus flux coefficients are very close to the values of phosphorus export coefficients.

Perkins (1997) points out that phosphorus export coefficient ratios for different land uses are slightly more consistent across studies than are the absolute phosphorus export coefficient values. Therefore, we attempted to determine the mean phosphorus export from one hectare of a typical western Cascade forest parcel. From this value, ratios can be used to determine phosphorus export coefficients for other land use classifications. Table 2 in the appendix summarizes phosphorus export coefficients from studies that are most relevant to the Lake Washington watershed. Based on these data and group discussions, we have selected the following values as phosphorus export coefficients for the Lake Washington/Lake Sammamish watershed:

Land Use Type	Phosphorus Export Coefficient (kg/ha/yr)	Land Use Type: Forest Ratio
Forest	0.30	1.00
Agriculture	1.50	5.00
Grazing	1.30	4.33
Urban (avg.): Commercial HD Residential LD Residential	1.80 2.25 1.45 1.25	6.00 7.50 4.80 4.20

Table 1. Phosphorus Export Coefficients Selected for Lake Washington/Lake Sammamish Watershed

Compiling spatial data takes considerable time, expertise and sophisticated software and hardware systems. Many individuals have been of great assistance to this project in this respect. From within the University of Washington, these people include Marina Alberti of the Department of Urban Planning, Kristina Hill from the Department of Landscape Architecture, Miles Logsdon of the School of Oceanography, and Ulysses Hillard and James VanShaar of the Department of Civil and Environmental Engineering. Also very generous in lending time and advice was Greg Pelletier of Washington’s Department of Ecology.

The team encountered difficulties as model development progressed due to hardware and software limitations. Initial work was performed using ArcView on a WinNT system. However, the more powerful computational resources of the Unix based ArcInfo software and Solaris workstations were needed to reproject the coverages so that data from multiple sources could be integrated. The remainder of the project tasks, which included the actual modeling, took place within the Windows-based ArcView program.

The applicable GIS data was obtained from several sources. Miles Logsdon supplied the USGS digital elevation model (DEM) encompassing the entire Lake Washington-Sammamish watershed. This data, in the form of ArcInfo coverages, is the base layer for other PRISM work being done at the University of Washington. Logsdon also supplied land cover data for the 1991 Lake Washington/Lake Sammamish watershed. This landcover from the 1991 U.S. Fish & Wildlife Service GAP analysis was classified for the purposes of habitat preservation. The classifications are less effective for hydrologic modeling because the classifications in the urban areas are too coarse to model land use changes (Figure 1). For example, the image could easily delineate between areas of medium density coniferous shrubs and areas of medium density coniferous trees but it agglomerated the entire area from Kirkland to Renton and Mercer Island as medium density mixed use. Instead, the team decided use landcover data produced by Dr. Kristina Hill specifically for hydrologic modeling. This data, which distinguishes urban areas more finely, included 1991 and 1998 classified Landsat images for the study area (Figure 2). Additional data and consultation were supplied by Greg Pelletier.

The basic GIS-based steps in the development of the model are outlined below, followed by a discussion of each one. This process is identified within the Soranno article.

1. Create distance layer, which assigns a value to each cell based on the elevation distance in pixels to the nearest water body.
2. Assign P-flux and transmission coefficients to existing land cover classifications.
3. Apply model summation equation to distance and coefficient layers.

The overall goal of developing the "distance" layer is to determine for each pixel in the watershed the distance from that pixel along its downstream flowpath to the nearest stream or lake. First, a digital elevation model (DEM) of the watershed was created in 10 x 10m grid cells with elevations at 1-cm increments with ArcInfo. The DEM was artificially smoothed with an ArcInfo function to remove any sinks in the watershed and ensure all water would drain. While this procedure removes any natural depressions in the landscape where water would otherwise collect and percolate down into the ground water, the ArcInfo and ArcView programs can not account effectively for any water sink other than a lake or similar large water bodies. The exact same "smoothing" procedure is being performed in the DHSVM portion of PRISM and was done in the Sorrano study of Lake Mendota.

The next step involved "burning in" the stream networks onto the DEM. This commonly-used procedure was developed by the University of Texas-Austin’s Department of Civil Engineering to force the computer model to follow known stream paths. The DEM was overlaid with a Washington State Department of Natural Resource (DNR) stream network map and wherever DEM was not a stream, the elevation was raised 40,000 m. The end result looks like the normal landscape except all the streams would be very deep and narrow canyons in which water would flow.

The DEM was used to develop a flow direction model, which codes pixels according to the cardinal flow direction using interpolated hill slopes and directions. A flow accumulation layer was then developed based on the DEM and flow direction layers. This sums the number of pixels upland of an individual pixel. This is analogous to standing at a stream bank and determining the land area that is being drained through that point.

The flow direction layer was then used in conjunction with a layer that contained the stream network within the watershed to develop the distance layer. The stream network layer was a coverage that simply assigned all pixels that represented streams a value of 40,000.

A series of avenue scripts (or basic code) is used to perform these functions, a few of which are described here. First, the script outlined in Figure 4 merged the stream layer with the flow direction layer.

Next, the script outline in figure 5 defines all pixels that represent streams (which were previously given an artificial value of 40,000) as null points.

The final step, outlined in Figure 6, is to calculate the flow length from each pixel to the nearest surface water body.

At this point, the distance layer has been created. Next, appropriate P-flux coefficients based on literature reviews must be assigned to the land cover classifications. Table 2 illustrates the difference between the classifications used in the literature and the classifications of the land use/land cover spatial data. Phosphorus export coefficients in the literature are based on land use categories that do not directly relate to the land cover categories based on satellite detected light reflectance. The assigned values column in Table 2 displays the values the team chose to assign the land cover classifications for this study. Figure 7 outlines the script used to assign P-flux values to land classifications.

After all these necessary layers are complete, they can be multiplied based on the Soranno formula (L = å å f*A*T^p) to derive a phosphorus loading value for each pixel (A) based on pflux (f) and phosphorus transport (T) at distance (p). The resulting grid layer has a value assigned to each pixel which is equal to the amount of phosphorous contributed by that pixel (kg P/yr). The second script in Figure 7 outlines this process. Figure 8 shows an example of such an output grid layer. Then the loading for each sub-watershed is summed in a spreadsheet.

The model was run on several sub-watersheds with some interesting results. Most effort was placed upon investigating Issaquah Creek. An initial guess at modeling the P-load due to this sub-watershed resulted in a P-load very close to that calculated by the P-export team. P-flux and transmission coefficients were then modified to determine what effect upon the overall result these variables might have. The following table shows these results.

Note that the resultant P load is very sensitive to the classification scheme. The classification of Test Number 1 results in a value very close to the P-export group’s average calculation for Issaquah Creek from 1991 to 1995 (5,725 kg/yr). When the same P-flux and transmission coefficients are used for Juanita Creek, however, the model results do not match up well with the P-export average value of 670 kg/yr.

Table 4 summarizes the results of the model runs made for Issaquah, Juanita and Bear Creeks and the P-export group’s calculations.

To make more sense out of these results, a tabulation of the land coverage for these creeks for the two years of investigation was made:

Land Cover Type	Juanita Ck 1991	Juanita Ck 1998	Lower Bear Ck 1991	Lower Bear Ck 1998	Issaquah Ck 1991	Issaquah Ck 1998
Forested Urban	46	47	34	33	14	15
Grassy Urban	15	25	6	9	3	5
Paved Urban	19	11	5	3	5	3
Forest	12	14	46	50	74	74
Ag/Grass	8	3	9	5	4	3

There are many interesting issues to look at concerning the above data. First, it is quite obvious that the land cover classification scheme is suspect. For instance, the greatest percent change in land cover for the Juanita Creek sub-watershed involves an increase in "grassy urban", but also an associated decrease in "paved urban" of 8%. This does not make sense, seeing as that area of the watershed has seen tremendous growth in the past 10 years. Similarly, percent forested areas appear to increase with time in both Juanita and Lower Bear Creek sub-watersheds. Again, this is something that has not been observed in terms of land use in these regions over the past decade. There is something incorrect in the satellite image interpretation that is leading to this kind of conclusion.

Aside from land cover classification issues, it is apparent that the P-flux coefficients are not correct. It is only coincidental that the modeled output for Issaquah Creek matches up with the P-export group’s calculations. It must be remembered that no calibration was done to arrive at this outcome. The P-flux numbers are suspect because while the P load appears to match up well for Issaquah Creek, the model output for Juanita Creek grossly overestimates P load, as compared to the P-export group’s output.

All of these inconsistencies between the modeled and calculated P loads are most likely due to multiple variables. Not only the land cover classification and P-flux values, but any inherent errors on the P-export group calculations. At this point, it difficult to tease apart where most of the error lies. Better refinement of the P-flux values and higher confidence in the accuracy of land use coverages will move this project toward its goal.

Model calibration would be the next step in this process. The spatial modeling team did not have the time to perform this work, which would mainly involve the adjusting of transmission coefficients so as to force model output to agree with actual data. The work that the P-export group accomplished could be used in such a calibration; however, time did not allow for full examination of this. However, our team did conduct an initial analysis of stream flow and King County TP data to evaluate relative trends in Phosphorous loading within the watershed. Three sub-watershed were chosen to reflect near-pristine, urban, and urbanizing conditions. Graphs included in the appendix demonstrate relative trends.

In an attempt to determine the P-flux coefficient for each sub-watershed, we took a closer look at the stream flow and TP data supplied by King County. Average yearly P-export values (kg/ha-yr) were determined by averaging all of the data points for a given year (except those that obviously represented storm events) and dividing by the watershed area. Although it is well known that a considerable amount of P will enter a stream during storm runoff, the data that coincide to such events in this case were excluded due to the inconsistent sampling schedules. From this data, there is no way to determine the storm event duration or the amount of P transported to the stream that can be attributed to the event. In this case, the incorporation of the storm data would only serve to skew the computation of a baseline average P-export coefficient. To calculate the average yearly "dry season" and "wet season" P-export of sub-watersheds, data was averaged over the months July-September and October-June, respectively. This was done to better define seasonally differentiated baseline conditions.

Phosphorous export decreased for the Juanita Creek watershed from 1975 to 1998, reflecting a change in land-use patterns from predominantly rural (mainly dairy farms) to its current high density suburban state. The average yearly P-export value throughout the 1990's is roughly 0.18 kg/ha-yr.

Similar to Juanita Creek, Issaquah Creek area has been undergoing land use changes from rural to suburban. However, Issaquah creek has experienced a general increasing trend in P export concentrations since 1995, when a baseline average was roughly 0.02 kg/ha-yr.

Though consistent TP concentrations were seen over the entire 20 years of available data, the flow data proved to be unusable. Flows recorded at the sampling site were artificially controlled by releases at the Landsburg Dam by the City of Seattle to meet consumer water demands during summer dry periods. The resulting hydrograph is inconsistent and saw-toothed, rather than sinusoidal. Thus, the multiplication of the TP concentrations with the hydrograph produced skewed data results. Therefore, the resulting phosphorous export could not be used to describe a natural environment.

• Refinement of P-flux values and land cover classification schemes.
• Delineation of individual subwatersheds for the entire study area.
• Comparison of model results of the different watersheds from the different land use/land coverage images.
• Transmissivity calibration of model with actual P-export data.

Figure 9. Juanita Creek annual phosphorus export coefficient values. Note decreasing phosphorus export coefficient values in the 1990s, as sub-watershed becomes urbanized.

Figure 10. Issaquah Creek annual phosphorus export coefficient values. Note the U-shape in the data.

Dillon, P.J. and Kirchner, W.B. 1975. The Effects of Geology and Land Use on the Export of Phosphorus from Watersheds. Water Research. 9(2):135-148.

Hakanson, L. and Nilsson, A. 1992. Relationships Between Drainage Area Characteristics and Lake Water Quality. Environmental Geology and Water Science. 19(2):75-81.

Haygarth, P. 1997. On Phosphates in Europe. Scope Newsletter. 21:1-14.

Johnes, P., Moss, B., and Phillips, G. 1996. A Model for Conservation and Water Quality Management. Freshwater Biology. 36:451-473.

Paschal, J.E., and Sherwood, D.A. 1987. Relation of Sediment and Nutrient Loads to Watershed Characteristics and Land Use in the Otisco Lake Basin, Onondaga County, New York. USGS, Ithaca, New York.

Perkins, W.W, Welch, E.B., Frodge, J., and Hubbard, T. 1997. A Zero Degree of Freedom Total Phosphorus Model; 2. Application to Lake Sammamish, Washington. Lake and Reservoir Management. 13(2):131-141

Reckhow, K.H., Beaulac, M.N., and Simpson, J.T., 1980. Modeling Phosphorus Loading and Lake Response Under Uncertainty: A Manual and Compilation of Export Coefficients. EPA, Washington, D.C.

Saunders, William K., and Maidment, David R. A GIS Assessment of Nonpoint Source Pollution in the San Antonio-Nueces Coastal Basin.

http://civil.ce.utexas.prof/maidment/GISHydro/saunders/report.html.

Soranno, P.A., Hubler, S.L., Carpenter, S.R., and Lathrop, R.C. 1996. Phosphorus Loads to Surface Waters: A Simple Model to Account for Spatial Pattern of Land Use. Ecological Applications. 6(3):865-878.

Turner, G.W. 1997. Environmental Condition Indicators. Environmental Management. 21(4): 623-634.