Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Geoderma. Author manuscript; available in PMC 2010 October 15.
Published in final edited form as:
PMCID: PMC2951733

Geostatistical modeling of the spatial distribution of sediment oxygen demand within a Coastal Plain blackwater watershed


Blackwater streams are found throughout the Coastal Plain of the southeastern United States and are characterized by a series of instream floodplain swamps that play a critical role in determining the water quality of these systems. Within the state of Georgia, many of these streams are listed in violation of the state’s dissolved oxygen (DO) standard. Previous work has shown that sediment oxygen demand (SOD) is elevated in instream floodplain swamps and due to these areas of intense oxygen demand, these locations play a major role in determining the oxygen balance of the watershed as a whole. This work also showed SOD rates to be positively correlated with the concentration of total organic carbon. This study builds on previous work by using geostatistics and Sequential Gaussian Simulation to investigate the patchiness and distribution of total organic carbon (TOC) at the reach scale. This was achieved by interpolating TOC observations and simulated SOD rates based on a linear regression. Additionally, this study identifies areas within the stream system prone to high SOD at representative 3rd and 5th order locations. Results show that SOD was spatially correlated with the differences in distribution of TOC at both locations and that these differences in distribution are likely a result of the differing hydrologic regime and watershed position. Mapping of floodplain soils at the watershed scale shows that areas of organic sediment are widespread and become more prevalent in higher order streams. DO dynamics within blackwater systems are a complicated mix of natural and anthropogenic influences, but this paper illustrates the importance of instream swamps in enhancing SOD at the watershed scale. Moreover, our study illustrates the influence of instream swamps on oxygen demand while providing support that many of these systems are naturally low in DO.

Keywords: Sediment oxygen demand, Geostatistics, Sequential Gaussian Simulation, Blackwater streams, Dissolved oxygen, Organic Carbon

1. Introduction

The study of riverine ecosystems has been furthered by extending a principle of landscape ecology to aquatic ecosystems—namely that streams should be viewed as a diverse assemblage of multiple habitat patches that act as a mosaic across the landscape (Pringle et al., 1988). The River Continuum Concept was one of the first attempts to differentiate waterways into more discrete patches and compare different lotic environments based on location within a longitudinal gradient and across stream order (Vannote et al., 1980). Aquatic ecologists have further subdivided these reaches into relatively homogenous units such as riffles, pools, and runs for ease of study and comparison (Pringle et al., 1988). These designations have proved useful in characterizing different aspects of stream function and processes that occur within them. The ability to quantify and map these patches or the features within them on a reach scale could be of particular benefit.

The changing viewpoint of a stream as an integral, variable force operating as a connected entity to the surrounding catchment has modified the ways streams are studied and perceived (Bencala, 1993). Today, a stream is less often viewed primarily as a pipe transporting water, solutes and materials downstream, but rather as a dynamic system, involving bidirectional connections between the stream and landscape as a whole (Bencala, 1993). One such system where this mosaic of patches and dynamic bidirectional connections are readily evident is within blackwater river watersheds of the southern Atlantic Coastal Plain. Within Georgia, blackwater rivers have not been impacted by large impoundments and in many areas maintain intact riparian buffers over the majority of their length (Smock and Gilinsky, 1992; Katz and Raabe, 2005). These streams are characterized by a marked seasonal change in discharge and areal extent that are driven largely by precipitation events which cause extensive floodplain inundation lasting from winter to early spring, followed by drying (sometimes completely) when evapotranspiration becomes dominant (Wharton, 1978; Wharton et al., 1982). Streams up to at least 5th order flow intermittently.

Because of this pronounced change in discharge and areal extent, a critical and abundant feature of blackwater streams are large instream swamps which form a vital link between the aquatic and terrestrial environment. Because of a lack of topographical relief, floodplains can become inundated for hundreds of meters in width during periods of high discharge. This marked change in inundation between periods of high and low discharge can lead to drastic differences in wetted perimeter, width/depth ratio and hydraulic roughness (Hupp, 2000). Blackwater watersheds are characterized by low gradients, but the slight changes in elevation that are present cause a variety of hydrologic, soil and vegetational patterns (Conner and Buford, 1998; Burke et al., 2003).

Many blackwater streams are also characterized by having extremely low dissolved oxygen (DO) levels, one of the many measures of water quality subject to regulation under the Clean Water Act. Considered to be an excellent indicator of stream biological activity and health and among the most important methods when investigating instream habitats (Wetzel and Likens, 2000), many of the streams draining the Georgia Coastal Plain are listed in violation of the state’s DO criteria. Within the state of Georgia, the standard for DO is defined as: “A daily average of 5 mg L−1 and no less than 4 mg L−1 at all times for waters supporting warm water species of fish” (GADNR, 2005). Of the river segments listed as impaired within the four main blackwater river basins in 2000–2001, 90% of the listings were a result of failure to meet the DO standard (GADNR, 2000–2001). In the 2004 303(d) list, the number of segments still listed as impaired for the DO standard made up 82% of all listed segments.

Despite the high prevalence of impaired streams due to low DO, evidence suggests that low DO may be a natural phenomenon in these systems and not a sign of pollution or impairment. Levels of DO below the 5 mg L−1 limit have been shown in several studies, even those in largely forested watersheds or where there is extensive riparian vegetation (Joyce et al., 1985; Benke and Meyer, 1988; Ice and Sugden, 2003). One of the hypothesized reasons for relatively low DO is the extended contact time with organic sediments.

Sediment oxygen demand (SOD) is defined as the rate which DO is consumed from the water column due to the decomposition of organic matter in bottom sediments (Hatcher, 1986). Organic materials can come from either outside the system (allochthonous material such as leaf litter and wastewater point discharges) or from within the system (autochthonous material such as decomposing plant materials). Once organic matter reaches the sediment matrix, SOD is influenced by two different phenomena: (1) the rate at which oxygen diffuses into the sediments and is then consumed; and (2) the rate at which reduced organic substances are conveyed into the water column and then oxidized (Bowie et al., 1985). The literature is consistent in describing SOD as the combination of two processes: (1) Biological respiration of benthic organisms residing in the sediment and (2) chemical oxidation of reduced substances found within the sediment matrix (Bowman and Delfino, 1980; Hatcher, 1986; Chau, 2002).

SOD has been shown to play a major role in instream water quality, as it can be a critical sink of DO within a river system (Wu, 1990). In some rivers, SOD has been shown to be over half of the total oxygen demand (Rutherford et al., 1991; Matlock et al., 2003). Measures of SOD within blackwater streams, and especially within instream swamps, are often higher than estimated values for southeastern sandy bottomed rivers and streams (Cathey, 2005; Crompton, 2005; Utley et al., 2008; Todd et al., 2009). Despite this major influence on the total oxygen demand within a system, SOD is often assumed (or estimated) in modeling exercises since measuring this parameter in the field is problematic due to limited access in many areas, the need for specialized equipment and training, and its time consuming nature (Hatcher, 1986; Matlock et al., 2003). Together, estimating or assuming SOD rates in models can lead to inaccurate water quality models with a potential of great biological and financial cost.

Along with the challenges listed above, in situ measurements of SOD are localized in nature. Potentially, the sediment physical and chemical properties that influence this rate change markedly across a stream bottom, thereby making extrapolation and characterization of SOD’s true impact on the reach or watershed scale problematic (Bowie et al., 1985). In an effort to relate SOD to a more easily measured field parameter, Todd et al. (2009) found that SOD was most significantly correlated with total organic carbon (TOC) in the top 5 cm of the sediment (Fig. 1, r2=0.358; p=0.002). Using this relationship between SOD rate and TOC within the top 5 cm, we investigate two sites subject to different flooding regimes and predict the distribution of SOD values at a reach scale based on geostatistical analysis. Using GIS analysis, we also map and characterize the importance of these areas at a watershed scale. In total, we utilize these findings to identify areas of possible intense oxygen demand within a reach and address the importance of instream swamps in contributing to high DO demand on a watershed scale.

Fig. 1
Relationship between total organic carbon (TOC) and sediment oxygen demand (SOD). From Todd et al. (2009).

2. Methods

2.1. Study Site

Research was conducted in part of the Little River Experimental Watershed (LREW), a 334 km2 research watershed located in the western headwaters of the Suwannee River Basin (Fig. 2). The LREW, a characteristic blackwater river system, was instrumented by the Southeast Watershed Research Laboratory of the USDA Agricultural Research Service for the measurement of rainfall and streamflow beginning in 1967 and has been designated as representative of the soils, topography, geography and land use within the southern Coastal Plain (Sheridan and Ferreira, 1992). Land usage across the watershed is 50% forested, 31% agricultural, 10% pasture, 7% urban, and 2% open water (Bosch et al., 2006). While land use is primarily agriculture and forestry, riparian vegetation remains largely intact along portions of the river, with swamp hardwood communities consisting of a closed canopy and thick undergrowth. The most common tree species found in these systems are red maple (Acer rubrum var. trilobum), Ogeechee tupelo (Nyssa ogeche), water oak (Quercus nigra), swamp black gum (N. sylvatica var. biflora), and bald cypress (Taxodium distichum).

Fig. 2
Map of Little River Ecological Watershed with selected sampling sites.

The area is characterized by low gradients of less than five percent and extensive inundation of wooded floodplains and instream swamps when streamflow is moderate to high (Cathey, 2005). Designation of stream channels may be difficult, with identification only possible during low flow conditions when streamflow is confined to braided channels.

The areas selected for measurement are in two locations. The first is a 270 m long stretch of stream located in the headwaters of the LREW above the gauging site designated as Station J (Figs. 2 and and3)3) (Sheridan and Ferreira, 1992). The stream at this point is 3rd order, with a bankful main channel 5–15 m in width and subject to overbank flooding during high flow events where floodplain width can extend to over 100 m. The second location is a 1600 m long stretch of river located in the lower part of the LREW above the gauging site designated Station B (Figs. 2 and and4).4). The stream at this point is 5th order and can be as wide as 350 m during periods of complete inundation. At this site, inundation of the floodplain usually begins in December with complete inundation until April or May. During the summer months, flow may stop along with complete drying of the river channel at both locations.

Fig. 3
Map of 3rd order sampling site.
Fig. 4
Map of 5th order sampling site. The white areas represent uplands that are only flooded during very high flow events.

2.2. Soil Collection and Preparation

Soil samples were taken at the 5th order site in summer 2006 and the 3rd order site in fall 2006 when the channel was completely dry. Samples at the 5th order site were taken every 25 m along five transects spaced roughly 150 m apart (Fig. 4). The 3rd order site had a total of 9 transects spaced roughly 35 m apart (Fig. 3). Due to the narrower area of inundation at this location and to capture areas that were flooded most often, a sample was taken at the center of the channel, 5 m on either side of center and then 10 m apart thereafter. Samples were collected at a total of 80 points for the 5th order site and 68 points were collected for the 3rd order site. One soil core (5.7 cm diam.) was taken from each quadrant surrounding a sampling point. Soil samples were taken to a depth of 15 cm with each set of four cores separated into two depth classes (0–5 cm and 5–15 cm) in the field. All four cores from each depth class were composited into a single Ziploc bag. Upon return to the lab, soil samples were stored in a refrigerator until analyzed. They were dried in an oven at 35 °C and were then crushed and sieved with a #10 sieve (<2 mm). Material that passed the sieve was considered soil matter with that remaining above as litter/residue.

2.3. Organic Carbon Determination

From sieved samples, approximately 20 g of soil was pulverized using a rolling table until the sample was completely crushed. Three 3 g subsamples were taken from the pulverized sample and organic matter was determined by loss on ignition (LOI) (Nelson and Sommers, 1996). After soil organic matter determination by LOI, the organic matter content of each sample was converted to soil organic carbon (SOC) by assuming a 0.5 proportion of C by weight (Nelson and Sommers, 1996). Carbon content in litter residues (LOC) was assumed to be 0.4 by weight (Schlesinger, 1997). The litter and soil organic carbon pools were added together to estimate total organic carbon (TOC). Established conversion factors (Nelson and Sommers, 1996; Schlesinger, 1997) were used to arrive at final measures of SOC and LOC due to equipment limitations that prevented the experimental measure of organic carbon directly.

2.4. Geostatistical Analysis

Geostatistics provides a set of statistical tools that incorporate the spatial coordinates of observations into the description and modeling of spatial patterns, predictions at unknown locations, and assessment of the uncertainty associated with those predictions. This section highlights the sequence of steps adopted for the present analysis; a more complete discussion of each geostatistical procedure is given by Goovaerts (1997; 2001).

First, Sequential Gaussian Simulation (SGS) was used to generate realizations of the spatial distribution of TOC at both experimental sites. Each realization exactly matches the sampled data (conditional simulation) while reproducing their spatial pattern as modeled from the semivariogram that measures the average dissimilarity between data as a function of their separation distance (Goovaerts, 2001; Vann et al., 2002). In total, 10 simulated maps of TOC were generated at both sites, using a spacing of 2 m at the 5th order site and 1 m at the 3rd order site. Cross-validation was used to assess the prediction accuracy and the quality of the model of local uncertainty.

The second step was to transform each TOC map into a SOD map, accounting for the uncertainty attached to the coefficients of the linear regression model. Let y(u) and z(u) be the SOD and TOC values at a given location u, respectively. Using the significant linear relationship (Fig. 1, p =0.002) developed between these two attributes, the expected SOD value at any unsampled location ui can be predicted as y(ui)= f(z(ui))=a + bz(ui). The uncertainty about the predicted y-value arises from the uncertainty attached to: 1) the z-value (i.e. TOC), and 2) the regression coefficients a and b. The uncertainty associated with the TOC value was modeled through the generation of 10 TOC maps, {z(l)(ui), i =1,…, N; l=1,…, L′}, with N the number of simulated nodes and L′ the number of realizations, using Sequential Gaussian simulation. The parameter uncertainty can be modeled by building the probability distributions of each parameter a and b as well as their joint distribution. For linear regression, these distributions are Gaussian and fully characterized by the parameter estimates and standard errors, as well as the covariance between parameters. Both sources of uncertainty can then be incorporated numerically by sampling randomly the two parameter distributions and combining the simulated coefficient values {(a(l), b(l)), l=1,…, L} with the simulated TOC values in the regression equation to retrieve the corresponding simulated SOD values y(l)(ui) = f(z(l)(ui))=a(l)+b(l)z(l)(ui). In this paper 1,600 pairs of regression coefficients were generated using the method described by Rossel et al (2001), yielding a total of 16,000 SOD maps after combination with the 10 simulated TOC maps.

The third step was to average the TOC and SOD simulated values across all N nodes (area-wide averages) and across all L′ realizations (E-type estimate). The distribution of 16,000 area-wide averages provides a numerical model of the uncertainty attached to the mean SOD value over the site A. Extreme scenarios, such as the map of simulated SOD values out of the 16,000 possible realizations that yields the lowest or the highest area-wide average can be identified easily. These two realizations will be referred to as the minimum and maximum area-wide averages respectively. The E-type estimate was calculated as the mean of 10 simulated TOC maps or 16,000 simulated SOD maps.

All geostatistical analyses were conducted using Stanford Geostatistical Modeling Software (SGeMS) (Remy, 2004) and the variograms modeled using TerraSeer Space-time Information System (STIS). All statistical computations were performed using SAS for Windows (Version 9.1, SAS Institute Inc., Cary, North Carolina, 2002–2003).

2.5. Floodplain Soils Calculation

Floodplain soils are a common feature of the LREW and the experimental locations described above are only an example of the extent of instream swamps within the greater LREW. An estimate of areas subject to elevated SOD was calculated using existing geographical data layers. The extent of the basin was delineated as described by Sheridan and Ferreira (1992). The soils data for the watershed were retrieved from the USDA-NRCS Soil Survey Geographic (SSURGO) database with all datasets at the 1:24,000 scale, but publication dates range from 1999–2004 for individual counties (USDA-NRCS, 2008). After combination of datasets, hydric soils meeting either of two criteria were extracted: 1) Those soils having a landform designation of “Floodplains” and/or 2) Those soils having a hydric criterion of 3 or 4. A designation of 3 is defined as “soils that are frequently ponded for long or very long duration during the growing season” and a designation of 4 is defined as “soils that are frequently flooded for long or very long duration during the growing season” (USDA-NRCS, 2008). A hydrological layer with all streams within the LREW from the USGS National Hydrography Dataset (NHD) at the 1:24,000 scale was used and individual stream segments delineated by the Strahler stream order classification method. Each stream segment was buffered to encompass all extracted hydric soil along that segment and area of floodplain soil along each segment calculated and summed by stream order. For the purposes of this paper, all streams 1st through 3rd order were considered representative of the 3rd order site and headwater streams. Streams 4th and 5th order were considered representative of the 5th order site. Additionally, the gradient was computed in each watershed by downloading digital elevation models at the 1:24,000 scale from the USGS National Elevation Dataset (NED).

All geographical analysis was completed using ESRI ArcMap (Version 9.2, ESRI Inc., Redlands, California, 1999–2006) and ArcView (Version 3.3, ESRI Inc., Redlands, California, 1991–2000).

3. Results

3.1. Geostatistical Analysis

Each site showed a distinct spatial correlation as illustrated by differences between the TOC variograms (Fig. 5a and b). At the 3rd order stream site, the experimental variogram, created using 10 lags of 14 m, was fitted by an exponential model with a nugget of 0.41, a range of 27.3 meters, and a total sill of 1.02 (Fig. 5a). The nugget effect is much smaller for the variogram computed at the 5th order stream site (9.92*10−8) which was created using a series of 20 lags of 20 m (Fig. 5b). The model is also exponential with a sill of 1.03 and a range of 64.46 m.

Fig. 5
Experimental variograms (dotted line) and variogram models (full line) of the total organic carbon (TOC) normal scores used for Sequential Gaussian Simulation for analysis of 3rd order site (a) and 5th order site (b).

The E-type estimates of TOC at the 3rd order stream site ranged between 22.68 and 225.91 mg gtotal1 with a mean of 73.47 mg gtotal1 (Table 1). In contrast, the E-type TOC estimates at the 5th order stream site were on average nearly double that of the 3rd order stream site averaging 145.61 mg gtotal1 and ranging between 15.11 and 275.28 mg gtotal1. For SOD, the E-type estimates had a mean of 4.90 g O2 m−2d−1 at the 3rd order stream site(Table 1; Fig. 6a). The area-wide minimum and maximum realizations had average SOD rates of 2.32 g O2 m−2d−1 and 7.49 g O2 m−2d−1 respectively (Table 1; Fig. 6b and c). The average of all E-type SOD estimates at the 5th order site was higher than the 3rd order site with a rate of 7.13 g O2 m−2d−1 (Table 1; Fig. 7a). The area-wide minimum and maximum averages of SOD were also higher than their 3rd order counterparts showing rates of 4.37 g O2 m−2d−1 and 11.13 g O2 m−2d−1 respectively (Table 1; Fig. 7b and c). The realizations also show a difference in distribution patterns and patchiness between the two locations. The higher rates of SOD tend to occur away from the center of flow (dashed line) in the 3rd order site while at the 5th order site the higher rates are generally found in areas following the center of flow (Figs. 6a and and7a7a).

Fig. 6
E-type estimate (a), area-wide minimum (b), and area-wide maximum (c) realizations of sediment oxygen demand (SOD) (g O2 m−2d−1) at 3rd order site. SOD classes constitute equally populated intervals of total data set at 3rd order location. ...
Fig. 7
E-type estimate (a), area-wide minimum (b), and area-wide maximum (c) realizations of Sediment oxygen demand (SOD) (g O2 m−2d−1) at 5th order site. SOD classes constitute equally populated intervals of total data set at 5th order location. ...
Table 1
Summary statistics (mean and range) of total organic carbon (TOC) and sediment oxygen demand (SOD) rate computed for: 1) the averaged realization (E-type estimate calculated as the mean of 10 simulated TOC maps or 16,000 simulated SOD maps), and 2) the ...

The accuracy of the geostatistical model was assessed using the “leave-one-out” or cross-validation approach, whereby one observation is removed at a time and estimated from the remaining data. This validation was conducted in an estimation framework (i.e. kriging) since it would be too time-consuming to perform in a simulation framework. To be consistent with the stochastic simulation approach, kriging was conducted on the normal scores and results were back-transformed using the procedure described in Saito and Goovaerts (2000). The scatterplots of estimated versus observed TOC values indicate that the prediction is unbiased for both sites (Fig. 8, top graphs). However, the large nugget effect of the variogram for the 3rd order site causes a strong smoothing effect and a weak correlation between predicted and sampled values.

Fig. 8
Results of the cross-validation approach for the two sites: scatterplot of estimated versus observed total organic carbon (TOC) values, and accuracy plots that show the fraction of true TOC values falling within probability intervals of increasing size. ...

The quality of the model of local uncertainty was explored using the accuracy plot (Goovaerts, 2001). At each location, we computed a series of symmetric p-probability intervals (PI) bounded by the (1−p)/2 and (1+p)/2 quantiles of the normal distribution defined by the normal score kriging estimate and variance. For example, the 0.5-PI is bounded by the lower and upper quartiles. A correct modeling of local uncertainty would entail that, for example, there is a 0.5-probability that the actual TOC value falls into the 0.5-PI or, equivalently, that over the study area 50% of the 0.5-PI include the true TOC value. The “accuracy plot” in Fig. 8 (bottom graphs) allows one to visualize departures between observed and expected fractions as a function of the probability p. These plots indicate an accurate model of uncertainty where the expected fractions are generally exceeded, in particular for the second site. This good agreement is confirmed by the large value of the goodness statistic (Goovaerts, 2001).

3.2. Floodplain Soils Extent

The LREW is characterized by an extensive floodplain network along all stream orders (Fig. 9). The three soil types meeting the hydric soil criteria defined above were Grady sandy loam (Gr), Kinston and Osier soils (KO), and Olustee sand (Os) with the majority being KO soils. There are a total of 553.46 km of streams within the LREW with 89% of that length classified as headwater streams. Within the entire watershed, nearly 2,700 ha are identified as floodplain soils with increasing area of floodplain per stream km with increasing stream order (Table 2).

Fig. 9
Map of floodplain soils within the Little River Ecological Watershed.
Table 2
Total length of designate stream orders, floodplain area, and floodplain area per length of stream by stream order in the Little River Experimental Watershed.

4. Discussion

Due to the large influence of SOD on watershed scale DO models and overall water quality, especially within the LREW (Cathey, 2005), analyzing its effect across a landscape is of vital importance. Nevertheless, extrapolating point in situ measurements of SOD across a reach or watershed are problematic. The results of this study show that relating SOD to a less time intensive and more easily measured variable (TOC), along with geostatistical analysis, allows for not only the mapping of SOD across a landscape but also the ability to quantify the importance of instream swamps.

The mapping of TOC (and by relation SOD rates) at our two sites shows distinct differences. While concentrations of TOC are high at both sites, the 5th order site has a TOC concentration that is nearly double that of the headwater, 3rd order site (Table 1). Highly organic soils are not an uncommon feature of many blackwater Coastal Plain environments. For example, Burke et al. (2003) reported floodplain organic matter levels on a South Carolina blackwater river ranging from 1.1 to 8.1% depending on the flooding and vegetation regime. Higher levels of organic content are possible given that organic matter levels as high as 35% are reported in backswamps characterized by water tupelo stands (Wharton et al., 1982; Sharitz and Mitsch, 1993). In our study, percentages of organic matter content in some sediment samples were greater than 50% at both sites, indicating levels of organic carbon higher than previously reported for blackwater systems.

The use of geostatistical analyses allows for not only the characterization of sediment properties on a reach scale, but also the identification of distributional patterns and patchiness (Figs. 6 and and7).7). In the headwater, 3rd order site, lowest SOD values (dark green colors) are found in conjunction with the stream center line, while highest SOD values are in floodplain areas (Fig. 6a). This is especially evident in the area-wide maximum realization (Fig. 6c) that shows high SOD rates on the eastern floodplain (brown colors). In contrast, the distributional pattern of SOD rates in the larger, 5th order site is reversed (Fig. 7). At this location, highest SOD rates (orange to red colors) are primarily found in the middle of the stream in line with the delineated stream center (Fig. 7a), while the lowest SOD rates (blue to dark green colors) are located on the margins of the stream. Again this phenomenon is shown in the most extreme using the area-wide maximum realization (Fig. 7c) since the center of the stream is dominated by high SOD (red to brown colors) with lower SOD (blue to green colors) on the periphery. Additionally, along the flow line at the 5th order location there are alternating patches of high and low values while the values along the flow line at the 3rd order site are consistently low. The alternating high and low patches of SOD in the 5th order site are likely a result of localized depressions found within the stream channel that promotes the deposition of organic matter.

The entire LREW has a low watershed gradient of 0.15%, typical of many Coastal Plain systems, but the majority of that elevation change occurs in the headwaters area. Both experimental sites have intact riparian buffers and are completely covered by overhead canopy, but with different hydrological patterns as a result of their watershed position and gradient. The gradient for the 3rd order watershed is 0.39% and the gradient from the outlet of the 3rd order site to the outlet of the 5th order site is 5.3 times less with a gradient of 0.07%. Over 63% of the total elevation change in the entire LREW is located within this 3rd order watershed. These differences in watershed gradients play a direct role in producing different flow regimes at these two locations. The 3rd order site receives overbank flooding during high flow events which causes water to spread over the floodplain and cover much of the sampled area. However, such events are typically short in duration and typically do not persist for more than a few days. Normally, the stream is within a more confined channel roughly the extent of the 10 m width designated in Fig. 3 leaving much of the mapped area dry and not an active part of instream oxygen demand. In contrast, the 5th order site is completely inundated across its entire width for a period of months at a time with high discharge events leading to increased connection of the side channels (Fig. 4). However, even during baseflow these areas remain inundated.

These differences in flow permanence and velocity likely lead to the observed spatial differences in organic carbon accumulation both within an individual reach and between sites. Prior to data collection, we hypothesized that areas of higher organic carbon (and hence SOD) would be located on the periphery at both locations since there tends to be a marked reduction in flow velocity in Coastal Plain streams as water exits the main channel and travels onto the more hydraulically dynamic floodplain (Hupp, 2000). This phenomenon is supported at the 3rd order site since the floodplain is intermittently flooded, allowing organic material to fall out of suspension. In contrast, during high discharge events the channel is subject to higher velocities allowing the flushing of organic material out of the deeper, more defined channel. Meanwhile, the 5th order site is flooded for extended periods of time, has a less defined channel, and is characterized by extensive surface roughness as a result of debris dams and multiple overstory trees growing in the flooded area. Regardless of flow intensity, the high surface roughness and lower watershed gradient likely prevents large increases in velocity even during storm events and retards the movement of organic material out of the flooded area including the stream channel. This allows for the accumulation and development of highly organic sediments. This combination of low gradients and heavy vegetative interference (surface roughness) has been seen in other blackwater ecosystems containing swamp-stream complexes. Mulholland (1981) states, “Organic carbon loading potential in Creeping Swamp is very large, primarily because of its great width and complete canopy; however its low gradient and dense vegetation enhance organic carbon retention especially of coarse particulate organic carbon by maintaining low water velocities with little erosive force, tortuous flow pathways and debris dams.”

Previous work highlighted the importance that SOD could have in determining water column DO concentrations, especially within instream swamps as they were shown to have much higher SOD rates than previously reported for southeastern Coastal Plain streams (Todd, 2008; Todd et al., 2009). This study extends the principle of point measurements of SOD to a reach scale. When looking at the average of all E-type estimates, the 5th order site had an average SOD rate of 7.13 g O2 m−2d−1(Table 1; Fig. 7). If one considers the entire mapped area being flooded it leads to an oxygen demand of 1.56 kg O2 d−1 m−1 length of stream. In contrast, the 3rd order site has a smaller area of inundated sediment available for oxygen demand. If one considers the entire mapped area being flooded, the average E-type estimate has an average SOD rate of 4.90 g O2 m−2d−1 and total oxygen demand of 0.34 kg O2 d−1 m−1 of stream length (Table 1; Fig. 6a). If only the area that is most often flooded (10 m wide strip identified in Fig. 3) is considered, the average SOD rate in the averaged realization drops slightly to 4.54 g O2 m−2d−1, but total inundated area decreases by 86%, decreasing oxygen demand by nearly 96% to 0.015 kg O2 d−1 m−1 of stream length. While the average rates of oxygen demand are high regardless of stream order, the rates per meter of stream length are over two orders of magnitude higher in the 5th order stream than in the 3rd order stream due to its greater area of constant inundation (1.56 kg O2 d−1 m−1 vs. 0.015 kg O2 d−1 m−1 of stream length). Even if one allows for complete inundation at the 3rd order site (a rare hydrologic occurrence), the 5th order site has a rate of oxygen demand per meter of stream length over four times that of the 3rd order site.

Although these two experimental locations only give examples of the effects of SOD on the watershed scale, the mapping of floodplain soils (and by relation highly organic areas most prone to high SOD) across the entire LREW shows that these locations are very common (Fig. 9). In fact, their prevalence and scale becomes more pronounced as the stream system gets larger, a situation most closely resembling the 5th order location. For instance, 89% of all stream length within the watershed is found in the smaller headwater streams (1st-3rd order). Meanwhile, the total area of floodplain soil is split approximately evenly between the headwater streams (48%) and the larger 4th-5th order stream segments (52%). However, with only 11% of the total stream length in the watershed located in 4th-5th order segments, the floodplains are disproportionately more abundant on these larger order stream segments having 7.4 times more area per stream km (21.15 ha km−1 to 2.87 ha km−1). When combined with the higher rates of oxygen demand (per unit length of stream) and longer and more sustained periods of inundation, locations such as the 5th order site appear to be common and a major driver of oxygen dynamics on the watershed scale.

5. Conclusions

Blackwater streams of the Georgia Coastal Plain are often listed as impaired due to chronically low DO levels. Previous research has shown that high SOD values, a hypothesized cause of lowered DO within these waters, are significantly positively correlated with TOC within the stream sediments (Todd, 2008; Todd et al., 2009). SOD measurements in the previous study were point measurements, making it difficult to characterize SOD values at the reach and watershed scale. However, the use of geostatistics and SGS allowed for the characterization and spatial depiction of SOD across a reach in two study locations through its relationship with TOC. The results showed TOC to be spatially autocorrelated at both experimental locations. The corresponding distribution and patchiness of SOD differed between the sites as a result of dissimilar hydrological regimes. The 5th order site, with a larger, more persistent area of inundation had higher average rates of oxygen demand both in terms of area-wide average and as a function of stream length when compared to the 3rd order site.

While only measured at two experimental locations, the mapping of floodplain soils on the watershed scale showed that areas subject to inundation are common in this watershed and that these areas are more expansive per unit stream length in larger order streams. The greater area per unit stream length in the larger order streams demonstrates the importance of areas such as the 5th order experimental area. This study highlights the importance of instream swamp areas in coastal blackwater streams and further illustrates their importance to oxygen dynamics on a watershed scale. Additionally, this research provides support for the hypothesis that many blackwater streams draining Georgia’s Coastal Plain are naturally low in DO as a result of elevated SOD.


This work has been supported by a grant from the USDA-CSREES Integrated Research, Education, and Extension Competitive Grants Program – National Integrated Water Quality Program (Award No. 2004-5113002224), Hatch and State funds allocated to the Georgia Agricultural Experiment Stations, and USDA-ARS CRIS project funds. The third author’s contribution was funded by grand R43-CA135814-01 from the National Cancer Institute. The views stated in this publication are those of the author and not necessarily represent the official views of the NCI. We would like to thank Wynn Bloodworth, Chris Clegg, Susan Crow, Herman Henry, Rodney Hill, Lorine Lewis, Andrew Mehring, and Brenda Ortiz for field, laboratory, and technical assistance. We would also like to thank Zach Aultman, Pat Pailey, and James Walker for generously allowing access onto their land.


  • Bencala KE. A perspective on stream-catchment connections. Journal of the North American Benthological Society. 1993;12 (1):44–47.
  • Benke AC, Meyer JL. Structure and function of a blackwater river in the southeastern U.S.A. Verhandlungen des Internationalen Verein Limnologie. 1988;23:1209–1218.
  • Bosch DD, Sullivan DG, Sheridan JM. Hydrologic impacts of land-use changes in coastal plain watersheds. Transactions of the ASABE. 2006;49 (2):423–432.
  • Bowie GL, et al. Rates, constants, and kinetics formulations in surface water quality modeling, EPA/600/3-85/040. 2. Office of Research and Development; Athens, GA, USA: 1985.
  • Bowman GT, Delfino JJ. Sediment oxygen-demand techniques - A review and comparison of laboratory and insitu systems. Water Research. 1980;14 (5):491–499.
  • Burke MK, King SL, Gartner D, Eisenbies MH. Vegetation, soil, and flooding relationships in a blackwater floodplain forest. Wetlands. 2003;23 (4):988–1002.
  • Cathey AM. Master’s Thesis. University of Georgia; Athens, GA: 2005. The calibration, validation, and sensitivity analysis of DOSag, an in-stream dissolved oxygen model.
  • Chau KW. Field measurements of SOD and sediment nutrient fluxes in a landlocked embayment in Hong Kong. Advances In Environmental Research. 2002;6 (2):135–142.
  • Conner WH, Buford MA. Southern deepwater swamps. In: Messina MG, Conner WH, editors. Southern Forested Wetlands: Ecology and Management. Lewis Publishers; Boca Raton, FL: 1998. pp. 261–287.
  • Crompton BJ. Master’s Thesis. University of Georgia; Athens, GA: 2005. Effect of land use on sediment oxygen demand dynamics in blackwater streams; p. 109.
  • GADNR. Water quality in Georgia, 2000–2001. Georgia Department of Natural Resources; Atlanta, GA, USA: 2000–2001.
  • GADNR. Rules and regulations for water quality control, Chapter 391–3–6. Georgia Department of Natural Resources; Atlanta, GA, USA: 2005.
  • Goovaerts P. Geostatistics for natural resources evaluation. Oxford University Press; New York, NY: 1997.
  • Goovaerts P. Geostatistical modelling of uncertainty in soil science. Geoderma. 2001;103 (1–2):3–26.
  • Hatcher KJ. Introduction to part 1: Sediment oxygen demand processes. In: Hatcher KJ, editor. Sediment oxygen demand: Processes, modeling and measurement. Institute of Natural Resources, University of Georgia; Athens, GA, USA: 1986. pp. 3–8.
  • Hupp CR. Hydrology, geomorphology and vegetation of Coastal Plain rivers in the south-eastern USA. Hydrological Processes. 2000;14:2991–3010.
  • Ice G, Sugden B. Summer dissolved oxygen concentrations in forested streams of northern Louisiana. Southern Journal of Applied Forestry. 2003;27 (2):92–99.
  • Joyce K, Todd RL, Asmussen LE, Leonard RA. Dissolved-oxygen, total organic-carbon and temperature relationships in southeastern United-States coastal-plain watersheds. Agricultural Water Management. 1985;9 (4):313–324.
  • Katz BG, Raabe EA. Suwannee River basin and estuary: An integrated watershed science program. 2005-1210. United States Geological Survey; 2005.
  • Matlock MD, Kasprzak KR, Osborn GS. Sediment oxygen demand in the Arroyo Colorado river. Journal of the American Water Resources Association. 2003;39 (2):267–275.
  • Mulholland PJ. Organic carbon flow in a swamp-stream ecosystem. Ecological Monographs. 1981;51 (3):307–322.
  • Nelson DW, Sommers LE. Total carbon, organic carbon, and organic matter. In: Sparks DL, editor. Methods of soil analysis. Part 3. Chemical Methods: SSSA Book Ser. 5. SSSA. Madison, WI: 1996. pp. 1000–1010.
  • Pringle CM, et al. Patch dynamics in lotic systems - the stream as a mosaic. Journal of the North American Benthological Society. 1988;7 (4):503–524.
  • Remy N. Geostatistical earth modeling software: User’s manual. Stanford University; 2004.
  • Rossel RAV, Goovaerts P, McBratney AB. Assessment of the production and economic risks of site-specific liming using geostatistical uncertainty modelling. Environmetrics. 2001;12 (8):699–711.
  • Rutherford JC, Wilcock RJ, Hickey CW. Deoxygenation in a mobile-bed river: I. Field studies. Water Research. 1991;25 (12):1487–1497.
  • Saito H, Goovaerts P. Geostatistical interpolation of positively skewed and censored data in a dioxin contaminated site. Environmental Science & Technology. 2000;34 (19):4228–4235.
  • Schlesinger WH. Biogeochemistry: an analysis of global change. Academic Press; San Diego, CA: 1997.
  • Sharitz RR, Mitsch WJ. Southern floodplain forests. In: Martin WH, Boyce SG, Echternacht AC, editors. Biodiversity of the southeastern United States: Lowland terrestrial communities. John Wiley & Sons, Inc; New York, NY: 1993. pp. 311–371.
  • Sheridan JM, Ferreira VA. SEWRL Research Report No. 099201. SEWRL; 1992. Physical characteristics and geomorphic data for Little River watersheds, Georgia.
  • Smock LA, Gilinsky E. Coastal plain blackwater streams. In: Hackney CT, Adams SM, Martin WH, editors. Biodiversity of the southeastern United States: aquatic communities. John Wiley and Sons, Inc; New York, NY: 1992. pp. 271–313.
  • Todd MJ. PhD Dissertation. University of Georgia; Athens, GA: 2008. Instream swamps and their effect on dissolved oxygen dynamics within blackwater streams of the Georgia coastal plain: Role of hydrology and sediment oxygen demand; p. 156.
  • Todd MJ, Vellidis G, Lowrance RR, Pringle CM. High sediment oxygen demand within an instream swamp in southern Georgia: Implications for low dissolved oxygen levels in coastal blackwater streams. Journal of the American Water Resources Association. 2009;45 (6):1493–1507.
  • USDA-NRCS. Soil Survey Geographic (SSURGO) Database. 2008.
  • Utley BJ, Vellidis G, Lowrance R, Smith MC. Factors affecting sediment oxygen demand dynamics in blackwater streams of Georgia’s coastal plain. Journal of the American Water Resources Association. 2008;44 (3):742–753.
  • Vann J, Bertoli O, Jackson S. An overview of geostatistical simulation for quantifying risk. Proceedings of Geostatistical Association of Australasia Symposium “Quantifying Risk and Error”; 2002. pp. 1–12.
  • Vannote RL, Minshall GW, Cummins KW, Sedell JR, Cushing CE. River Continuum Concept. Canadian Journal of Fisheries and Aquatic Sciences. 1980;37 (1):130–137.
  • Wetzel RG, Likens GE. Limnological analyses. Springer; New York: 2000. p. 429.
  • Wharton CH. Bulletin 114. Georgia Department of Natural Resources; Atlanta, GA, USA: 1978. The natural environments of Georgia.
  • Wharton CH, Kitchens WM, Sipe TW. The ecology of bottomland hardwood swamps of the southeast: a community profile. 81/37. U.S. Department of the Interior; Washington, DC, USA: 1982.
  • Wu RSS. A respirometer for continuous, insitu, measurements of sediment oxygen-demand. Water Research. 1990;24 (3):391–394.