Extent of Stream Burial and Relationships to Watershed Area, Topography, and Impervious Surface Area

Abstract

Stream burial—the routing of streams through culverts, pipes, and concrete lined channels, or simply paving them over—is common during urbanization, and disproportionately affects small, headwater streams. Burial undermines the physical and chemical processes governing life in streams, with consequences for water quality and quantity that may amplify from headwaters to downstream receiving waters. Knowledge of the extent of stream burial is critical for understanding cumulative impacts to stream networks, and for future decision-making allowing for urban development while protecting ecosystem function. We predicted stream burial across the urbanizing Potomac River Basin (USA) for each 10-m stream segment in the basin from medium-resolution impervious cover data and training observations obtained from high-resolution aerial photography in a GIS. Results were analyzed across a range in spatial aggregation, including counties and independent cities, small watersheds, and regular spatial grids. Stream burial was generally correlated with total impervious surface area (ISA), with areas exhibiting ISA above 30% often subject to elevated ratios of stream burial. Recurring patterns in burial predictions related to catchment area and topographic slope were also detected. We discuss these results in the context of physiographic constraints on stream location and urban development, including implications for environmental management of aquatic resources.

1. Introduction

Stream burial—the routing of streams through culverts, pipes, and concrete lined channels, or simply paving them over—is common in urbanized areas [1,2,3,4,5]. In many regions, the majority of buried reaches are headwater (1st and 2nd order) streams, including ephemeral and intermittent reaches, with removal of up to 70% of headwater stream length in some areas [1,3]. While the cumulative effects to ecosystem processes wrought by stream burial have important ramifications for regulation under the U.S. Clean Water Act [6], the extent of burial has not been assessed for most urban areas. Knowledge of the extent of burial events is critical for effective resource management, including preservation of remaining intact streams, restoration of urban streams, and assessing cumulative impacts of urbanization to water quantity and quality.

As with other forms of stream modification, burial alters the primary physical, chemical, and biological processes in headwater systems, contributing to a state of degradation commonly referred to as “urban stream syndrome” [7]. With the percentage of the world’s population living in urban areas continuing to grow [8], a greater emphasis has been placed on understanding the structure and function of urban streams, and associated impacts to human health, and that of downstream ecosystems. The identification of stream burial as a critical and pervasive driver of the “urban stream syndrome” [7], and the recognition that key questions remain in regards to ecosystem structure and function within piped ecosystems [9], has led to increased research activity, and the explicit consideration of stream burial in formulation of long-term research agendas (e.g., the Baltimore Long Term Ecological Research (LTER) site; Sujay Kaushal, personal communication). Subsequent attention has focused on extensive stream burial as a key contributor to “urban stream deserts” [4], riverless urban areas that evolved primarily due to stream burial for human development and population growth [10], and the phenomenon of convergent surface water distributions across urban areas of the United States [11], whereby there is a decrease in variance in the density of surface flow channels—streams and rivers—with increasing development intensity, largely due to burial of stream channels.

While not extensive to date, studies of stream burial have demonstrated a consistent set of interrelated ecosystem impacts, including modified flow velocities, altered carbon and nutrient inputs, and amplified nitrogen transport [9], loss of habitat and decreased nutrient subsidies [12], and barriers to dispersal of aquatic organisms [13]. More recently, Kaushal and Belt [14] have recognized stream burial as part of an “urban stream continuum”, whereby extensive engineering of headwater systems has expanded natural flowpaths (“urban karst”), leading to increased hydrologic connectivity within watersheds, thereby influencing the flux and transformation of nutrients, contaminants, and energy across both space and time. Leveraging this novel conceptual framework, subsequent research has focused largely on biogeochemical cycles within buried headwater streams, documenting significant reductions in nitrogen (N) uptake, gross primary production (GPP), and ecosystem metabolism (ER), with potential to influence watershed nutrient exports to downstream waters [15,16,17].

Research has demonstrated a consistent negative effect of increasing levels of urbanization on various indicators of stream health [7,18,19]. Most studies have relied on total imperviousness (TI; the proportion of a watershed that is covered in impervious surface) as the primary measure of urbanization impacts on freshwater ecosystems as TI is viewed as an integrative and comprehensive indicator [20] that can be readily incorporated into land use planning [12,21,22]. However, impervious cover alone has proven an insufficiently sensitive measure of river health [23], as significant aquatic assemblage degradation has been observed across a wide range of watershed imperviousness [24].

Recognition that the spatial configuration of impervious cover relative to stream channels may be an important moderator of the magnitude of stream ecosystem response to urbanization [25,26] has led to development of alternative metrics for measuring urbanization effects on stream ecosystems. Measures of effective imperviousness (EI, impervious cover directly adjacent to a stream channel [27,28]) and directly connected impervious area (DCIA, impervious surfaces that route stormwater runoff directly to streams via stormwater pipes [29]) have been shown to better integrate the multiple stressors of urban development, relative to TI. However, these methods also have shortcomings; they either fail to explicitly capture piped and concrete-lined stream channels (e.g., EI) or necessitate detailed information on stormwater conveyances to determine runoff routing and specific on-lot drainage patterns (DCIA). Most importantly, neither approach directly quantifies the impact of urbanization on stream habitat, and instead relies on indirect measures such as changes in sedimentation and hydrology.

Stream habitat is most directly impacted when impervious surface completely covers the stream channel. Previous research has shown a relationship between stream size and the probability of burial, with the smallest, headwater streams of the urbanized Gunpowder-Patapsco watershed, Maryland, USA, exhibiting disproportionately high rates of burial in relation to larger streams [2]. Whether this pattern remains consistent in other watersheds or across broader geographic scales is unknown. Local topographic patterns, such as slope, are also known to affect the probability of urbanization [30,31], by making some places inaccessible or unstable for building [32]. It remains unclear whether these same physiographic constraints may limit, or necessitate, the burial of streams. Insights into both the spatial and temporal patterns of stream burial, particularly with respect to stream size and topographic slope, are critical for gauging the effectiveness of land-use policies meant to foster development, while protecting the health of stream ecosystems. Historical patterns of stream burial also provide insight into the characteristics of streams that remain, information that is potentially useful for describing and understanding patterns of watershed exports, ecosystem functions, and remaining aquatic biodiversity [33,34].

Recent advancements in stream mapping, remote-sensing of impervious cover, and predictive models now make it possible to predict stream burial at a relatively high level of detail and accuracy across large areas. To enhance our understanding of the phenomenon of stream burial, we developed a novel analytical approach using improved headwater stream maps [35], moderate resolution impervious cover data [36], and recursive partitioning models [37] to map the extent and magnitude of burial across an urban gradient in the Mid-Atlantic United States. We then analyze the spatial patterns of stream burial within the context of watershed area, topography, and impervious surface area. We expect these new burial maps and analysis of burial response to physiographic parameters to provide a fresh perspective into land use decision-making processes, and the development pressures facing critical headwater stream ecosystems, past, present, and future.

2. Materials and Methods

2.1. Study Area

The study region spans the Potomac River Basin (PRB; Figure 1), which is the second largest tributary to the Chesapeake Bay (supplying 17% of the Bay’s mean annual freshwater input), and the fourth largest river along the U.S. Atlantic Coast [30,38]. Draining an area of approximately 38,000 km² across five geological provinces (Appalachian Plateau, Ridge and Valley, Blue Ridge, Piedmont and Coastal Plain), the PRB constitutes a landscape continuum from mountains to sea characterized by gradients in topography and climate that support an array of aquatic ecosystems and a diverse biota [38]. Spanning portions of 4 states (Maryland, Pennsylvania, Virginia, and West Virginia) and the entire District of Columbia, the PRB is also home to more than 5.3 million people [39]. The long history of cultural and economic development in the basin has rendered a diversity of land uses, ranging from rural forested and agricultural landscapes, to high-density, completely urbanized municipalities. The distinctive physiographic, ecological, and socio-economic characteristics of the PRB make it a model system for studying the cumulative impacts of stream burial on aquatic ecosystem function in urbanizing landscapes.

2.2. Mapping Potential Stream Burial

Our predictions of stream burial for the entire PRB used recently generated “potential” stream maps (10-m resolution, [35], henceforth referred to as “streams” or “stream maps”), which represent a much more complete coverage of streams (including ephemeral and intermittent streams) than is represented by the more commonly used National Hydrography flow lines. The term “potential streams” arises from the fact that these stream maps were generated from topographic modeling and trained from observations of streams in forested settings. Therefore, they show our best understanding of where streams would be if the entire basin had the same land use history as modern forested lands. Therefore, the stream maps provided by Elmore et al. [35] allow identification of the potential stream network in areas that have been urbanized for many decades, including those streams that were buried prior to any available stream maps. Details on the methodology used to generate potential stream maps are available, including a detailed analysis of classification omission and commission error [35].

A remote-sensing derived impervious surface area (ISA) product was acquired from the 2006 National Land Cover Dataset [36] and was used to represent hard, urban surfaces. The overall approach we used was to identify buried streams as stream segments that intersect ISA above a defined threshold. The purpose of predictive modeling (next section) is to identify the threshold that most accurately separates the training data into buried and intact stream segments. In practice, the stream and ISA layers are created at different resolutions and geo-location accuracy. It is entirely possible for a stream to appear to flow through a high-ISA pixel, but to in fact be protected by a riparian buffer that is represented by an area of low ISA pixels a short distance away. To address these and similar spatial co-registration issues between stream and impervious datasets, we calculated a selection of statistics intended to represent the spatial variability (e.g., mean, standard deviation, range, minimum, maximum) of ISA surrounding each 10 m stream segment. The results of these statistical calculations were associated with each stream segment for use in subsequent burial modeling. Additionally, the flow accumulation area (FAC, discussed below) for each segment was calculated for use with predictive modeling, which serves as a proxy for stream size.

To build a training dataset for mapping stream burial, representative, high spatial resolution (30-cm) aerial photographs (circa 2001, 2006), were obtained for 17 sample areas (mean area ~50 km², sum >2% basin area) across the PRB, chosen to capture a diversity of environments, including the full gradient in land use from forest to urban. From the stream segments located within the bounds of the aerial photos, a stratified-random set of reaches (n = 1620) was selected, with 10% of the points located within each 10% increment of ISA. Additional streams in the highest ISA category (90%–100%) were selected to ensure the training data included roughly equal numbers of buried and intact streams. The aerial photos were examined to determine the status of each selected stream reach as either buried or intact. If a stream was visible in the air photograph the stream was labeled as intact. If the stream flowed through a forested area, but the trees obscured the stream, the stream was also labeled as intact. However, if the stream was not visible due to urban land cover (pavement, buildings, or graded bare earth) the stream was labeled as buried. This technique is not sensitive to other types of stream burial in urban environments (e.g., buried streams piped and covered by vegetation under lawns or agricultural lands, or inundated by man-made impoundments). When confronted by these situations, streams were designated as intact, likely leading to a slight underestimate of the total length of affected stream across the study area.

A recursive partitioning decision tree (‘Party’ package [37]) within the R statistical programing language [40] was built using the training data (buried or intact) as the dependent variable. The ISA neighborhood statistics for the 8-neighboring pixels to each stream pixel and the flow accumulation area were used as independent variables. The fraction of buried streams in each terminal node of the resulting decision tree was used as the probability of burial for every 10 m stream segment matching the ISA statistics in that node. Each split in the tree structure was determined by conditional inference and maximized the explanation of variance in the data (p < 0.05). Comprehensive maps of each ISA statistic used in the decision tree were then used to project stream burial across the entire study region.

2.3. Accuracy Analysis of Burial Classification

Accuracy assessment of the burial probability classification was performed using components of the “party” package, and a series of custom functions, in R [40]. Two sets of data (training and validation) were created from the observations of stream burial developed from 17 areas across the PRB (described in detail, above), consisting of a stream status class (intact or buried) and impervious cover statistics for 2216 stream reaches. First, a series of 100 recursive partitioning decision trees were generated by iteratively and randomly selecting 70% of the data to train each model, with the remaining 30% of the data held in reserve for model validation. Modeled accuracy of the trees was compiled and averaged for an overall accuracy score for the burial probability classification. After accuracy statistics were calculated and it was established that the final tree structure did not depend strongly on the subset of training data used, a second accuracy assessment was conducted that utilized the same functions to create a single decision tree, using the full dataset for model training and validation. A receiver operating characteristic (ROC) analysis was then performed on this classification tree to identify the burial probability threshold that minimized false positive and maximized true positive fractions resulting from the classification.

2.4. Analysis of Stream Burial Patterns across the PRB, in Relation to Slope and Catchment Area

The completed stream burial probability layer was imported into a GIS [41], and stream burial predictions were evaluated in the context of urbanization intensity (total impervious surface area in the 8-pixel neighborhood), catchment area (FAC; flow accumulation area), and local topographic slope. Counties and independent cities, the jurisdictional level where most development decisions are made, were used to represent the largest analysis units. Urban development across the PRB largely began in and around Washington, DC, and has since radiated outward in concentric rings from the urban center [42]. A series of 16 counties and independent cities were selected to capture the full development gradient, from rural to densely urban (Figure 2), and their boundaries (National Boundary Dataset 2014; http://nationalmap.gov/boundaries.html) used to quantify the extent of predicted burial and total impervious cover [36] within each. For the purposes of comparison, counties were assigned to one of four “development tiers”, based on a combination of their level of total impervious cover (%) as of 2006 [36] and relative distance from Washington, DC (Figure 2). Mean values of impervious cover area for development tiers 1–4, were 32.9%, 13.9%, 3.0%, and 2.4%, respectively.

Flow direction and flow accumulation (FAC) were calculated for each county from 10-m National Elevation Data [43] using the TauDEM [44] suite of tools. Values for FAC represent the areal sum of all 100 m² pixels upstream from each stream pixel. Slope values for each stream pixel were derived directly from the DEM using the average landscape gradient (rise/run) in a nine-cell window surrounding and including each stream pixel.

Data for county-level burial rates, catchment area, and slope were examined to determine the relationship between burial and each of the two other variables. For comparisons with predicted stream burial, catchment area and slope were each individually treated as the independent variable. To summarize the large variability in predicted stream burial rates, stream reaches were grouped into ‘bins’ of equally spaced catchment area and the mean stream burial was calculated for streams in each bin. Mean predicted stream burial was then plotted against catchment area and slope and the resulting relationships were analyzed. To understand collinearity in the independent variables (catchment area and slope), we also analyzed slope against binned catchment area.

Initial examination of the plotted relationships revealed a consistent pattern between average predicted burial rates and catchment area, with burial exhibiting a characteristic “hump-shaped” distribution (e.g., Figure 3). For each county, a local burial ‘maximum’ is evident in the mid-range of catchment area, with a characteristic ‘ascending limb’ and ‘descending limb’ in the burial distribution to either side. To analyze how the shape of the distribution of values on either side of this maximum changed with the level of development, we subset the data, and ran separate regressions for each limb. High rates of collinearity between slope and catchment area across all levels of development precluded reliable analysis of the interaction between the two, so separate analyses were run to determine the relationship between predicted burial rates and the two independent variables individually within each limb. Due to similarities in development history and resulting patterns in burial predictions, and to ensure sample sizes for analysis across catchment area bins, the 5 independent cities of Virginia were analyzed as a single unit. All statistical analyses were performed in the statistical programming package R [40].

To investigate the potential effects of analysis scale on the relationship between stream burial and impervious surface area (ISA), a series of 3 additional analysis units covering the study area were devised at varying scales smaller than counties (Figure 4): (1) Subwatersheds; (2) 45 km² grid cells; and (3) 22.5 km² grid cells. Subwatersheds consisted of 12-digit hydrologic catalog units (HUCs) from the USGS Water Boundary Dataset (http://nhd.usgs.gov/wbd.html), and were selected as a representation of naturally-derived units. Scales for the two, grid-based analysis units were determined by halving the average area of the subwatersheds (~90 km²) to create 45 km² units, and again, halving the area of these units to create finer-scaled, 22.5 km² units. Grid cell layers were generated in ArcGIS 10.1 [41] at the extent of the county layer, and along with the subwatershed polygons, were overlaid with 2006 burial estimates and ISA data (as described above) to calculate burial and impervious cover statistics at each scale of analysis.

3. Results

Average predicted stream burial for all streams within the counties and independent cities of the study area ranged from 1.6% to 51.1% (μ = 20.5%;

Table 1. Area, stream length, impervious cover and stream burial statistics for the 16 study units. Development ‘tier’ designations correspond with those shown in Figure 2, above.

	Counties and Independent Cities	Area (km2)	Potential Stream Length (km)	Impervious Cover (%)	Stream Burial (%)
Tier 1	Washington, DC	177.0	325.8	37.6	47.3
	Arlington County, VA	67.3	170.9	33.6	39.4
	City of Alexandria, VA	39.9	95.7	42.8	51.1
	City of Fairfax, VA	16.3	40.9	28.9	31.3
	City of Falls Church, VA	5.1	11.2	26.6	32.4
	City of Manassas, VA	25.8	49.1	32.0	31.7
	City of Manassas Park, VA	6.5	16.4	29.1	22.9
Tier 2	Prince George’s County, MD	1291.0	2984.5	16.4	18.5
	Montgomery County, MD	1313.5	3830.8	10.2	11.0
	Fairfax County, VA	461.5	2803.9	15.2	14.3
Tier 3	Frederick County, MD	1728.4	4146.2	3.3	4.4
	Loudon County, VA	1053.1	3555.3	5.5	4.9
	Prince William County, VA	902.3	2566.7	8.0	7.5
Tier 4	Washington County, MD	1210.9	2151.0	3.7	4.4
	Jefferson County, WV	548.0	733.0	2.5	4.5
	Clarke County, VA	1349.6	658.4	0.9	1.6

). The City of Alexandria, Virginia, exhibited the greatest predicted burial (51.1%), followed closely by Washington, DC (47.3%), and then Arlington County, Virginia (39.4%). Predicted burial rates generally decrease with distance from the Washington, DC urban center (Figure 5), with the exception of the independent cities of Virginia. Clarke County, Virginia, by far the most rural and relatively undeveloped county, has experienced the least amount of stream burial, with only 1.6% of its streams predicted as buried in 2006.

3.1. Burial Prediction Accuracy

The average of all iterative validation models, each reserving a different 30% of the training data for validation, resulted in prediction accuracy scores of 92.7% (87%–98%) and 55.8% (36%–72%) for intact and buried stream segments, respectively, with a mean model accuracy of 83.1% (

Table 2. Accuracy values for iterative and full burial models.

Model	Accuracy (%)
	Intact	Buried	Overall
Iterative	92.7	55.8	83.1
Full	88.0	71.0	83.0

). Accuracy scores for the full model, utilizing all available data for both training and accuracy, were 88% and 71% for intact and buried segments, respectively, with a similar overall accuracy of 83%. A plot of the receiver operating characteristic identified a burial probability threshold of 0.35 as the optimum level that minimizes false positive predictions while maximizing true positive predictions.

3.2. Relationship between Predicted Burial Rates and Impervious Cover

Predicted, county-level stream burial increases linearly with total impervious cover (Figure 6), in an approximately 1:1 relationship (Y = 1.135 × X − 0.572; R² = 0.99; p < 0.001), with a few exceptions (i.e., deviations from the 1:1 line). Several Tier 1 counties and independent cities exhibited a higher proportion of stream burial than expected based on a 1:1 relationship with impervious cover. The largest discrepancies exist in the oldest municipalities with the greatest levels of urban development, including Washington, DC (1.26), the independent cities of Falls Church, VA (1.22), and Alexandria, VA (1.19), and Arlington County, VA (1.17). Conversely, counties that developed more recently, but still exhibit high levels of development (mostly Tier 2), were found to have lower levels of predicted burial than expected based on the 1:1 relationship. The Virginia independent cities follow the same pattern as counties, where cities that developed earlier exhibited higher predicted rates of burial per unit of impervious cover than did the more recently developed cities. All the Virginia cities exhibit very high predicted burial rates per unit imperviousness, with the exception of Manassas, VA, which falls directly on the 1:1 line, and Manassas Park, which is the only independent city to show a significantly lower burial to impervious ratio (0.79).

3.3. Effects of Scale on Burial/ISA Relationships

Relationships between percent stream burial and ISA for the 4 scales of analysis (

Table 3. Statistics for four scales of analysis unit across the study area.

Analysis Units	N	Mean Area (Range)	Area with >30% ISA
Counties and Independent Cities	16	632.1 km2 (5.19–1728.50)	292.01 km2
Subwatersheds	166	89.9 km2 (33.75–170.44)	221.48 km2
45 km2 grid cells	283	45 km2	450.24 km2
22.5 km2 grid cells	534	22.5 km2	517.41 km2

) are shown in Figure 7. In each case, the data exhibit a roughly 1:1 relationship below a level of approximately 30% ISA, above which most, if not all analysis units exhibit a greater amount of stream burial per unit ISA. This general relationship between burial and ISA is consistent for all analysis units, regardless of scale, across the full range in ISA cover, though not surprisingly, analyses with smaller-scaled units yield a larger number of units above the 30% ISA barrier.

At the county-level, units above 30% ISA (Figure 7A) include Washington, DC, Arlington County, and the independent cities of Alexandria and Manassas, Virginia. These regions of the study area (Figure 8) have the longest history of urban development, spanning back to mid-late 18th century, and have experienced continuous development and redevelopment to the present time. Of these units, only Manassas has maintained a roughly 1:1 ratio of burial to ISA, while the others exhibit significantly higher levels of burial per unit ISA. Only two subwatersheds (Figure 7B) exhibit greater than 30% ISA, their area corresponding with the most densely-developed regions of Washington, DC, and its first tier suburbs (Figure 8). In general, burial-ISA relationships for subwatersheds vary little from 1:1.

At an analysis scale of 45 km² units, those grid cells with the highest ratio of burial to ISA (Figure 7C, orange ovals) still center on Arlington County, Alexandria, and the heart of Washington, DC (Figure 8). In addition, exhibiting >30 ISA, and elevated burial/ISA ratios at this analysis scale are significant areas of the first-tier DC suburbs in Maryland and Virginia. Interestingly, there are two grid units with just over 30% ISA (Figure 7C, green oval) that exhibit significantly lower ratios of burial to ISA, centering on the cities of Herndon and Sterling in the western Suburbs of Virginia, and in suburban Prince George’s County, Maryland, to the southeast of Washington, DC These areas are characterized by dense residential and commercial development centered on the local transportation corridors of VA Route 28 and the Dulles Airport Access Road (Rt. 267), and on the Interstate 95/495 corridors, respectively (Figure 8).

Finally, at the smallest analysis scale of 22.5 km², the greatest number of analysis units exhibit increased burial/ISA ratios above the threshold of 30% ISA (Figure 7D). Those with the greatest deviation from 1:1 again include grid units spanning areas of Arlington, Alexandria, and the core of Washington, DC (orange ovals). There is also a large cluster of units above the 1:1 line, consisting of additional portions of DC and close suburbs with a longer history of development and infill. Those cells falling just above, and on the 1:1 line include large portions of the Potomac River, and other park areas with significant undeveloped area. Several additional areas in the suburbs of Maryland and Virginia appear at this analysis scale, largely centered on more recent suburban and commercial development around significant transportation corridors (Figure 8), and exhibiting a wide range in deviation from the 1:1 line.

One grid cell, centering on Hagerstown, MD (Figure 8), falls just above the 30% ISA threshold, with recent and rapid development around the intersection of Interstates 70 and 81 leading to a very high ratio of burial per unit ISA (Figure 7D, yellow oval). Additional cells, including Springfield, VA (I-95 corridor), Manassas, VA (I-66 corridor), and in Prince George’s County, MD (I-495 corridor), fall on the 1: line, having managed to maintain significant natural land cover (Figure 8) despite centering on significant transportation corridors. A grid cell in Prince William County, VA, adjacent to Manassas, and one cell centered on the I-270 corridor in Montgomery County, MD, exhibit burial/ISA ratios just below 1:1. Two additional grid cells with >30% ISA fall well below the 1:1 line (Figure 7D, green oval), and include the areas around Herndon, VA, and suburban Prince George’s County, MD, that exhibited similarly low ratios at the 45 km² scale, as discussed above.

3.4. Relationship between Burial and Catchment Area

Stream burial-catchment area relationships exhibit a generally consistent pattern across all counties, regardless of development stage (Figure 9 and Figure 10). Predicted burial rates are elevated for the very smallest streams, decreasing sharply to a local minimum (corresponding with maximum slope values, discussed below) around catchment areas of 0.1–0.2 log₁₀ ha, after which they increase linearly to a maximum burial rate in the mid-range of catchment area (~0.8–2.5 log₁₀ ha, depending on development stage). After this local maximum, burial rates generally decrease with increasing catchment area, until approximately 3.5 log₁₀ ha, where predicted burial rates tend to level off and remain very low (~zero) as catchment area continues to increase.

The magnitude of the maximum predicted burial rate, and thereby the slope of the distribution of burial data on either side of the local maximum, becomes steeper (more positive to the left, and more negative to the right of the maximum) with increasing urban development (

Table 4. Characteristics of the distribution of burial prediction in relation to catchment area for the sixteen study counties. Results include analyses for both the ascending and descending limbs of the burial distribution. Results with daggers (^†) are for models where the slope of the ascending or descending limbs was not significant.

Counties		Max Burial % (FAC)	Local Minimum (FAC)	Local Maximum (FAC)	Ascending Slope	Ascending Intercept	Descending Slope	Descending Intercept
Tier 1	Washington, DC	79.4 (3.0)	0.1	1.2	37.20	10.10	−9.28	65.21
	Arlington County, VA	47.4 (3.5)	0.1	1.3	21.75	18.34	−13.30	65.60
	Independent Cities, VA	78.0 (3.4)	0.2	0.8	36.09	17.82	2.48 †	34.96 †
Tier 2	Prince George’s County, MD	15.0 (1.3)	0.2	1.2	7.36	6.90	−3.71	19.59
	Montgomery County, MD	12.9 (0.8)	0.2	0.8	11.06	4.75	−4.48	16.88
	Fairfax County, VA	18.6 (1.3)	0.25	1.3	9.95	6.64	−6.93	26.89
Tier 3	Frederick County, MD	9.1 (3.5)	0.35	1.4	0.89	3.35	0.22 †	4.10 †
	Loudon County, VA	7.3 (0.0)	0.25	1.3	3.61	1.88	−1.62	8.30
	Prince William County, VA	9.0 (1.0)	0.15	1.0	5.25	4.34	−2.63	11.83
Tier 4	Washington County, MD	8.2 (2.3)	0.1	2.2	1.74	1.42	−2.90	13.46
	Jefferson County, WV	12.6 (0.1)	0.2	0.8	1.78 †	3.19 †	−0.08 †	4.45 †
	Clarke County, VA	3.8 (2.2)	0.7	2.1	0.85	−0.01	−1.12 †	5.16 †

, Figure 9 and Figure 10). In the cities and counties with the highest levels of development (e.g., Tier 1, Figure 9), rather than decreasing sharply as catchment area increases past the initial maximum rate, predicted burial remains high across a wide range, before decreasing sharply at catchment areas approaching 3.5 log₁₀ ha. The maximum predicted burial generally occurs at approximately 1.0 log₁₀ ha (range 0.8–1.4 log₁₀ ha, with the exception of two counties—Washington, MD and Clarke, VA—where the maximum is shifted far to the right (~2.5 log₁₀ ha), and Clarke County, where maximum predicted burial occurs at 2.1 log₁₀ ha catchment area.

The relationship between predicted burial rates and catchment area was highly significant for the majority of counties, and across both ascending and descending limbs of the burial curve (

Table 5. Results from linear models consisting of Burial~FAC, and Burial~SLOPE.

Counties		Limb	FAC	R2	SST	SLOPE	R2	SST
Tier 1	Washington, DC	Ascending	***	0.91	2183.87	***	0.95	2183.9
		Descending	*	0.18	5624.6	---	-	-
	Arlington County, VA	Ascending	***	0.65	1335.17	***	0.96	1335.17
		Descending	**	0.38	4672.4	***	0.62	4672.4
	Independent Cities, VA	Ascending	**	0.81	450.05	---	-	-
		Descending	---	-	-	---	-	-
Tier 2	Prince George’s County, MD	Ascending	***	0.86	69.312	***	0.95	69.313
		Descending	***	0.75	210.571	***	0.71	210.57
	Montgomery County, MD	Ascending	***	0.96	35.493	***	0.99	35.492
		Descending	***	0.97	379.15	***	0.87	379.16
	Fairfax County, VA	Ascending	***	0.86	126.375	***	0.99	126.374
		Descending	***	0.96	503.26	**	0.45	503.27
Tier 3	Frederick County, MD	Ascending	***	0.86	1.01895	---	-	-
		Descending	---	-	-	---	-	-
	Loudon County, VA	Ascending	***	0.89	16.044	***	0.96	16.0435
		Descending	***	0.68	39.097	***	0.56	39.097
	Prince William County, VA	Ascending	***	0.92	18.042	***	0.98	18.0419
		Descending	***	0.88	114.503	***	0.60	114.502
Tier 4	Washington County, MD	Ascending	***	0.81	33.0171	***	0.64	33.017
		Descending	**	0.46	41.182	---	-	-
	Jefferson County, WV	Ascending	---	-	-	---	-	-
		Descending	---	-	-	---	-	-
	Clarke County, VA	Ascending	***	0.60	3.4154	---	-	-
		Descending	*	0.35	10.0001	---	-	-

Notes: Significance codes: 0 ‘***’, 0.001 ‘**’, 0.01 ‘*’, no significance ‘---’.

). In some cases, e.g., the Independent Cities of Virginia and for Frederick County, MD, catchment area was not significantly correlated with predicted burial rates for the descending limb of the burial distribution. Catchment area did not explain predicted burial rates at any point in the distribution for Jefferson County, West Virginia.

3.5. Relationship between Burial and Slope

The very smallest watersheds in all counties exhibit moderate levels of predicted burial (compared to predicted burial calculated across larger watersheds), and these areas correspond with low to moderate slope values of up to 7–8 degrees. Conceptually, these are relatively flat areas located on tops of hills supporting only the smallest of headwater streams. Predicted burial rates decrease sharply to a local minimum as slope increases to its maximum (range = 6.2–24.16 degrees, μ = 9.47, Figure 9 and Figure 10), which occurs in all counties at approximately 0.2 log₁₀ ha catchment area. Conceptually, these are areas just downstream of the channel head where stream gradients are largest. Predicted burial rates then increase sharply to their maximum, as slope values decrease from their maximum, with maximum burial rates occurring in the range of 2.37–6.36 (μ =3.38) degrees of slope across all counties. Beyond the point of maximum burial rate, both predicted burial and slope decrease linearly until both approach values of 0 at approximately 3.5 log₁₀ ha catchment area (Figure 9 and Figure 10).

The degree to which slope explained the variability in predicted burial rates across counties and development tiers was mixed, except for the Tier 2 counties, where slope was significant for both ascending and descending limbs in all counties (Table 5). Where slope was significant, it was highly so (p-values < 0.001), but there were several counties across the development gradient for which slope was not related to predicted rates of burial for either the ascending or descending limbs, including the Virginia Independent Cities (Tier 1), Frederick County, MD (Tier 3), and both Jefferson County, WV, and Clarke County, VA, in Tier 4). Slope was a significant factor in predicted burial rates for the ascending burial limbs in both Washington, DC, and Washington County, MD, but insignificant for the descending limbs in these two areas.

4. Discussion

Stream burial is an emerging environmental issue associated with urbanization locally, regionally, and globally [2,3,4]. This paper suggests that headwater stream burial is strongly related to watershed characteristics including topography, catchment area, and impervious surface cover, with repeated, consistent patterns developing across both space and time. Our findings improve our understanding of: (1) spatially explicit patters of burial from local to regional scales; (2) evolution of drainage networks with urbanization; (3) identification and recognition of the widespread occurrence of buried streams as eco-hydrological units within landscapes; and (4) an improved understanding of evaluating land use/management decisions along sociopolitical boundaries to prevent stream burial. We discuss some of these aspects in further detail below.

4.1. Relationship of Predicted Stream Burial to Impervious Surface Area (ISA) across Multiple Scales

In general, it appears that streams are buried in proportion to ISA when ISA is below 30%. Above this threshold, it is common to find areas where streams are buried at a higher rate per unit ISA. This phenomenon is most likely due to the lack of remaining developable land that is not directly adjacent to streams, with continued development (“infilling”) leading to stream impacts. For example, the ratio of predicted burial to ISA for any given analysis unit (county, watershed, or regularly-spaced grid cell) did not exhibit a dependence on the spatial area of the analysis unit. Initial observations at the county-level (Figure 6) suggested an approximate 1:1 relationship between predicted burial and ISA, until an ISA level of roughly 30% is achieved (Figure 7A). County units with greater than 30% ISA exhibited greater ratios of burial/ISA, and included 5 independent cities in Virginia, Arlington County, VA, and Washington, DC We hypothesize that this is due principally to the long time since development in these areas (since the mid-late-18th century), and their long term roles as commercial and government centers, leading to a greater degree of predicted stream burial through sustained and dense development during a period when there were few governmental policies constraining development.

Observations suggest areas with elevated burial/ISA ratios tend to be older developments, either in the DC core, or near-suburbs that have been developing and redeveloping for decades, though there are a couple of newer-developing areas with elevated ratios (e.g., Hagerstown, MD). Areas above 30% ISA that fall along the 1:1 line include a mixture of both old and new development, and tend to be areas that include significant non-urban cover areas, such as parks or agriculture, or are centered on large bodies of water such as the tidal portions of the Potomac River (Figure 8). Areas with lower burial/ISA ratios tend to be new suburban and commercial centers, including those surrounding important transportation hubs and corridors. Whether this is due to policies more effective at preventing stream burial, or whether they have simply not developed to their full density, is as yet unclear.

A larger fraction of small landscape units generally exhibit ISA values above the 30% ISA threshold (Table 5). The simple explanation is that because of the larger number and smaller size of the analysis units, the units are more efficient in capturing portions of the study area with a higher percentage of ISA, and therefore potential stream burial. Nevertheless, a value of 30% ISA, which corresponds with an elevated burial/ISA ratio at all analysis scales, appears to be a useful threshold value for separating areas that seem to have slightly different behavior.

4.2. Predicted Burial Patterns by Development Stage, and Relation to Catchment Area and Slope

Headwater streams across the study area are predicted to be buried more extensively than are larger streams at all levels of urban development, with maximum predicted burial rates for all but the most urbanized counties peaking at catchment areas of approximately 1.0 log₁₀ ha, after which burial rates decrease linearly with increasing catchment area. In the more urban counties, streams are predicted as buried at high rates (spatially distributed evenly) between approximately 1.0 and 2.5 log₁₀ ha catchment areas, after which predicted burial rates decrease fairly rapidly as catchment size increases. Both of these findings are consistent with those of Elmore and Kaushal [2] from the Gunpowder-Patapsco watershed of Maryland, suggesting similar initial constraints to stream burial across basins and management entities. In the case of heavily urbanized areas (e.g., Tier 1, above), it appears the initial burial constraint consistent with the 1.0 log₁₀ ha catchment area is overcome. This “soft constraint” is most likely related to some aspect of the landscape affecting the ‘favorability’ of a particular area to urban development (e.g., topographic slope, discussed below), and becomes less of a constraint once an area becomes developed enough that the most favorable (e.g., low slope) lands are unavailable for further development. The upper limit to maximum predicted burial, occurring at roughly 2.5 log₁₀ ha catchment area appears to be universally observed, suggesting a ‘hard constraint’, beyond which burial is largely impossible. As suggested by Elmore and Kaushal [2], streams above this threshold may simply be too large for burial to be economically feasible, or the limit may be set by regulatory statute, such as the floodplain development restrictions imposed by the Federal Emergency Management Agency (FEMA).

Time-series data on stream burial rates and patterns are currently unavailable, therefore, our analysis utilizes a space-for-time substitution [45] across the study counties, whereby the urban-rural gradient, as indicated by total impervious cover (IC), serves as a proxy for understanding trends in stream burial over time. Space-for-time substitutions are often used in urban ecology to assess the consequences of increasing urbanization intensity [46,47,48,49,50], with more heavily urbanized areas serving as a reference for the future of other, less-developed places [51]. However, this analysis approach assumes that spatial and temporal variation across sampling units are equivalent [45]. Given the long history (200+ years) of urban development in the region, it is unlikely that these assumptions are fully upheld; in other words, the relationship between stream burial, catchment area, and slope has almost certainly changed over time. Our results adequately incorporate these changes (and were the subject of our discussion), but do not speak to the potential for change in the future and therefore our results should not be used to make future predictions.

4.3. Mapping Approach and Uncertainties

The combination of new methods developed here, including improved stream maps, moderate resolution impervious area maps, and non-parametric modeling allowed for the successful mapping of burial probability for every stream reach within the PRB, overcoming obstacles presented by previous mapping attempts requiring intensive, on-the-ground investigations [3]. There are several instances, however, where error may have been introduced into the burial classification. For example, the improved stream maps of Elmore et al. [35] were modeled using training data gathered from forested watersheds as a reference, and using digital elevation models to place streams where they should occur naturally, based on the topographic flow of water across the landscape. These methods may be less accurate in heavily urbanized landscapes, where land cover and topographic relief have been highly altered, and streams may have never existed or been diverted significantly from their original course. In many of these areas, there may be no maps to confirm the pre-development status and location of streams, thereby preventing confirmation of their previous existence, or lack thereof.

The methods for calculating the probability of burial for stream reaches cannot account for slight variations in stream path caused by natural shift of channels over time, or other co-registration errors between the stream and impervious cover data layers, potentially leading to inaccuracies in the calculation of values for the average level of impervious cover in the neighborhood of each stream pixel, which was used to classify burial potential. Additionally, while creating the training data for classifying burial potential, only burial related to impervious cover was considered. In reality, other forms of burial can also be prevalent across a watershed, including streams routed into pipes, and buried under residential lawns, in agricultural fields [52], or inundated by man-made impoundments. As we were primarily interested in urban-related stream burial, these cases were not classified as being buried, leading to the probability that stream burial is underestimated. Despite these uncertainties, we see robust, repeated patterns in the data.

4.4. Future Work

While this work examines the phenomenon of stream burial to a greater degree than previous efforts, numerous open questions remain regarding the evolution and extent of burial across the landscape, and its effects on ecosystem function. Analysis of changes in stream network structure due to burial, and potential impacts to habitat connectivity for aquatic organisms is ongoing. Further examination of trends in burial distribution, and relationships with a broader range of landscape and socio-economic factors is warranted, to develop a more complete understanding of the complex interplay between topography, human economics and demographics, and decisions to bury streams. Newly available land cover data and updated remote sensing products allow for broader analysis across both time and space. Detailed study of burial dynamics across the development gradient in biomes with differing precipitation regimes would be particularly interesting, comparing stream network and ecosystem response in systems with variable natural network architecture and hydrologic regimes. At smaller scales, detailed stream burial data can inform critical analyses of specific development policies across urbanizing counties, to identify more effective protection policies for headwater streams, and guide restoration efforts to maximize positive ecosystem outcomes within heavily urbanized stream networks.

5. Conclusions

Stream burial is a spatially pervasive phenomenon, with predicted burial increasing linearly with total impervious cover across all levels of development, bringing into question the efficacy of existing stream protections. The close relationship between stream burial and total impervious surface area suggests that the two measures provide similar, and perhaps redundant, information. However, stream burial maps provide a spatially-explicit measure of potential stream-specific impacts, accounting for the effects of impervious surface area immediately adjacent to and covering stream channels, including direct habitat loss, and the probable effects of contiguous impervious cover on physical and hydrologic regimes in stream ecosystems. Predicted stream burial data could be used to identify high impact watersheds for targeting restoration, to address riparian and network connectivity issues, and to integrate effects of hydrological change into efforts to manage downstream water quality, including the Federally-regulated TMDL (Total Maximum Daily Load) process.

We know that loss of in-stream and riparian habitat holds potential implications for aquatic organisms, and their ability to move both within and between headwater systems. Future work should consider the effects of stream burial on network geometry (the size and spatial orientation of remaining stream reaches), and the effects on habitat connectivity within and between headwater systems on biodiversity patterns in aquatic communities. Future work might also include examining predicted burial rates across time, to better discern how burial has proceeded in relation to physiographic and policy constraints, and related effects to ecosystem structure and function across large, developing watersheds.

Headwater stream burial is prevalent across the study area, even within watersheds with very little urban development. Both slope and catchment area combine to limit stream burial during development, but these constraints were largely overcome in the most intensely urbanized jurisdictions. Headwater stream systems are critical to the maintenance of downstream water quality and hydrologic regimes [53], and yet, continue to be disproportionately affected relative to larger streams. This understanding might be used to justify more rigorous and uniform protection policies and other strategies to reduce the impacts of burial and preserve the ecological function of these vital ecosystems.