Home | Archives | About | Login | Submissions | Notify | Contact | Search
 E&S Home > Vol. 12, No. 2 > Art. 25

Copyright © 2007 by the author(s). Published here under license by The Resilience Alliance.
Go to the pdf version of this article

The following is the established format for referencing this article:
Pijanowski, B., D. K. Ray, A. D. Kendall, J. M. Duckles, and D. W. Hyndman. 2007. Using backcast land-use change and groundwater travel-time models to generate land-use legacy maps for watershed management. Ecology and Society 12(2): 25. [online] URL: http://www.ecologyandsociety.org/vol12/iss2/art25/

Research, part of Special Feature on Crossing Scales and Disciplines to Achieve Forest Sustainability

Using Backcast Land-Use Change and Groundwater Travel-Time Models to Generate Land-Use Legacy Maps for Watershed Management

Bryan Pijanowski 1, Deepak K. Ray 1, Anthony D. Kendall 2, Jonah M. Duckles 1 and David W. Hyndman 2

1Purdue University, 2Michigan State University


We couple two spatial-temporal models, a backcast land-use change model and a groundwater flow model, to develop what we call “land-use legacy maps.” We quantify how a land-use legacy map, created from maps of past land use and groundwater travel times, differs from a current land-use map. We show how these map differences can affect land-use planning and watershed management decisions at a variety of spatial and temporal scales. Our approach demonstrates that land-use legacy maps provide a more accurate representation of the linkage between land use/cover and current water quality compared to the current land-use map. We believe that the historical signatures of land-use impacts on current water quality should be considered in land-use planning and watershed management.

Key words: land use; legacy maps; groundwater travel times; backcast land change model


The current global rate of land-use change is unprecedented (Foley et al. 2005). Land use directly influences hydrologic processes such as evapotranspiration and overland flow (Tang et al. 2005a, Dale et al. 2005), which in turn affect the spatial and temporal distribution of effective recharge flux to groundwater (Harbor 1994, Wayland et al. 2002, Jayawickreme and Hyndman 2007). In the Great Lakes Basin, several large-scale land-use change trends over the last 150 yr have probably had significant impacts on water cycle fluxes (Changnon 1992). In the late 1800s, many of the old-growth forests in Michigan and Wisconsin were logged to provide building materials for the region’s large cities, most notably Chicago and Detroit. Much of this land was converted to agricultural use in the early 1900s with another agricultural expansion in the 1960s that caused a large loss of wetlands and forests. In the last 30 yr, many marginal agricultural lands have been converted to residential and commercial uses because of urban sprawl (Brown et al. 2005, Pijanowski et al. 2006a). Agricultural decline in this region has also led to significant aforestation in rural areas (Brown et al. 2005).

Studies of the impact of land use on surface water quality have focused mostly on how overland flow, particularly runoff, affects water biogeochemistry (Niehoff and Brownstert 2002, Fitzpatrick et al. 2007). However, during the stream base-flow conditions prevalent during dry summer months in Michigan, more than 85% of stream water volume is derived from groundwater (Boutt et al. 2001). Stream inputs from groundwater are temporally averaged signals from historic land uses. Groundwater travel times from the source of recharge water to surface water can range from days in areas close to streams to hundreds of years in regions near watershed boundaries (Freeze and Cherry 1979, Senger and Fogg 1990, Wayland et al. 2002). Groundwater travel times are influenced by (1) the distance groundwater travels from the point of recharge to the outlet, (2) the physical properties of the medium through which the water flows, and (3) the gradient in hydraulic heads in the aquifer that drive flow through the subsurface, which is affected by recharge rates and hydraulic properties.

Investigations of land use and surface water quality for the purposes of watershed management should consider the potential impact of past land uses on current water quality. Creating a map of the locations and types of land uses that influence the state of current water quality requires the use of tools that estimate the spatial distribution of past land uses in relationship to groundwater travel times. Such maps, which we call land-use legacy maps, represent the integration of complex spatial-temporal interactions of land-use change and hydrology.

The objective of this paper is to introduce the concept of land-use legacy as it is influenced by groundwater travel times. We quantify spatial differences between land-use legacy maps and current land use/cover using geographic information systems. We describe how we developed a legacy map using (1) a backcast land-use change model based on historical census data and a neural network tool, and (2) a groundwater travel-time model. Our discussion focuses on the importance of using a legacy map in watershed management and planning. We also summarize important issues surrounding the accuracy of land-use change and groundwater travel-time models used to create land-use legacy maps.


This study was conducted in the Muskegon River Watershed (MRW) located in the west-central part of the Lower Peninsula of Michigan, USA (Fig. 1). The 7057-km² MRW is currently dominated by forest in the north, agriculture in the central portion, and urban areas in the south. The watershed was heavily logged from the late 1880s through the 1890s, mostly in efforts to rebuild Chicago after the 1871 fire (Alexander 2006), and agricultural expansion occurred in the 1910–1930s and then again in the 1960s. Twenty percent of the watershed is in public ownership, most as state forest land.

Twelve counties fall within the MRW (Fig. 1). Land-use planning occurs at a subcounty, i.e., minor civil division or MCD, level of government. In the MRW, there are 124 MCDs that fall wholly or partially within the watershed. Each is required to develop land-use master plans and create ordinances and zoning laws that control land-use decisions.

The main river, the Muskegon River, runs the length of the watershed from its headwaters at two large inland lakes, Houghton Lake and Higgins Lake, to its outlet in Muskegon Lake. The city of Muskegon, the largest city in the MRW with a population of about 40,000 in 2000, is located along the southern shore of Muskegon Lake. The peak flows in this system tend to be in the spring and associated with either snowmelt or large rain events.

The average annual rainfall from 1899 through 2001 was 83 cm, with lower average precipitation of 75 cm from 1920–1940, followed by an increase to the recent average of approximately 90 cm. The central portion of the watershed is dominated by interlobate tills, whereas the southwestern and northeastern portions are predominantly sand and gravel outwash and lacustrine deposits with a few scattered areas of peat.


We used a variety of GIS-based data sets to construct the backcast land-use change and groundwater travel-time models. Details of these data and how they were processed are contained in Appendix 1.

Backcast land-use change model

Five different submodels (Fig. 2) of our backcast model were needed to create an annual time series of land-use/cover maps dating from presettlement to 2006.

The earliest portion of our backcast model is represented as a presettlement map developed from a mid-1800s General Land Office survey of vegetation, acquired in digital form from the Michigan Department of Natural Resources (MDNR). Land-cover codes from the map were placed into a modified Anderson level III classification, which was then reclassed to the major vegetation types, i.e., forests, grasslands, wetlands, open water, and barren. Land uses/covers for dates prior to 1900 are assumed to be fixed for the purposes of our analysis.

The four other submodels of our backcast land-use change model consist of specially parameterized land transformation models (LTMs). The land-use/cover maps from 1978 and 1998 were used as inputs to these submodels. We followed the procedures of Pijanowski et al. (2002a, 2005, 2006) in developing GIS inputs and using a neural network to determine how spatial drivers such as distance to city center are fit by a nonlinear model to land-use/cover change maps created using the 1998 and 1978 land-use/cover maps. The technical details of the model can be found elsewhere (Pijanowski et al. 2005, Pijanowski et al. 2006), but relevant information is provided here. LTM modeling consists of four steps: (1) training and testing of the neural network to generate a spatial ranking map of cells transitioning to urban, forest, and shrub/grassland classes; (2) using a battery of goodness-of-fit tests, e.g., the area under the receiver operating characteristics curve (see Pijanowski et al. 2006 for details), to select the best neural network model; (3) developing annual allocation tables, i.e., the number of cells that need to transition between time periods, for urban, agriculture, forest, and shrub/grasslands classes using historical data; and (4) creating the annual land-use/cover maps using the allocation tables and the spatial ranking maps. Two sets of historical data, land in farms and year-built housing information, are used as part of the allocation subroutine to adjust the amount of agriculture and urban development back in time from a baseline 1978 land-use/cover map.

The second backcast submodel is described below. For the third backcast submodel, two historical databases were used to calculate the amounts of urban and agricultural transitions between 1939 and 1978. First, the “year built” variable of the 2000 U.S. Census 100% SF-1 file database on housing statistics was used at the subcounty, i.e., minor civil division (MCD), level to determine the proportion of urban cells in each MCD that should exist during each census period. We assumed that the proportion of housing stocks existing in each census period reflected the proportion of ages for all urban categories, e.g., commercial, industrial, etc., within an MCD. U.S. Department of Agriculture historical Land In Farms (LIF) data were used to determine the proportional amounts of agricultural land use that should be allocated to each county based on our 1978 base year. We used proportional amounts so that the amount of agriculture in the 1978 MiRIS database was set to 1.0 in the LIF database; agricultural amounts for previous years were then adjusted proportional to the base amount for that county. Thus, if a county had 10,000 cells of agriculture in the 1978 MiRIS map and the LIF database difference between 1978 and 1974 showed a 20% gain in agriculture as we moved backward in time, then the model allocated 12,000 cells of agriculture in that county for 1974. Because LIF data were available at approximately 5-yr intervals from 1910 to the present, we interpolated to annual values, by county, using a cubic spline. R² = 1.0 for all cubic spline models, so true values from the census were not adjusted to fit the trend.

For the second and third backcast submodels, a specially parameterized LTM was used to create a spatial ranking map for areas likely to transition at any given time step moving backwards from 1978 to 1900. To backcast urban areas, we applied spatial ranking maps to areas that were urban in 1978 and used the year-built statistic to determine the proportion of cells to take away from the urban class as we moved back in time. We used a similar procedure to develop spatial ranking maps for forest and shrub transitions.

The second backcast submodel, spanning 1900 to 1939, uses two different extrapolations. Annual allocations of urban areas were determined using a linear extrapolation of housing data from 1950 to 1900 because housing data prior to 1939 are lumped into a single category by the U.S. Census. For agricultural allocation, LIF data are reported back only to 1910, so a cubic spline of LIF was used to extrapolate agricultural land-use amounts by county back to 1900. The spatial ranking maps for the 1978 and 1998 training sessions were used to create annual maps between 1900 and 1939 as part of submodel 2.

For our fourth backcast submodel, we used the neural net model from the third submodel to assign spatial rankings from 1978 to 1998 in known areas of urban and forest transitions. For the allocation routine of the LTM, we assumed a linear rate of transition of urban areas and new forest to create annual land-use/cover maps between 1978 and 1998.

Finally, our fifth submodel addresses the need to develop annual maps from 1998 forward to 2006. This submodel was developed using a separate neural net model (cf. Pijanowski et al. 2005, 2006) using 1978 to 1998 change in urban area and forest to create the spatial ranking map of post-1998 change. Population estimates from the Michigan Historical, Arts and Libraries for Muskegon River Watershed (MRW) counties between 2000 and 2005 were used to estimate annual urban change (see Pijanowski et al. 2002b for details). Pijanowski (2006) found that the aforestation rate was 64% of the urbanization rate in the MRW between 1978 and 1998. New forests were thus added in our 1999 to 2006 forecast submodel in an amount that was 64% of new urban areas.

For all submodels, we assumed that land-use/cover transitions did not occur in public lands and that, once urbanization occurred, that location remained urban. In addition, we also made the assumption that shrub/grassland can remain in that category for only 5 yr, after which it transitions to forest if it does not urbanize.

Groundwater travel-time model

The groundwater travel-time (GWTT) model is based on a MODFLOW-2000 simulation of the MRW and surrounding region. To avoid edge effects, the geographic domain of this model was extended beyond the MRW to significant hydrologic barriers such as large lakes or rivers. The active domain of the model is described using approximately 3.34 million cells that are 106.3 m on a side, with two layers in the vertical direction. Thus, 16 cells of the backcast land-use change model fit into one cell of the GWTT model.

Data from several thousand oil and gas wells were interpolated to provide a bedrock bottom for the model that was assumed to be impermeable. The geologic materials within the model were assigned using digital maps of quaternary sediments (Farrand and Bell 1982). Hydraulic conductivities for these zones were calibrated to minimize both the differences between modeled and observed water-table levels at thousands of residential wells throughout the region and modeled and observed stream flows during late-summer base-flow conditions. We used the groundwater model to delineate the contributing groundwater to the surface water bodies in the MRW; the coverage of the groundwater watershed differs from that of the surface water watershed (Boutt et al. 2001).

As part of a larger set of research projects, approximately 10 stream gages were in operation from 2002 through 2004 in the MRW along with five long-term flow records from USGS gages. Approximately 150 streamflow measurements collected across the watershed during base-flow periods in 2003 and 2004, along with thousands of water-table elevations collected at the time of drinking-water well installations, were used to calibrate the groundwater models for the MRW.

The extended watershed model is able to identify regions of the landscape that are not within the surface watershed of a river or stream, but still contribute to streamflow via groundwater flow paths. We call this region the “groundwatershed,” a distinction that can be important in identifying bodies of water in which the management of water quality is of significant concern. In some areas with significant topography or lateral variability of hydraulic conductivity, the regions that are only in the groundwatershed can contribute significant water and solute fluxes to streams and lakes (Boutt et al. 2001, Hyndman et al. 2007).

The steady-state simulated water-table levels were then imported into ArcGIS along with a raster grid of the hydraulic conductivities within the groundwater model. Groundwater flow directions were then calculated using the ArcGIS D8 FLOWDIRECTION algorithm. The groundwater flow velocity of each cell was calculated using Darcy’s law, as the product of the hydraulic conductivity and the gradient of the water table, divided by the porosity of the sediments within the cell. An estimate of 30% porosity was applied for the entire model because no direct porosity measurements were available. Given the flow-direction grid, the travel-time grid was calculated using the ArcGIS FLOWLENGTH algorithm (cf. Maidment 2002), cost-weighted by the inverse of the groundwater flow velocity.

The water table is a dynamic feature that responds to seasonal cycles and long-term trends in recharge. Here we have assumed that a steady-state condition is an adequate representation of the true groundwater travel times and that there is no spatial variability in recharge, which was fixed at 30.5 cm/yr. The simplicity of this GWTT model was the primary motivation for this assumption because it allows rapid parameter sensitivity analysis and scenario evaluation. In addition, simulated steady-state transport times are expected to provide a reasonable representation of the average transport time through these systems for a multidecadal simulation, because transient water-table responses to seasonal recharge variations have been shown to be on the order of three years for this watershed (Kendall and Hyndman 2007).

Creating legacy maps

A C program was written to use the map containing groundwater travel times in years to select the land-use/cover class for each location in the historical time series of land-use/cover maps generated from the backcast model, and to place that land-use/cover class into the legacy map. Thus, if a cell requires 35 yr for water to travel from the recharge point to a surface water body, then the land-use/cover class assigned to that location 35 yr ago, i.e., 2006 - 35 = 1971, is placed in the legacy map. The resultant legacy map then contains land-use/cover classes during years influenced by groundwater travel times.

Analysis of land-use/cover, groundwater, and legacy maps

Cross tabulations between the 2006 land-use/cover maps and the legacy maps were made using the TABULATEAREA command in ArcGIS 9.1. Cross tabulations involve a cell-to-cell map comparison that generates a table with areas of similarity, i.e., cells that have the same land-use/cover class, and dissimilarity, i.e., cells that have different land-use/cover classes in each of the maps. Adding the percent area of similarity of all land-use/cover classes provides a metric of total map agreement that we call total similarity. We also examine the level of map agreement for each land-use/cover class by calculating the percent area similarity of the same land-use/cover class that matches in both maps divided by the total percent area occupied by that land-use/cover class in the current map. We refer to this calculation as land-use/cover class similarity. Finally, we also calculate a measure of land-use/cover class dissimilarity. This is accomplished by taking the percent area of land-use/cover class j, e.g., forest, that is contained in the legacy map but is currently occupied by land-cover class i, e.g., urban area, and dividing this by the total percent area that land-use/cover class i occupies in the current map.

Map overlays of surface watershed boundaries and groundwatershed boundaries were also created to compare areas of overlap between these two important hydrologic divides. We calculated the total area occupied by both hydrologic boundaries and areas unique to each. An analysis of the distribution of current land-use/covers within the 0-age groundwater travel-time maps are in what we refer to as “high-impact groundwater zones,” i.e., areas in which land use impacts surface water quickly. This region is also the most likely to be influenced by concentration inputs from overland flow, because high-impact zones are near surface water bodies.

Finally, we performed the same land-use cover/class similarity analysis described above but grouped these calculations by MCDs that fell within the groundwatershed. This provides a measure of the amount of variability in the similarity between current and legacy maps at the government level responsible for land-use planning.


Groundwater travel-time model

In Fig. 3, the surface watershed and the groundwatershed share a common area of 6230 km² (orange), whereas the areas draining to the Muskegon River solely via surface and near-surface flowpaths cover 829 km² (yellow). A 322-km² area (green) lies outside the surface watershed but is included in groundwater drainage to the Muskegon River. The most notable areas of nonoverlap of watersheds are the regions in which groundwater drains to the Muskegon River and surface water does not. These are an area in the north Black-Macatawa Watershed located south of the Muskegon River Watershed (MRW) and an area several kilometers north of Houghton Lake. The three most significant areas in which surface water drains to the Muskegon River and groundwater does not are the areas just west of Lake Muskegon, the northern portion of the watershed bordering the Manistee Watershed, and a region in the north bordering the Kawkawlin Watershed. To distinguish the areas within our analyses from the MRW proper, the area within the groundwatershed is hereafter simply referred to as the “study area.”

The spatial distribution of groundwater travel-time (GWTT) dates is highly variable (Fig. 4). Note, for example, that few GWTT dates that are pre-1990 are located in the lower portion of the watershed, i.e., around Muskegon County. In Mecosta County, GWTT dates range from the late 1990s to current along the mainstem of the Muskegon River to several large areas of pre-1900 GWTT dates located in the center of the county.

We also found (Fig. 5) that a large proportion of the cells, about 17% of the entire area, contain locations in which GWTT is less than 1 yr (shown as the number of cells for year 2006). More than one-third (34%) of the cells are located in pre-1900 GWTT areas. Another important characteristic of the distribution of these travel times is the decreasing trend in the proportion of cells with older GWTT.

Backcast land-use change model

Decadal time frames from our backcast model are shown in Fig. 6. Presettlement maps were used to represent pre-1900 land covers; 1978 and 1998 are shown here as well because they are used as inputs to four of the backcast submodels. Note that most of the presettlement map is composed of forest; approximately 4% of the map is grassland. In 1900, we estimate that a large portion of the watershed was agricultural. Agriculture is the dominant land-use/cover through to about 1960, when forests become the dominant cover. Also note the increasing amount of urban cover over the 106-yr time period. By 2006, more than 10% of the watershed was urban. Although difficult to visualize because of its scattered presence in the study area, the amount of shrub/grassland increased from the 1970s through the early 2000s, which is consistent with the patterns of land-use change over the last 20 yr.

Legacy maps

Overall, there is a 72.5% total similarity, i.e., sum of diagonals, of all land-use/cover classes between the legacy and current maps (Table 1). Note that the urban class makes up 2.11% of the study area in the legacy map and 8.73% in the current map; this discrepancy represents the largest difference in total area of any land-use/cover class between the two maps. More than 10% of the area is composed of agriculture, which is present in both maps. There is more forest (67.0%) in the legacy map than in the current map (54.2%), with slightly more than half (50.6%) of the area comprising forest present in both maps. Shrub/grassland makes up more than 2% of the cross-tabulated area, but it is more than 8% in the current map and more than 5% in the legacy map. Wetlands and barrens, both represented in small amounts in the watershed, have relatively similar amounts in both land-use/cover maps.

Note that only 22.3% (1.95%/8.73%) of current urban cells are located in the legacy map. About a third (3.74%/8.73%) of the current urban cells are forest in the legacy map, and another third of the current urban cells are agriculture in the legacy map; both are values of dissimilarity. There is a 56% similarity in agriculture between the two maps. However, there is a significant dissimilarity; nearly 8% of the area is composed of forest in the legacy map but agriculture in the current map. In fact, more than 40% (7.77%/18.7%) of current agriculture was forest in the legacy map. Few areas (2.05%) were in the transitional category of shrub in both the legacy and current maps. Two land-use/cover classes, agriculture and forest, had greater measures of dissimilarity than the measure of shrub land-use cover-class similarity (23.6%).

A zoom into one area illustrates the spatial heterogeneity of the GWTT produced by the GWTT model. Note how the southern portion of this area mostly contributes pre-1900 (dark blue) groundwater (Fig. 7). Areas bordering surface water, lakes, rivers, and streams contain short GWTT (brown). Note that most urban locations border lakes and rivers. As expected, there is also less urban area in the legacy map compared to the current land-use/cover map.

Distribution of current land uses/covers in high-impact groundwater recharge areas

Cross-tabulating land-use/cover distributions by backcast land-use change submodels shows that the most common land-use/cover and GWTT classes are agriculture (11.7%) in the pre-1900 GWTT class and forest (11.6%) in the 0-age GWTT class. As seen in Table 2, several land-use/cover-GWTT cross tabulation classes are few in number, representing less than 0.5% of the total area; these include many of the wetland GWTT classes and the agriculture class in the 1998–2005 GWTT class.

Approximately 10.7% of the total area within the 0-age GWTT map is urban; forests make up more than 44% of this area. A little more than 20% of the current agriculture is located in the 0-age GWTT class. However, comparing the percentage of a land-use/cover class within the 0-age GWTT class and the 2006 land-use/cover class underscores how disproportionately land-use/covers are located with respect to high-impact recharge areas. For example, a greater proportion of urban cells is located in the 0-age class (10.7% vs. 7.9%) than in the 2006 map because urban areas tend to develop next to water bodies; this is also true for wetlands (11.6% in 0 age class vs. 5.3% in current map). The latter class is most frequently found in the 0-age class because many of the study area’s wetlands are located along rivers and lakes.

Interestingly, only 24% of the study area’s forests are located in these high-impact groundwater recharge areas. Alternatively, more than a third (35.6%) of the current urban cells are situated in high-impact groundwater recharge areas with GWTT of less than a year. Many homes are located along the study area’s lakes and rivers.

Almost half (48.1%) of the lands in the groundwatershed currently devoted to agriculture are located in the pre-1900 GWTT class. On the other hand, 22% of current agriculture is located in the high-impact 0-age GWTT class.

Spatial distribution of similarity between current and legacy land-use/cover maps

The distribution of total similarity in these current and land-use legacy maps, which are shown in light blue, varies considerably across the study area (Fig. 8). In the minor civil division (MCD) labeled A (Mecosta Township), nearly all of the locations in the current map contain the same land-use/cover class as in the legacy map. More than half of the land-use/cover classes shown as red in Marion Township, labeled B, are dissimilar between the current and legacy maps because of long GWTT and significant transitions from presettlement vegetation into other land-use/cover classes.

We quantified current vs. legacy map variability across 53 of the study area’s MCDs but selected only those MCDs that had 95% or more of their total area located within the groundwatershed. We used the GIS to calculate total and land-use/cover class map similarities for each MCD. We then counted the number of MCDs that had total and land-use/cover class similarities binned into 10% categories and plotted these counts of MCDs vs. cumulative percent similarity. Figure 9A displays the number of MCDs vs. cumulative total similarity (block line) and land-use/cover class similarity (colored lines). These plots show, for example, that 11 MCDs had 90% or greater total map similarity for all land-use classes; 51 MCDs had 50% or greater total map similarity. Examining land-use/cover class similarity values illustrates how the maps differ. Only four MCDs had a 90% or greater similarity (Fig. 9A) for the urban class, whereas 49 MCDs had a 90% or greater similarity in the forest class. Land-use/cover class similarities varied considerably at the 50% or greater cumulative similarity category. In general, the shrub/grassland class had the least amount of similarity compared to the urban, forest, and agriculture classes, indicating that previous land-use/cover classes for most MCDs were more different in the current shrub/grassland class than in any of the other land-use/cover classes.

Figure 9B plots the cumulative percentage of dissimilarity values for land-use/cover class in those areas that were only urban in the current map. We can see from this plot that the forest class is the most common land-use/cover class in the legacy map that is now marked as urban in the current map. In fact, there are 11 MCDs that have 50% or more percent land-use/cover class disagreement that fall in the urban(current) x forest(legacy) portion of the cross tabulation table. The other three nonurban classes (agriculture, shrub, and wetland) have 20% or less land-use/cover class dissimilarity. Note that the wetland class is not present in the current maps for 14 MCDs, so this land-use/cover class does not sum to 53 for the > 0% cumulative category.


There are several issues related to model calibration and model assumptions that are important in our analysis but are beyond the scope of the current paper. For example, we are currently assessing the spatial and temporal accuracy of each land-use change submodel. The backcast model between 1939 and 1978 is being assessed (D. K. Ray and B. C. Pijanowski, unpublished manuscript) using aerial photography. Current maps are being developed using recent high-resolution aerial photography. The accuracy of land-use/cover prior to 1910 is challenging, because few records provide reliable information on the spatial distribution of land uses/covers in this watershed. However, less than 1% of the total study area is located within these 100- to 120-yr-old groundwater travel dates, although more heterogeneous models of the geologic composition of the subsurface may increase this amount.

A second important issue that needs to be addressed is how each model assumption impacts the differences between legacy and current land-use/cover maps. We made several assumptions about the land -use/cover transition pathways occurring over time, such as assuming that an agriculture → shrub → forest transition pathway was common throughout all time periods, with shrub appearing for only a 5-yr period. We based this assumption on an analysis of just two dates (1978 and 1998) of land use/cover. It is quite possible that shrub could be a managed landscape that appears for more than 5 yr in the watershed.

We also acknowledge that combining several sources of data from different spatial scales is likely to introduce errors in model output. Interpolating between data reporting periods, especially decadal census data as represented by the housing data, requires assumptions about the nature of temporal data that are likely to be oversimplified. Nevertheless, our approach does make a preliminary effort to combine data at different spatial and temporal scales that are usually available in some form in most developed countries. Further work is planned that will examine how significant these assumptions are to model output.

According to the groundwater travel-time model, travel times range from 0 to approximately 500 yr in some areas, although the vast majority of cells have travel times of less than 100 yr. Areas near streams and areas with high topographic or groundwater recharge gradients tend to have much lower travel times than areas far from connected surficial hydrologic features and located above somewhat flat regions of the water table. As expected, areas with high hydraulic conductivity and high recharge have relatively shorter travel times than those with low conductivity and low recharge. The simulation in this paper does not incorporate changes in recharge that would occur as land covers change, because the effect this has on recharge is still an open question and a subject of our continued research (e.g., Jayawickreme and Hyndman 2007, Kendall and Hyndman 2007, Hyndman et al. 2007).

We should also consider the assumptions made to develop the groundwater travel-time model. For example, the hydraulic conductivities of sediments across a watershed greatly influence travel times to surface water; unfortunately, high-resolution estimates of conductivities over such as large area are not available. Incorporating a greater degree of heterogeneity of subsurface sediments will tend to shorten travel times in some areas while increasing them significantly elsewhere. Thus, varying the conductivities, as well as recharge, over this time period could generate different travel-time map estimates. Furthermore, vegetation patterns vary over time and, as such young forests eventually mature to older forests, will alter hydrologic processes such as evapotranspiration and runoff/recharge in significant ways.

Finally, land management over time has also changed so that agriculture from, e.g., 70 yr ago has different chemical inputs than current agriculture. Thus, historical land management practices need to be studied and introduced to fully characterize the interactions of land use/cover, groundwater travel times, and current water quality measures.

There are several differences that exist between the legacy and current maps that significantly impact land-use planning and watershed management. Coupling a backcast land-use change model and a groundwater travel-time model has allowed us to generate a land-use/cover legacy map that could be used to estimate the geochemical signature of historic land use rather than relying solely on a current map of land use/cover. First, the amount of area occupied by several land-use/cover classes differs greatly between the two maps, with the most significant difference in urban areas. Across the entire watershed, less than one-third of current urban areas were present in the legacy map. Instead, those areas of the legacy map tended to contain agriculture or forests. The water quality signal between urban and forest differs greatly (Wayland et al. 2003). Interpretation of stream water-quality signals should take into account the disparity between the originating land use/cover of stream water and current land use/cover. Thus, legacy maps provide a more accurate representation of the linkage between land use and water quality. We also show that past agriculture might contribute as much to the current water-quality signal as current urban regions. Agriculture often requires nutrient and pesticide/herbicide inputs that can be carried to streams via overland or groundwater flow (Harbor 1994). Urban areas, on the other hand, may contain non-nutrient chemicals that are generated from manufacturing, transportation, and residential use activities (Fitzpatrick et al. 2007, Tang et al. 2005b).

A second important finding of our analysis is that there was a considerable amount of spatial variability in the similarity of legacy and current maps. In some cases, we found that there was a large amount of agreement, as much as 90% or greater, in land-use/cover classes between legacy and current maps. In other cases, however, there was a significant amount of dissimilarity between these two maps, including locations in which a legacy map would be far more important for land-use planning and watershed management. Our analysis of the amount of similarity of land-use/cover classes between the maps grouped by minor civil division (MCD), which is the level of government responsible for planning in Michigan, showed that many MCDs could benefit from using a legacy map as opposed to a current or recent map of land use/cover. Using a legacy map of land use/cover derived from groundwater travel time would also be more appropriate in situations in which water quality samples are taken during the summer months, when as much as 85% of the water volume at stream base flow in the Muskegon River Watershed (Hyndman et al. 2007) is being contributed by groundwater.

Another important finding from our development of a legacy map is that urban areas are found in large proportions in high-impact groundwater travel-time areas, i.e., < 1 yr. This is to be expected, because people prefer to place their homes near amenities such as lakes, rivers, and streams. Historically, manufacturing was also located near significant surface water sources such as lakes because of the need for large volumes of water in production processes. In contrast to urban areas, we found that most of the legacy agriculture was located in what might be considered low-impact groundwater areas in which simulated groundwater travel times to surface water exceed a century. These areas are important for planning and watershed management because current land-use management decisions have a very long legacy on water quality.

This special issue of Ecology and Society focused on characterizing spatial-temporal scales in ecological processes that are affected by or in turn affect society. Our legacy maps characterize land-use/cover changes in a large regional watershed over a long period from presettlement to 2006. We show how the legacy map varies in its importance across units of management such as townships, villages, and cities. We believe that land-use legacy maps provide a more accurate representation of the linkage between land use/cover and current water quality.


Responses to this article are invited. If accepted for publication, your response will be hyperlinked to the article. To submit a response, follow this link. To read responses already accepted, follow this link


The authors would like to thank the following funding sources for supporting this work: the Great Lakes Fisheries Trust, National Science Foundation WCR 0233648, EPA STAR Biological Classification, EPA STAR Multiple-Stressor, the Purdue Research Foundation, the Wege Foundation, the Great Lakes Fisheries Trust, and the NASA Land Use Cover Change Program. Amelie Davis and two anonymous reviewers provided comments on an earlier draft. The guest editor, Michael Papaik, provided very useful comments on a previous version of the manuscript that improved its readability.


Alexander, J. 2006. The Muskegon: the majesty and tragedy of Michigan’s rarest river. Michigan State University Press, East Lansing, Michigan, USA.

Boutt, D. F., D. W.Hyndman, B. C. Pijanowski, and D. T. Long. 2001. Identifying potential land use-derived solute sources to stream baseflow using groundwater models and GIS. Groundwater 39:24-34.

Brown, D. G., K. Johnson, T. Loveland, and D. Theobald. Rural land use change in the conterminous U.S., 1950-2000. Ecological Applications 15(6):1851-1863.

Changnon, S. 1992. Inadvertent weather modification in ruban areas: lessons for global cliamte change. Bulletin of American Meterological Society 73(5):619-627.

Dale, V., S. Archer, M. Chang, and D. Ojima. 2005. Ecological impacts and mitigation strategies for rural land management. Ecological Applications 15(6):1879-1892.

Farrand, W. R., and D. L. Bell. 1982. Quarternary geology of southern Michigan. University of Michigan Press, Ann Arbor, Michigan, USA.

Fitzpatrick, M. L., D. T. Long, and B. C. Pijanowski. 2007. Exploring the effects of urban and agricultural land use on surface water chemistry, across a regional watershed, using multivariate statistics. Applied Geochemistry: doi 10.1016/j.apgeochem.2007.03.047.

Foley, J., R. DeFries, G. Anser, C. Barford, G. Bonan, S. Carpenter, F. Chapin, M. Coe, G. Daily, H. Gibbs, J. Helkowski, T. Hollowy, E. Howard, C. Kucharik, C. Monfreda, J. Patz, I. Prentice, N. Ramankutty, and P. Snyder. 2005. Global consequences of land use. Science 22(309):570-574.

Freeze, R. A., and J. A. Cherry. 1979. Groundwater. Prentice-Hall, Englewood Cliffs, New Jersey, USA.

Harbaugh, A. W., E. R. Banta, M. C. Hill, and M. G. McDonald. 2000. MODFLOW-2000, the U.S. Geological Survey modular ground-water model; user guide to modularization concepts and the groundwater flow processes. Open-File Report 00-92. USGS, Reston, Virginia, USA.

Harbor, J. 1994. A practical method for estimating the impact of land use change on surface runoff, groundwater recharge and wetland hydrology. Journal of American Planning Association 60:91-104.

Hyndman, D. W., A. D. Kendall, and N. R. H. Welty. 2007. Evaluating temporal and spatial variations in recharge and streamflow using the Integrated Landscape Hydrology Model (ILHM). Michigan State University Press, Ann Arbor, Michigan, USA.

Jayawickreme, D. H., and D. W. Hyndman. 2007. Evaluating the influence of land cover on seasonal water budgets using NEXRAD and streamflow data. Water Resources Research 43: W02408, doi:10.1029/2005WR004460.

Kendall, A. D., and D. W. Hyndman. 2007. Examining watershed processes using spectral analysis of hydrologic time series. Michigan State University Press, Ann Arbor, Michigan, USA.

Maidment, D. 2002. ArcHydro: GIS for water resources. ESRI Press, Redlands, California, USA.

Niehoff, U. F., and A. Bronstert. 2002. Land-use impacts on storm-runoff generation: scenarios of land-use change and simulation of hydrological response in a meso-scale catchment in SW-Germany. Journal of Hydrology 267(1/2):80-93.

Pijanowski, B., K. Alexandridis, and D. Mueller. 2006. Modeling urbanization in two diverse regions of the world. Journal of Land Use Systems 1(2):83-108.

Pijanowski, B. C., D. G. Brown, G. Manik, and B. Shellito. 2002a. Using neural nets and GIS to forecast land use changes: a land transformation model. Computers, Environment and Urban Systems 26(6):553-575.

Pijanowski, B., S. Pithadia, K. Alexandridis, and B. Shellito. 2005. Forecasting large-scale land use change with GIS and neural networks. International Journal of Geographic Information Science. 19(2):197-215.

Pijanowski, B. C., B. Shellito, and S. Pithadia. 2002b. Using artificial neural networks, geographic information systems and remote sensing to model urban sprawl in coastal watersheds along eastern Lake Michigan. Lakes and Reservoirs 7:271-285.

Senger, R. K., and G. E. Fogg. 1990. Stream functions and equivalent freshwater heads for modeling regional flow of variable-density groundwater. 1. Review of theory and verification. Water Resources Research 26(9):2089-2096.

Tang, Z., B. A. Engel, B. C. Pijanowski, and K. J. Lim. 2005a. Forecasting land use change and its environmental impact at a watershed scale. Journal of Environmental Management 76:35-45.

Tang, Z., B. Engel, K. Lim, B. Pijanowski, and J. Harbor. 2005b. Minimizing the impact of urbanization on long-term runoff. Journal of the Water Resources Association 41(6):1347-1359.

Wayland, K., D. Long, D. Hyndman, B. Pijanowski, and S. Haack. 2002. Modeling the impact of historical land uses on surface water quality using groundwater flow and solute transport models. Lakes and Reservoirs 7:189-199.

Wayland, K., D. Long, D. Hyndman, B. Pijanowski, S. Woodhams, and S. Haack. 2003. Identifying relationships between baseflow geochemistry and land use with synoptic sampling and R-mode factor analysis. Journal of Environmental Quality 32:180-192.

Address of Correspondent:
Bryan Pijanowski
Department of Forestry and Natural Resources
Purdue University
West Lafayette, Indiana 47907 USA

Home | Archives | About | Login | Submissions | Notify | Contact | Search