Introduction

Throughout Chesapeake Bay fluvial sediments have been rapidly building deltas at the heads of subestuarine tributaries (Gottschalk, 1945; Marcus and Kearney, 1991; Pasternack and Brush, 1998). DeFries (1986) found that 75 per cent of the fine sediment carried by Port Tobacco River (located in the lower Potomac River tributary in Maryland) was deposited within the upper third of the associated subestuary. Jordan et al. (1986) reported an average 20 tons of sediment deposited per week in the upper Muddy Creek branch of Rhode River in Maryland during rainy weeks.

Stratigraphic and palaeoecologic analyses of estuarine delta deposits have been used to reconstruct the 300 year history of human-induced landscape change as well as prehistoric natural variability. In the Furnace Bay tributary of Chesapeake Bay, located on the eastern flank of the mouth of Susquehanna River at the estuary’s head, pre-European settlement sedimentation rates were < 1 mm·yr^-1. After five periods of deforestation and modern quarrying of sand and gravel, sedimentation rates peaked at 18 mm aÿ1 (Brush, 1989). Similarly, in the Patuxent River tributary on the western shore of the estuary in Maryland, prehistoric sedimentation rates were < 1 mm·yr^-1, but now exceed 4 mm·yr^-1(Khan and Brush, 1994).

While these studies provide detailed environmental histories of deltaic and estuarine sedimentation, they do not accurately estimate watershed sediment supply and transport. Reconstructing watershed sediment supply based on delta sedimentation rates is an inverse boundary-value problem: knowing the effect (sedimentation rates) one deduces the cause (sediment supply). Mathematically, an equation governing delta progradation is needed along with output from the solution domain.

Purpose

In this study we developed and presented an approach for combining palaeoecological reconstructions of sedimentation history with physically based models of delta evolution to quantify the impact of historic land-use change on sediment delivery to an estuarine delta. The approach yields simulations of the morphology and evolution of a delta that developed as a result of humaninduced high sediment supply and provides estimates of watershed sediment supply rates corresponding to observed delta sedimentation rates. Using additional field-derived relationships between elevation and plant associations, simple deterministic simulations of plant succession are also obtained.

Theory Behind Diffusion-Based Delta Modeling

To model the evolution of an estuarine delta, it is necessary to consider growth on the scale of decades at a minimum, and more likely centuries, depending on sediment supply and deposition rates. Explicit long-range modelling of processes that strongly influence sedimentation at short time scales, such as subestuarine hydrometeorology (Pasternack, 1998) and vegetation growth cycles (Pasternack and Brush, 1998), would require significant parameterization and would be computationally prohibitive. Also, sedimentary, biogeochemical and ecological processes that play a small role in delta growth over the short term may have a significant impact over the long term, and thus would be missed in a reductionistic sedimentation model.

For 30 years geomorphologists have been modelling landscape evolution with simplified forms of mass and momentum conservation (Hirano, 1968; Carson and Kirkby, 1972). The most common equation used is the diffusion equation. According to this equation, the time rate of elevation change is linearly proportional to the change in sediment transport over space, and sediment transport is linearly proportional to slope (Hanks et al., 1984; Martin and Church, 1997). Paola et al. (1992) gave a theoretical derivation of the parameters from a process-model perspective.

Diffusion has been used extensively in erosion studies of hillslope retreat (Carson and Kirkby, 1972; Saunders and Young, 1983; Dietrich et al., 1995), bluff erosion (Nash, 1980), fault scarp degradation (Colman and Watson, 1983; Hanks et al., 1984), wave-cut terrace/shoreline retreat (Hanks et al., 1984; Andrews and Bucknam, 1987; Rosenbloom and Anderson, 1994) and alluvial channel incision (Begin et al., 1981; Begin, 1988). The few studies that have used this approach in depositional settings include marine delta progradation (Kenyon and Turcotte, 1985; Rivenaes, 1992), alluvial basin filling (Syvitski, 1989; Paola et al., 1992) and foreland basin filling (Flemings and Jordan, 1989). An advantage of the approach over full processbased modelling is that there is just one parameter to calibrate: the diffusion coefficient. This parameter can characterize long-term holistic environmental interactions without requiring explicit accounting of reductionistic individual processes.

The primary problem with applying the diffusion equation to landscape evolution is the lack of evidence showing that the equation is correct for any of the above systems (Syvitski, 1989). Several investigators cite the good fit of model profiles with actual landscape profiles as evidence, but such a comparison is not a physically based test. In fact many equations, such as polynomials, exponentials or sine curves, could be made to fit landscape profiles. A true test of the validity of the diffusion equation would involve either showing that sediment transport in the setting of interest is actually linearly proportional to slope or showing that the time rate of elevation change is linearly proportional to landscape curvature. Andrews and Bucknam (1987) tested the first proposition and found that sediment transport was not linearly proportional to slope with respect to the process of shoreline scarp retreat. Instead, a cubic relation was more realistic. Heimsath et al. (1997) used data from California hillslopes to test the second proposition and found it to be valid for hillslope erosion.

Tests of the diffusion equation have not been performed for any depositional settings, but there is a strong conceptual basis for assuming the equation is appropriate there. As long as sediment transport is driven by gravity, no sediment moves when slope equals zero. As slope approaches 90 °, sediment transport increases unbounded, with significantly different trajectories for different transport relations. When slope is less than 45 °, a linear approximation is suitable regardless of the actual transport relation.

Simplified Endmembers of Delta Diffusive Growth

The duffision equation can be solved analytically for a flat-bottomed, semi-enclosed, linear, quiescent subestuary subject to any one of a number of sediment input functions. One simple solution is available for the condition when the delta head is fed by a source that maintains a constant elevation there.  The overriding problem with this simple model is the assumption of a constant delta head elevation, which asserts that most deforestation-derived sediment was transported to the head of the delta during the earliest decades of human activity. This assumption is incorrect (Brush, 1984; Ferguson, 1997). It causes a significant overestimation of delta progradation in early decades and underestimation in later periods.

Another simple solution is available for the conditions when the delta head elevation grows linearly through time. In this case, the delta progrades almost the same distance after 300 years, ending with almost the same final longitudinal profile as in the case with the assumed constant delta head elevation, but with a significantly different evolution pathway. To apply this approach it is necessary to specify (or know from site-based stratigraphic reconstruction) the forcing function governing the increase in delta head elevation per time interval. The rate of increase need not be constant.

New Methods for Calibrating the Diffusion Parameter

Methods for estimating the one diffusion-equation parameter (D) described earlier are not well suited for historic delta progradation because of the typical lack of historical surveys of delta longitudinal profile. Alternative approaches were developed for this study using additional sources of information on delta growth for the basin which provide constraining equations that close the mathematical problem. Each approach has its strengths and limitations, and when all are combined they independently constrain the value for D.

The first method is the ‘minimum D’ approach. This approach recognizes that for a given observed sedimentation rate in the first time interval of delta growth, the diffusion equation requires a balance between D and initial delta head elevation, h₀ (which controls slope). Decreasing D forces an increase in h₀, and vice versa. However, the initial delta head elevation is fundamentally limited by the known topography of the basin such that the maximum h₀ cannot be greater than the present-day watershed base elevation, which can be easily surveyed. Consequently, it is possible to establish a minimum D value with a high degree of confidence.

The second method is the ‘reservoir analogy’ approach. Many artificial reservoirs have been built on rivers that drain into estuaries. In some instances, reservoir bathymetry has been surveyed through time, yielding high-resolution records of sedimentation rates. At a minimum, many reservoirs have been surveyed once in recent years and once when built. Where such information exists, it can be used to calibrate D by assuming that the initial years of estuarine delta progradation were analogous to the initial years of reservoir delta progradation. This assumption is reasonable because watershed sediment supply rates, in both instances, far exceed rates of sediment redistribution and export within the receiving basin, and accommodation space is not limiting. Several different specific constraints are possible using this approach, depending on available reservoir data.

The third method is the ‘land-use periods’ approach. Research on the history of land use in the Chesapeake Bay basin shows incremental land clearance over time, with estuarine sedimentation reflecting land-use history (Brush, 1984). In general, 0–20 % of land was deforested for agriculture in the early to mid-18th century, 40–50 % during the mid-18th to mid-19th centuries, up to 80 per cent in the mid-19th to early 20th century, and down to 40–60 % since the early 20th century (Davis, 1985; Brush, 1986). In addition to land clearance, human activities such as mechanized agriculture, urbanization and dam building have contributed to distinct sedimentation regimes for different time periods (Ferguson, 1997). Where palaeoecologically reconstructed sedimentation chronologies show distinct delta growth regimes associated with different land-use periods, D may be calibrated by using a single shift in delta head elevation for each land-use period rather than individual shifts for each model time step. Then D is optimized to yield a simulated core chronology that is as close as possible to the overall observed chronology of the actual multiple pulses within each landuse period.

OPC Diffusion-Based Delta Landform Evolution Model

The OPC delta model is available as a simple MATLAB code that is avaiable as a text file at this link.  Save the file with the .m extension and run using MATLAB.

Inverse Delta Modeling to Obtain Forcing Function

One way to use a delta landform evolution diffusion model to simulate the historical development of a real delta is to reconstruct sedimentation rates from sediment cores collected on the delta and then adjust the model's inputs to yield the same sedimentations rates in the model as observed at time intervals at core locations. This is exactly the approach that was used in this study.

OPC Historical Maps Showing Delta Progradation

The impact of land clearance and intensive land use on OPC is evident in historic records. Maps and aerial photos from 1836 to 1994 document the rapidly prograding river-mouth delta at OPC. European settlement in the Bush River basin began in the mid-17th century, and resulted in deforestation of up to 80 % of the landscape by the beginning of this century. Prior to European settlement, native populations had no significant environmental impact on the basin (Hilgartner, 1995). In 1836, OPC was an open estuary surrounded by farmland. The period from 1850 to 1910 was marked by mining of quartz and iron ore, cannery construction, major railroad construction and a steady increase in land clearance (Gardner et al., 1988; Hilgartner, 1995). Agricultural land use peaked during 1900–1910 in Harford County (Jacobson and Coleman, 1986). Historic maps show that from 1836 to 1922 the area of the delta grew by a factor of 3
.75. Despite subsequent afforestation due to farm abandonment during 1910–1950, maps show that the delta doubled in area between 1922 and 1951. Two water supply dams were built on Winters Run by the US Army towards the end of that period. These dams impede flow, block sediment transport and attenuate storm runoff. However, the capacities of both have been reduced by > 90 % due to reservoir sedimentation (NIER, 1996). The upper dam forming Atkisson Reservoir was found to be suitable for use in the ‘reservoir analogy’ approach for estimating D. Today, the Winter’s Run basin is 48 % forest, 23 % grassland, 21 % urban/developed, 7 % farmland and 1 % other (NIER, 1996). Grassland and farms are generally located in the upper portion of the watershed, whereas the middle and lower regions are in a rapidly urbanizing corridor.

OPC Paleoenvironmental Reconstruction

Hilgartner (1995) collected multiple cores from transects at OPC to assess plant succession and habitat development in that system. A core taken at the delta head only contained a small amount of sediment before impenetrable bedload-derived gravels were reached, precluding direct determination of sediment thickness at the delta head. Among the remaining cores, one located near the longitudinal axis on the forested delta plain (Auger 2) was used for this study, as it contained the longest undisturbed environmental record for the site. Auger 2 was taken in the centre of the riparian forest at a point 968 m down the longitudinal axis from the delta head. The boring was made with a soil auger, since a vibracorer could not penetrate the root systems there. The boring produced a series of sections 16 cm long and 7 cm in diameter. Samples from each 16 cm section were bagged separately in the field. Colour, stratification, sediment texture and plant content of each sample were observed and recorded.

The dates of four stratigraphic horizons in Auger 2 were determined using radiocarbon and pollen methods (Hilgartner, 1995). The surface of the core has a date of AD 1994. The sand and gravel at the base of the boring (287 cm depth) contained woody material and seeds that were assessed at Beta Analytic, Florida, and found to have a ¹⁴C date of 680 years BP. The oak to ragweed pollen ratio is known to be >10
.0 prior to the date of European settlement in an area. Based on the well-documented historical record from upper Chesapeake Bay, a date of AD 1730 was assigned to the depth where this ratio decreases to between 1
.0 and 5
.0, marking colonial land clearance of 5–20 %. A date of AD 1780 was assigned when the oak to ragweed ratio falls below 1
.0, marking the period of intensive agriculture when land clearance reached 40–50 %. In Auger 2 both of these pollen transitions were very abrupt. They occurred at 248 and 224 cm, respectively.

Palaeoecological reconstruction of the history of delta evolution at OPC shows that all of OPC was subtidal prehistorically, with lush gardens of submerged aquatic vegetation and sedimentation rates of < 1 mm·yr^-1. For c. 400 years prior to 1730 the stratigraphy of Auger 2 shows a mix of gravel, sand and silt. The macrofossils were from Zannichellia palustris, Najas sp., Elodea canadensis and Potamogeton diversifolius – all submerged aquatic plant species. The bed elevation, sediment stratigraphy and macrofossils from that period are all consistent with subtidal conditions. Beginning c. 1730, sediment stratigraphy changed to clay and silt, and sedimentation rates increased sharply by a factor of six from 0
8 to 5 mm·yr^-1. The maximum sedimentation rate of 34
6 mm yr^-1 occurred during the 20 year interval 1840–1860. After 1920 sedimentation decreased to nearly prehistoric levels, most likely due to the dams on Winters Run. The lack of any sand or gravel lenses interrupting the 270 years of sedimentation of fines since 1730 indicates that channel processes never interfered with the depositional record at the site of Auger 2. Based on changes in plant macrofossil assemblages over time, Hilgartner (1995) concluded that the vicinity of Auger 2 went through a succession of habitats from subtidal front in prehistoric times to middle marsh, shrub marsh and finally riparian forest. Taken together, the maps and sedimentary record reveal a rapidly prograding river-mouth delta with a succession of plant associations.

OPC Diffusion Constant Parameterization

Values for D were obtained using the three new methods described above. Using the ‘minimum D’ approach constrained by a maximum watershed base level of 6
1 m and an initial 20-year average sedimentation rate of 5 mm aÿ1, a D of 3,763 m²·yr^-1 was obtained. To use the ‘reservoir analogy’ approach, information on sedimentation in Atkisson Reservoir was obtained from the US Army. The average depth of the reservoir decreased by 3
3 m over the initial 54 years of its existence (1942 to 1996), yielding a long-term average sedimentation rate of 6
11 cm·yr^-1 (NIER, 1996). When this same long-term average sedimentation rate was imposed on the initial 54 years of OPC sedimentation, a D of 6,199 m²·yr^-1 was calculated. Finally, to use the ‘land-use periods’ approach, D was solved for using a scheme that optimized the coefficient of determination for the linear relationship between real deposition and model deposition in a run in which only three Dh0 values were used (Figure 10). The diffusivity (D) was found to be 5,709 m²·yr^-1.

To assess the viability of these D values prior to detailed model runs, the evolution of delta head elevation through time was assessed. Under the ‘minimum D’ approach 7
93 m of aggradation had to occur at the delta head in the initial 20 years (39.
65 cm·yr^-1) to achieve the observed 5 mm·yr^-1 at the location of Auger 2. However, the sedimentation rate at Auger 2 for the second time interval is the same as for the first, forcing the model to erode 5
5 m of sediment at the head. This degree of aggradation and degradation is not physically realistic. While the ‘minimum D’ approach provides an absolute constraint on the minimum rate of sediment spreading, it is not a reasonable estimate for the actual spreading rate. In contrast, both the ‘reservoir analogy’ and ‘land use periods’ approaches yield smaller changes in delta head elevation over time.

Results of Best OPC Delta Model

The evolution of the OPC delta was simulated using the diffusion equation model the ‘reservoir analogy’ D value to match the sedimentaiton rates at Auguer 2 using thirteen 20-year average delta head elevation change values. Four distinct periods (a 1740–1760 initial pulse, a 1760–1780 erosive/redistributive interval, a 1780–1920 growth period, and a 1920–present erosive/redistributive era) are evident in both simulations. The initial pulse is so large in part because it has to account for the difference between the actual initial bed elevation longitudinal profile and the assumed flat-bottomed profile. Given that the sedimentation rate in Auger 2 did not change between 1740 and 1780, the models show that significant erosion and export of material out of the system must have occurred to offset the spread of the initial pulse. Again, this may be influenced by the choice of initial conditions. After 1780 the model settles down and shows persistent growth until 1920. Atkisson and Van Biber dams were built during the 1920–1940 interval. Since then the sediment influx to OPC has been small enough that delta progradation has been dominated by redistribution of sediment from the head of the delta to the delta front. Field observations indicate that the active mechanism for this redistribution is rapid migration and entrenchment of previously ephemeral channels in the riparian forest brought on by sediment-starved flow.

OPC Historic Sediment Suppy Rates

Based on the simulated growth of the OPC delta it was possible to calculate the time series of 20-year average supply rate. When the supply rate is multiplied by corresponding sediment bulk densities reported for Auger 2 and divided by watershed area, minimum estimates for 20-year average watershed erosion are obtained. These values are minima because most of the coarse fraction of eroded sediment was never delivered to the delta and an unknown quantity of wash load passed over the delta into Chesapeake Bay. Peak sediment supply and watershed erosion occurred during 1840–1860, resulting in an erosion rate one order of magnitude greater than the pre-settlement rate. To estimate the suspended sediment concentration corresponding to sediment supply, the latter was multiplied by the long-term average Winters Run discharge observed at the USGS gauging station. While this approach is crude and does not account for climatological variations in 20-year average flow rates, the estimated suspended sediment concentration for 1980–2000 is within the range of values reported for Winters Run (NIER, 1996). Based on information from that study, the peak 20-year average concentration (1840–1860) is equivalent to the suspended sediment concentration during a present-day flood of c. 10
.2 m³·yr^-1, which has a 1.
25 year recurrence interval. This finding suggests that on any given day sediment loads prior to the Civil War (1861– 1865) were as high as those observed in higher magnitude discharges today.

Simulated Plant Association Succession

Historic succession of plant associations at OPC was extrapolated from model results based on an empirically derived relationship between elevation and plant species distributions. According to the empirical relationship, elevation accounts for 85 % of the spatial variability in habitat (Pasternack, 1998). From 1740 to 1840 no riparian forest was present. Driven by extreme sedimentation rates, the biggest change in habitats occurred in the interval 1840–1860 when the forest became established. It is possible that the additional surface roughness of the forest caused a positive feedback that further accelerated the rate of sediment accretion. After that time, the forest grew rapidly at the expense of the subtidal front. Meanwhile, intertidal habitats maintained a nearly constant proportionality of 1
8 units of low marsh, 1
6 units of middle marsh, 1
5 units of high marsh and 1
4 units of shrub marsh per unit of floating leaf habitat after 1840.

Model accuracy can be tested by comparing the predicted 1980–2000 relative habitat areas with those mapped in the field during 1990–1997. Actual relative habitat areas were calculated using Arc View Geographical Information System and an OPC habitat map. The habitat map was based on an April 1994 colour infrared digital orthophoto (1
2 m resolution) obtained from the Maryland Department of Natural Resources and vegetation data collected during 1990–1997 using standard ecological methods (Hilgartner, 1995; Pasternack, 1998). At the scale of the three major delta zones, subtidal, intertidal and forest, model predictions are within 14, 4 and 9 per cent respectively of the actual relative areas. While this result is very strong, the distribution of habitats within the intertidal zone shows significant deviation, with error ranging from a minimum of 20 per cent for the high marsh to a maximum of 76 per cent for the low marsh (Figure 15). The origin of this error stems from the narrow range in elevation spanned by plant habitats in the intertidal zone, and reflects the limits of the empirical equation as much as those of the physical model.

Discussion

The diffusion equation-based inverse boundary-value technique applied in this study yields a quantitative assessment of the impact of historic land-use change on sediment delivery to estuarine deltas. At the basin scale, the technique produces conservative data on sediment supply rates. Changes in deduced supply rates indicate real changes in basin functioning. The two possible causes of such changes are hydroclimatological variability and land-use impacts. For centuries before significant land-use impacts occurred in the Winters Run basin, hydroclimatology varied over a wide range including large magnitude events such as hurricanes and probably fires. However, all of the cores analysed by Hilgartner (1995) show no evidence of any impact of that variability in Otter Point Creek in sedimentation rates, pollen distributions or macrofossil profiles over a 1000 year period prior to settlement. At a minimum this indicates that the geomorphic recovery time of the subestuary to natural disturbance was shorter than the recurrence interval of any disturbance. Beyond that, it is possible that hydroclimatology did not induce any significant basin-scale impacts prior to European settlement. Regardless of which of these two theories is correct, it may be concluded that Otter Point Creek was a sustainable submerged macrophyte habitat under the natural conditions existing prior to European colonization.

Given the lack of any recorded geomorphic impact of hydroclimatology in the natural state of this system, land-use impacts must be the direct cause of recorded changes. Up to 80 per cent of the basin was deforested, and peak receiving-basin changes occurred during the period of most intense basin land use. Also, sedimentation rates in the subestuary dropped directly in response to reservoir construction which blocked sediment delivery. While Khan and Brush (1994) reported large impacts of hurricanes on the Patuxent watershed–subestuary system this century, the lack of pre-settlement impacts points to human land use precipitating these impacts by disturbing basins beyond critical recovery thresholds. Consequently, postsettlement geomorphic disturbance observed in Chesapeake Bay tributaries in apparent response to hydroclimatological events should be recognized as the direct consequence of human land use.

Beyond elucidating the role of human land use versus hydroclimatology in watershed disturbance, this study provides insight on the long-term dynamics of Chesapeake Bay subestuarine deltas. Calibrated values for the diffusion constant show that Chesapeake Bay sediment is dispersing downslope at about the same rate (10⁴ –10⁵ m²·yr^-1) as reported for other depositional systems, according to data summarized by Flemmings and Jordan (1989). Because delta sediment is moving downslope relatively quickly, the post-reservoir reduction in sediment supply will not protect Bush River or Chesapeake Bay from continued sediment impacts, which will now originate in the delta itself. Furthermore, vegetation and its roots may be reducing the rate of sediment dispersion, but the latter is still high enough to promote continued delta progradation.

From an applied standpoint, reconstructed sediment supply rates could be used within the US Environmental Protection Agency’s regulatory framework to establish reasonable estimates of sediment ‘total maximum daily loads’. Sediment is a known contaminant of physical habitats. At present, attempts are being made to promulgate regulations that limit sediment loads moving down many human-impacted rivers. Estimates of historic loadings based on the method used in this study show the range of natural and disturbed sediment loadings in a basin. The pre-settlement natural loading provides an objective basis for establishing non-impact loadings, if such a strict criterion is desired. Less strict standards could be based on observed intermediate loadings.

Finally, it is important to stress and further consider fundamental theoretical limitations of the inverse approach. The accuracy of the delta progradation simulations is dependent on the precision and resolution of receiving basin sedimentation chronologies. The pollen profile analysis methodology for reconstructing sedimentation rates has proven accurate in many past studies, and the small error in prediction of actual OPC delta zones independently corroborates the reconstructed history of Auger 2. However, it is possible that there are significant higher frequency fluctuations in sedimentation rate that were not recoverable using the coring technology at hand. Infinite inverse boundary value solutions exist at any frequency higher than the sampling frequency achievable with the Auger 2 chronology. At some threshold high frequency, the diffusion equation would no longer be a viable means for modelling delta progradation, so improved sedimentation reconstruction methods would only be useful to a point. At this time there is no theoretical basis for pinpointing a sedimentation or erosion interval that is optimal for modelling the diffusion of landscapes, but this may need to be investigated further if higher resolution assessments of basin scale dynamics are to be attempted.

Publications

• Hilgartner WB. 1995. Habitat development in a freshwater tidal wetland: a palaeoecological study of human and natural influences. PhD dissertation, Johns Hopkins University, Baltimore.

  Pasternack, G. B. 1998. Physical dynamics of tidal freshwater delta evolution. Ph. D. Dissertation. The Johns Hopkins University, 227pp, 5 appendices.

• Pasternack, G. B., Brush, G. S., and Hilgartner, W. B. 2001. Impact of Historic Land-Use Change on Sediment Delivery to an Estuarine Delta. Earth Surface Processes and Landforms 26:409-427.