Study Area

MRGB, a relatively undeveloped estuarine system located approximately 15 km north of Atlantic City in southern New Jersey (Figure 1), is one of the least disturbed estuaries in the northeastern corridor of the United States (Kennish and O'Donnell, 2002). The estuary lies entirely within the boundaries of the Jacques Cousteau National Estuarine Research Reserve (JCNERR), which is the 22nd program site of the National Estuarine Research Reserve System (Kennish, 2004). The Mullica River watershed is dominated by undeveloped forest with only 15% of the land developed or used for agricultural purposes (Zampella et al., 2001). Less than 1% of developed land occurs within JCNERR (Kennish, 2004). The estuary provides habitat and nursery grounds for several recreationally and commercially important finfish and shellfish species, as well as habitat for migrating waterfowl (Kennish, 2004).

MRGB is a shallow, well-mixed, coastal plain estuary. The estuarine system covers an area of 56 km², and it drains 87 km² of bordering salt marsh. The Mullica River drainage basin has an area of 1474 km², and it lies almost entirely within the Pinelands National Reserve. The estuary has a relatively simple geometry. The tidal freshwater region in Mullica River to the mouth of Great Bay is approximately 42 km in length (Durand, 1988). The bay is approximately 8 km long and 6 km wide, with a roughly circular shape. The average depth of the river and bay at mean low water is approximately 5 m and 2.5 m, respectively. The Mullica River drainage basin is the main source of freshwater to the bay with an annual mean discharge of approximately 29 m³ s⁻¹ (MacDonald, 1983). The bay exchanges water with the ocean through a single opening at Little Egg Inlet.

Water samples were collected at a total of 7 stations along a well-defined salinity gradient that began at the tidal freshwater region in Mullica River and ended at the mouth of Great Bay (Figure 1). In this study, the estuary was divided into four regions categorized as upper estuary (station 1–2, mean salinity = 0.80 ppt), mid-estuary (station 3–4, mean salinity = 15.4 ppt), lower estuary (station 5–6, mean salinity = 26.3 ppt), and coastal bay inlet (station 7, mean salinity = 30.3 ppt). Station 7, the coastal bay inlet, is composed of incoming flow from the immediate coastal area outside the estuary. The water composition at station 7 is similar to the coastal waters 24.1 km outside the estuary (MacDonald, 1983) and is considered representative of the nearshore coastal region.

Sampling and Analytical Methods

Physical (depth, water temperature) and chemical (salinity, pH, dissolved oxygen) water quality data from stations 2, 3, 5, and 6 (Figure 1) were obtained from the JCNERR SystemWide Monitoring Program (NERRS CDMO, 2004). Surface-water grab samples were collected at approximately 30-day intervals from February 2002 to March 2004. All grab samples were taken on the same day at slack low tide conditions (±2 hours before and after slack low tide). Replicate (n = 2) samples were collected with a bucket at an approximate depth of 10 cm. All samples were placed in amber high-density polyethylene (HDPE) bottles that were previously acid washed (15% H₂SO₄). Samples were placed on ice and processed immediately after sample collection was completed.

Sample processing for carbon and nutrient analyses entailed filtering samples through precombusted (500°C for 4 h) 47-mm Whatman glass fiber filters (GF/F) (0.7 μm) using an acid-washed, 47-mm Millipore glass filter unit with a vacuum system. Samples were placed in sterile, 15-mL Corning vials and then stored at −20°C until analyses were performed. For chlorophyll a analysis, sample processing consisted of filtering samples through precombusted (500°C for 4 h) 25-mm Whatman GF/F filters (0.7 μm) using an acid-washed, 25-mm Millipore glass filter unit with a vacuum system. Chlorophyll a sample filters were placed in acid-washed glass vials, covered with aluminum foil, and then stored at −20°C until analysis was performed.

Both DOC (Sharp et al., 1993) and TDN (Seitzinger and Sanders, 1997) were analyzed by high-temperature catalytic oxidation using a Shimadzu TOC-5000 (Shimadzu, Japan) coupled to an Antek Model 7000 Total N Analyzer with a nitric oxide chemiluminescent detector (Antek Instruments, Texas). The detection limit was 2 μM C for DOC and 1.9 μM N for TDN. DON was calculated by subtracting DIN from TDN. TDP was analyzed using the ash/hydrolysis method developed by Solórzano and Sharp (1980) with the inclusion of the 0.6 M HCl modification during the hydrolysis step (Lebo and Sharp, 1993). As recommended by Monaghan and Ruttenberg (1999), 1 mL of a 1 : 1 solution of 0.3 M MgSO₄ and 4.3 M NaCl was added to freshwater samples before the evaporation step to eliminate under-recovery of organic phosphorus. Concentrations were determined using a Lachat QuikChem AE automated ion analyzer (Lachat Inc. QuikChem Method 31-115-01-3-A, Lachat Instruments, Milwaukee, Wisconsin). The detection limit was 0.05 μM P for TDP. DOP was calculated by subtracting DIP from TDP.

Nitrite plus nitrate (NO₂⁻ + NO₃⁻, hereafter nitrate [NO₃⁻], achat, Inc. QuikChem Method 30-107-04-1-A), ammonium (NH₄⁺, Lachat Inc. QuikChem Method 31-107-06-1-A), and DIP (PO₄³⁻, Lachat Inc. QuikChem Method 31-115-01-3-A) concentrations were determined using a Lachat QuikChem AE automated ion analyzer (Lachat Instruments). The detection limit was 0.07 μM N for nitrate, 0.08 μM N for ammonium, and 0.05 μM P for DIP. DIN was calculated as the sum of nitrite plus nitrate and ammonium. Chlorophyll a was extracted from 25-mm Whatman GF/F filters with 90% acetone (EPA Method 445.0, 1997) and then measured using a Shimadzu Spectrofluorometer RF-1501 (Columbia, Maryland). Concentrations were calculated as described in Parsons, Maita, and Lalli (1984). The detection limit was 0.01μg L⁻¹ for chlorophyll a.

Biogeochemical Modeling

Carbon and nutrient budgets, fluxes, and metabolism were determined with the Land-Ocean Interactions in the Coastal Zone (LOICZ) biogeochemical model. The LOICZ model is a steady-state box model based on stoichiometrically linked water, salt, and nutrient budgets (Gordon et al., 1996). Nutrient budgets were developed using water and salt budgets to estimate the exchange of mass between each estuarine region. Nutrient sources, sinks, and fluxes are quantified from the deviation of each individual nutrient budget from expected conservative behavior (for details, see Gordon et al., 1996). Typically, LOICZ budgets are based only on dissolved nutrients because particle distributions in dynamic systems, like shallow estuaries, can be extremely heterogeneous, with potential to cause inaccurate flux estimates (Gordon et al., 1996). This LOICZ model approach is considered acceptable, given the conservation of mass law and the assumption that the nonconservative behavior of a nutrient represents either a net uptake or release from a particle (for details, see Gordon et al., 1996).

Biogeochemical processes regulating carbon, nitrogen, and phosphorus distribution, flux, and fate can be inferred from LOICZ model estimates of net ecosystem metabolism (NEM) and nitrogen metabolism, which are determined from the nonconservative flux of DIP and the C : N : P stoichiometry of organic matter (for details, see Gordon et al., 1996). In order to arrive at an accurate LOICZ model estimate of NEM and nitrogen metabolism, two assumptions must be met. First, the internal nonconservative flux of DIP must be proportional to the production and consumption of particulate organic matter. Second, abiotic processes such as DIP adsorption and desorption must be in equilibrium within the system. Consequently, if adsorption processes are significant, NEP estimates will be overestimated and denitrification will be underestimated.

For this study, the estuary was divided into a series of three lateral “boxes” or regions to account for the strong gradient in water composition from the tidal freshwater region to the estuary mouth. Each box comprised only one layer because of the weak vertical stratification in this estuary. The model boundary definition was as follows: box 1 station (1–3), box 2 station (3–4), box 3 station (4–6), and the “ocean” station (7) (Figure 1). To account for seasonal variation in river discharge, budgets were developed for each season and later averaged to create annual budgets. Seasons were defined as follows: Fall (September, October, November); Winter (December, January, February); Spring (March, April, May); and Summer (June, July, August). The estuary was assumed to be at steady-state within each season, where the change in volume (dV) over time (dt) is constant (dV/dt = 0).

The water budget was based on tributary and river discharge, direct precipitation and evaporation, and residual flow and exchange flow. Mean annual tributary and river discharge was determined using discharge data from gaged streams measured by the U.S. Geological Survey (USGS) (U.S. Geological Survey, 2004). Approximately 49% of the total watershed area is gaged. Mean annual discharge from the remaining ungaged area was estimated using gaged to ungaged watershed area ratios. Flow records for the following gaged stations were used to estimate the annual mean Mullica River basin discharge: Mullica River near Batsto (39°40′28″N, 74°39′54″W), Batsto River at Batsto (39°38′30″N, 74°39′01″W), West Branch of the Wading River near Jenkins (39°41′17″N, 74°32′53″W), Oswego River at Harrisville (39°39′48″N, 74°31′27″W), and East Branch Bass River near New Gretna (39°37′23″N, 74°26′29″W). Precipitation and temperature data (2002–2004) were obtained from the Atlantic City International Airport (40 km south of Tuckerton, New Jersey). Evaporation rate estimates were made using Hamon's Equation (Haith and Shoemaker, 1987; Hamon, 1961).

Carbon and nutrient budgets were based on inputs via tributary and river discharge, precipitation, and residual flow and exchange flow (Table 1). Concentrations of DOC, DON, DIN, DOP, and DIP in tributary and river discharge were determined from USGS 1991–2003 nutrient data. All data were categorized according to the season in which the sample was taken to obtain seasonal nutrient data for tributary and river input. Precipitation nutrient concentrations were determined from data in the literature. DOC and DON data were available from samples collected in the New Jersey Pinelands (S. Seitzinger, unpublished data). Seasonal nitrate and ammonium data were available for a station located near station 6 (Gao, 2002). Chesapeake Bay data indicate an atmospheric deposition ratio of 0.1 : 0.0 : 0.9 for dissolved inorganic, dissolved organic, and particulate phosphorus (Boynton et al., 1995). Based on this ratio, DOP atmospheric deposition was assumed to be zero. DIP seasonal data were estimated from New Jersey Pinelands atmospheric TP deposition data (Koelliker et al., 2004). The Chesapeake Bay phosphorus atmospheric deposition ratio was applied to the New Jersey Pinelands TP data to obtain an estimate of atmospheric DIP deposition. Finally, in the LOICZ model, the evaporation term for nutrient input is assumed to have a concentration of zero for all nutrient variables (Gordon et al., 1996).

Net ecosystem metabolism (NEM), defined as the difference between primary production and respiration (p-r), was estimated from the LOICZ model DIP flux multiplied by the particulate organic matter C : P ratio (Gordon et al., 1996). Nitrogen metabolism was represented by the difference between nitrogen fixation and denitrification (nfix-denit). The difference between these two processes was estimated from the sum of the LOICZ model DIN and DON flux minus the “expected” flux of nitrogen that would result from the production and consumption of organic matter. The “expected” nitrogen flux is determined from the sum of the LOICZ model DIP and DOP flux multiplied by the particulate organic matter N : P ratio (Gordon et al., 1996).

Particulate organic matter ratios used in this study to estimate (p-r) and (nfix-denit) are described below. Phytoplankton-based systems are assumed to have a C : N : P ratio of 106 : 16 : 1 (Redfield, 1958; Redfield et al., 1963). In systems with extensive freshwater marshes, brackish marshes, and salt marshes, a portion of total ecosystem primary production will result from marsh production. The upper estuary region is surrounded by both freshwater and brackish marshes (Kennish, 2004). Freshwater marsh vegetation is assumed to have a C : N : P ratio of 500 : 24 : 1 (Duarte, 1992). In the MRGB estuary, Spartina salt marshes extend upstream to station 2 (Kennish, 2004). The salt marsh species, Spartina alterniflora, in the MRGB estuary has an estimated C : P ratio of 550 and an N : P ratio of 25 (Huberty, 2005). The relative contribution of phytoplankton and marsh production to total ecosystem production in the upper and mid-estuary regions is unknown. As a result, a range for (p-r) and (nfix-denit) was calculated using C : P and N : P ratios from both a phytoplankton-based system and a marsh-based system.

In the lower estuary, ecosystem primary production is likely a function of both phytoplankton and macroalgal production. Dominant benthic macroalgae species in the lower estuary include Ulva lactuca, Enteromorpha spp., and Fucus sp.(Kennish, 2004). Systems dominated by Ulva and Enteromorpha have an assumed C : N : P ratio of 335 : 35 : 1 (Atkinson and Smith, 1983). The relative contribution of phytoplankton and macroalgal production to total ecosystem production in the lower estuary is unknown. As a result, a range for (p-r) and (nfix-denit) was calculated using C : P and N : P ratios from both a phytoplankton-based system and a mac-roalgal-based system.

Statistical Analysis

All statistical tests were performed using SigmaStat v.3.0 software. Nonparametric tests were applied when assumptions of parametric tests could not be met with either non-transformed or transformed data. Descriptive statistics (mean, median, range, standard deviation, variance) were determined for all physical, chemical, and nutrient variables. A one-way analysis of variance or Kruskal-Wallis test was applied to examine the individual effects of season and spatial location on physical and chemical variables, as well as dissolved organic and inorganic nutrient concentration. A Tukey's test or Dunn's test was performed to test for significant differences between seasons and spatial location. Correlations and regressions were used to examine associations between variables.

Water and Salt Budgets

Annual precipitation was 110.7 cm in 2002 and 122.7 cm in 2003. The highest precipitation rate was in spring (3.37 mm d⁻¹), and the lowest was in fall (2.91 mm d⁻¹). Annual evaporation was estimated to be 71.3 cm. The highest evaporation rate was in summer (4.09 mm d⁻¹), and the lowest was in winter (0.36 mm d⁻¹). The mean river flow for 2002 (18.9 m³ s⁻¹) was significantly less than the mean river flow for 2003 (32.2 m³ s⁻¹). River flow in 2002 and 2003 was 68% and 115%, respectively, of the annual mean river discharge (27.8 m³ s⁻¹) determined from long-term (75 years) USGS flow records. A seasonal flow pattern is apparent from long-term flow records where the highest river discharge occurs in spring and the lowest occurs in fall. River discharge data from 2003 follows the typical long-term seasonal pattern; however, in 2002, the seasonal pattern for river discharge exhibited an irregular pattern, with the highest river discharge in fall and the lowest river discharge in summer.

A significant difference in salinity existed between all estuarine regions and during all seasons, which is essential for accurate estimations of exchange flow between adjacent model boxes (Table 2). Residual flow and exchange flow were approximately equal on an annual time scale for the upper estuary. Annually, exchange flow dominated between the mid-estuary and lower estuary and between the lower estuary and the coastal bay inlet region. Residual flow and exchange flow were greatest in spring and lowest in fall for all regions as a result of seasonal variations in tributary and river flow, precipitation, and evaporation. Annual mean residence time in the upper estuary, mid-estuary, and lower estuary was estimated at 8.5, 3.4, and 6.2 days, respectively. All residence times were longer in fall and shorter in spring (Table 2).

Physical and Chemical Variables

Surface water temperature exhibited a seasonal cycle, with highest temperatures observed in summer and lowest in winter (Table 3). There was not a significant difference in annual mean temperature between stations (p > 0.05). Surface water salinity exhibited a strong spatial gradient of increasing salinity from station 1 to station 7. Lower salinities were typically observed at all stations during winter/spring and higher salinities during fall/summer (Table 3). Annual mean surface water salinities were higher during the low flow (18.9 m³ s⁻¹) regime in 2002 than the annual mean surface water salinities during the high flow (32.2 m³ s⁻¹) regime in 2003.

On an annual basis and during all seasons, pH increased with increasing salinity (Table 3). A significant difference in pH was found between all regions (p < 0.05), with the exception of pH between the lower estuary and coastal bay inlet (p > 0.05). Annual mean dissolved oxygen (DO) concentration was >8.9 mg L⁻¹ for all regions. No obvious spatial DO pattern was observed. The DO % saturation values increased with increasing salinity (Table 3). Annual mean DO % saturation for all regions was >90.1%.

Distribution of Organic and Inorganic Nutrients

Seasonal cycles were observed for DOC, DON, and DOP (Figure 2). In the upper estuary and mid-estuary, maximum DOC concentrations occurred in summer and minimum concentrations in late winter/early spring. A similar seasonal pattern was observed for the lower estuary and coastal bay inlet; however, these regions were less variable on an annual time scale, and the increase in DOC concentration during summer occurred to a lesser extent than the upper estuary and mid-estuary regions. The maximum DOC concentration for all regions was most often observed in June. The DON seasonal pattern was similar across all regions of the estuary, where maximum concentrations occurred in summer/fall and minimum concentrations in winter. Similar to DOC, DON in the upper estuary and mid-estuary exhibited more temporal variability than in the lower estuary and coastal bay inlet. The DOP seasonal pattern was similar across regions of the estuary, where maximum concentrations were observed in late summer/early fall and minimum concentrations in spring (Figure 2).

Spatial distributions of DOC, DON, and DOP varied greatly among each DOM component. DOC carbon decreased with increasing salinity during all seasons, with the highest concentrations observed in the upper estuary and the lowest concentrations in the coastal bay inlet (Figure 3). Distributions of DOC along the salinity gradient were generally more linear during the higher flow regime (2003) when compared to the DOC distribution during the lower flow regime (2002). The DON concentration was relatively constant along the salinity gradient during fall/winter. However, in spring/summer, DON was found to decrease with increasing salinity. The DOP concentration displayed no discernable distribution pattern according to season or region and was generally constant along the salinity gradient during all seasons (Figure 3).

Seasonal cycles were observed for nitrate, ammonium, and DIP (Figure 4). Nitrate exhibited a seasonal pattern in the upper estuary, where maximum concentrations occurred in late winter and minimum concentrations in summer/early fall. Nitrate in the mid-estuary, lower estuary, and coastal bay inlet displayed a seasonal pattern where lower concentrations were observed in spring/summer and higher concentrations in fall/winter. Ammonium in the upper estuary was relatively constant over temporal scales with no apparent seasonal pattern. The mid-estuary, lower estuary, and coastal bay inlet temporal ammonium concentration distribution was more variable compared to nitrate. Ammonium concentrations generally increased in fall, decreased in winter, and increased again to a lesser extent in early spring/summer. DIP in the upper estuary was less variable over temporal scales compared to the other regions, displaying a relatively constant pattern throughout an annual cycle. A distinct DIP seasonal pattern was observed for the mid-estuary, lower estuary, and coastal bay inlet, where maximum concentrations occurred in late summer/early fall and minimum concentrations occurred in spring. Maximum DIP concentrations were most often observed in August and September (Figure 4).

Spatial distributions of nitrate, ammonium, and DIP varied greatly among each inorganic nutrient species (Figure 5). Nitrate concentration decreased with increasing salinity down-estuary. Ammonium concentration distribution varied along the salinity gradient during most seasons with the exception of fall, when ammonium concentration generally increased with increasing salinity until it peaked in the lower estuary and then decreased in the coastal bay inlet. The concentration of DIP generally increased with increasing salinity along the salinity gradient. During fall, DIP exhibited a similar spatial pattern to ammonium (Figure 5).

Chlorophyll a

Chlorophyll a exhibited a seasonal pattern with maximum concentrations occurring during summer/early fall, and minimum concentrations during winter in all estuary regions (Figure 4). Spring exhibited episodes of high and low chlorophyll a concentrations. Chlorophyll a concentrations displayed no discernable distribution pattern along the salinity gradient (Figure 5). However, peaks of chlorophyll a most often occurred in the upper estuary at station 2, in the mid-estuary at station 4, in the lower estuary at station 5, and in the coastal bay inlet at station 7.

Nutrient Fractions and Ratios

The annual mean fraction of DON in TDN was 51% for the upper estuary, 59% for the mid-estuary, 68% for the lower estuary, and 76% for the coastal bay inlet. The fraction of DON in TDN increased with increasing salinity during winter/spring, whereas in fall/summer the DON fraction was relatively constant along the salinity gradient (Table 4). The annual mean fraction of DOP in TDP was 67% for the upper estuary, 44% for the mid-estuary, and 28% for both the lower estuary and coastal bay inlet. The DOP fraction dominated TDP in the upper estuary region during all seasons, ranging from 56% in winter to 82% in summer (Table 4). The annual mean DOC : DON ratio was 48 for the upper estuary, 34 for the mid-estuary, 24 for the lower estuary, and 19 for the coastal bay inlet. The annual mean DOC : DOP ratio was 2011 for the upper estuary, 1894 for the mid-estuary, 1316 for the lower estuary, and 1029 for the coastal bay inlet. The annual mean DON : DOP ratio was 46 for the upper estuary, 58 for the mid-estuary, and 54 for both the lower estuary and coastal bay inlet.

Nutrient Correlations

River discharge was significantly correlated with upper estuary DOC (r = 0.32, p < 0.05) and nitrate (r = 0.19, p < 0.05) concentration; however, no other nutrient variables were significantly correlated with river discharge. A significant negative correlation was also found between salinity and DOC (r = 0.78, p < 0.05). Chlorophyll a and river discharge displayed significant negative correlations in the upper estuary (r = −0.64, p < 0.05), mid-estuary (r = −0.61, p <0.05), and lower estuary (r = −0.52, p < 0.05). A significant positive correlation was determined between DOP and chlorophyll a at station 6 (r = 0.57, p < 0.05). No other significant correlations were found between chlorophyll a with either DOC, DON, or DOP at any other station. Ammonium and DON exhibited significant positive correlations at station 5 (r = 0.56, p < 0.05) and station 6 (r = 0.50, p < 0.05). Both DIP and DOP also exhibited significant positive correlations at station 5 (r = 0.74, p < 0.05) and station 6 (r = 0.71, p < 0.05).

Nutrient Budgets and Fluxes

The annual mean DOC input from the Mullica River watershed was 374 ×10⁶ mol y⁻¹ (0.25 mol m⁻² y⁻¹) (Table 5). The largest DOC input (1324 ×10³ mol d⁻¹) from the watershed occurred in spring, and the smallest input (741 × 10³ mol d⁻¹) occurred in fall. Annual export of DOC from the lower estuary to the nearshore coastal region was 679 × 10⁶ mol y⁻¹ (Table 5). The annual mean watershed contribution of DON and DIN to the estuarine system was 12 × 10⁶ mol y⁻¹ (0.01 mol m⁻² y⁻¹) and 27 × 10⁶ mol y⁻¹ (0.02 mol m⁻² y⁻¹), respectively. The largest watershed DON input occurred during summer (47 × 10³ mol d⁻¹), whereas the largest DIN input occurred during winter (105 × 10³ mol d⁻¹). Annual mean exports of DON and DIN from the lower estuary to the nearshore coastal region were 17.4 × 10⁶ mol y⁻¹ and 15.7 × 10⁶ mol y⁻¹, respectively (Table 5). The annual mean DOP and DIP input from the watershed was 0.36 ×10⁶ mol y⁻¹ (0.24 mmol m⁻² y⁻¹) and 0.15 × 10⁶ mol y⁻¹ (0.10 mmol m⁻² y⁻¹), respectively. Input of DOP from the watershed was greatest in spring (1.64 × 10³ mol d⁻¹), whereas DIP input was greatest in the spring/summer (0.46–0.48 × 10³ mol d⁻¹). The annual mean export of DOP and DIP from the lower estuary to the nearshore coastal region was 0.44 × 10⁶ mol y⁻¹ and 0.76 × 10⁶ mol y⁻¹, respectively.

Annually, the entire estuarine system was a net source of DOC, DON, and DIP, a net sink of DIN, and in approximate balance of DOP (Table 6). Overall, the upper estuary and mid-estuary generally acted as a net sink for most nutrients, whereas the lower estuary acted as a net source. Seasonally, the upper estuary and lower estuary were net sources of DOC during all seasons, with the largest net source in winter/summer and the smallest in fall/spring. The mid-estuary was a net source of DOC in fall/summer and a net sink in winter/spring. The upper estuary was a source of DON in fall/spring and a sink in winter/summer. The mid-estuary and lower estuary were net sources of DON, with the exception of winter in the lower estuary. The largest net source of DON in the mid-estuary occurred in fall, whereas the largest source of DON in the lower estuary occurred in summer. The upper estuary and lower estuary were net sinks of DOP in winter/spring and a net source during fall/summer. The largest source of DOP (∼0.05 mmol m⁻² d⁻¹) occurred in summer in both the upper estuary and lower estuary. The mid-estuary was a net sink of DOP in fall/summer and a net source in winter/spring. The largest DOP removal occurred in summer in the mid-estuary (Table 6).

During all seasons a net removal of DIN was observed in the upper estuary, mid-estuary, and lower estuary with the exception of summer in the mid-estuary and fall in the lower estuary (Table 6). The largest DIN sink occurred in the mid-estuary during winter (−11.30 mmol m⁻² d⁻¹). In the lower estuary, the winter, spring, and summer seasons exhibited a fairly constant sink of DIN in contrast to the large source of DIN during the fall. The upper estuary and mid-estuary were net sinks of DIP during all seasons with the exception of winter in the upper estuary and spring in the mid-estuary. The largest DIP removal occurred in the mid-estuary during summer (−0.250 mmol m⁻² d⁻¹) (Table 6). The lower estuary was a net source of DIP during all seasons, with the largest source of DIP (+0.113 mmol m⁻² d⁻¹) occurring in fall.

The LOICZ model results indicate the entire estuarine system was net heterotrophic. Annual mean NEM (p-r) estimates for the upper estuary and mid-estuary suggest that these regions were net autotrophic (Table 6). In contrast, the lower estuary appears to be net heterotrophic with annual mean NEM (p-r) estimates ranging from −1.8 to −5.8 mol C m⁻² y⁻¹ (−22 to −70 g C m⁻² y⁻¹). The LOICZ model results also indicate that the lower estuary increased in heterotrophy from winter to fall (Table 6). Annual mean nitrogen metabolism (nfix-denit) estimates for the entire estuarine system suggest that denitrification was greater than nitrogen fixation by −0.4 to −0.5 mol N m⁻² y⁻¹. Annual mean (nfix-denit) values indicate that net denitrification (−1.3 to −1.5 mol N m⁻² y⁻¹) occurs in the upper estuary, and net nitrogen fixation (+0 to +0.5 mol N m⁻² y⁻¹) occurs in the mid-estuary. Annual mean (nfix-denit) estimates for the lower estuary ranges from −0.2 to −0.6 mol N m⁻² y⁻¹, which suggests net denitrification occurs in this region.

Water and Salt Budgets

Water and salt budget results corroborate well with hydrological estimates from other studies on the MRGB estuarine system. Annual evaporation was estimated to be 71.3 cm for the study area. Rhodehamel (1998) and Dow and Dewalle (2000) have estimated annual evaporation in southern New Jersey at >66 cm and 71 cm, respectively. Annual mean Mullica River Basin discharge (27.8 m³ s⁻¹) during this study was similar to MacDonald's (1983) 29.4 m³ s⁻¹ estimate, which was calculated using the same watershed area ratio method used in this study. Annual mean residence time for the river section (upper estuary and mid-estuary) was 11.9 days, which is approximately equal to MacDonald's (1983) residence time estimate of 12 days, based on a freshwater fraction method. In contrast, the lower estuary annual mean residence time of 6.2 days differed from MacDonald's (1983) estimate of 9 days. The discrepancy between the lower estuary residence time estimates can be attributed to differences in the coastal area salinity value used to calculate the residence time; this study used a mean coastal area salinity of 30.3 ppt and MacDonald (1983) used 31.5 ppt. Using MacDonald's (1983) freshwater fraction method and a coastal salinity of 30.3 ppt, a lower estuary residence time of 6.4 days is calculated, which is approximately equal to the residence time estimated with the LOICZ model.

Nutrient Distribution, Budgets, and Fluxes

The MRGB estuarine system exhibited a DOM seasonal cycle to varying degrees in each region of the estuary. The entire estuarine system was a net source of DOC, DON, and DIP, a net sink of DIN, and in approximate balance of DOP (Table 6). The annual NEM estimate for the entire estuarine system indicates that this system is net heterotrophic, consuming 1.1–1.4 mol C m⁻² y⁻¹. The lower estuary acted as a net nutrient source and was determined to be net heterotrophic, consuming 1.8–5.8 mol C m⁻² y⁻¹. In contrast, the upper estuary and mid-estuary generally acted as a net sink for all nutrients, with the exception of DOC in the upper estuary and DON in the mid-estuary (Table 6). The NEM estimates indicate that the upper estuary and mid-estuary regions are net autotrophic, producing 1.1–5.7 mol C m⁻² y⁻¹ and 3.3–16.8 mol C m⁻² y⁻¹, respectively. The nutrient cycling trends exhibited by the entire MRGB estuarine system seem to follow Odum's (1980) often cited “outwelling” hypothesis, which asserts that marsh-estuarine systems produce more organic material than can be degraded or stored within the system, and consequently, the excess nutrients are exported to the coastal ocean. Additionally, nutrient cycling in this estuarine system also appears to follow the marsh-estuarine continuum theory proposed by Dame (1994), which suggests that net nutrient fluxes are a function of marsh-estuarine system geological age and developmental status, where the upper estuarine regions are assumed to be at an earlier stage of development and consequently function as net nutrient importers, whereas the lower regions are assumed to be at a more mature stage of development and function as net nutrient exporters. Other estuarine systems, such as the North Inlet Estuary (South Carolina) (Dame, 1994; Williams et al., 1992) and Fier d'Ars Estuary (France) (Bel Hassen, 2001), have been shown to follow a similar nutrient cycling trend, where net import occurs in the upper regions with extensive salt marshes and net export occurs from the coastal bays.

The following discussion addresses potential sources, sinks, and biogeochemical processes responsible for regulating the distribution, flux, and fate of dissolved organic and inorganic nutrients in the MRGB estuarine system. Evidence from past MRGB marsh-estuarine nutrient cycling studies (Durand, 1988), in addition to nutrient distribution, budget, and flux results from this study, indicate that the marshes (freshwater, brackish, and salt) may serve a significant role in the processes that are important to nitrogen and phosphorus cycling in the MRGB estuary.

Annual mean watershed DIN input (27 10⁶ mol y⁻¹) was approximately two times greater than DON input (12 × 10⁶ mol y⁻¹); however, export to the nearshore coastal region was approximately equal (15.7–17.4 10⁶ mol y⁻¹), indicating that a portion of the DIN in the TDN pool is converted to DON during transport through the estuarine system. In contrast, annual mean watershed DOP input (0.36 × 10⁶ mol y⁻¹) was approximately two times greater than DIP input (0.15 × 10⁶ mol y⁻¹). However, annual mean export of DIP (0.76 × 10⁶ mol y⁻¹) and DOP (0.44 × 10⁶ mol y⁻¹) was approximately 63% and 37%, respectively, of the total TDP exported from the lower estuary to the nearshore coastal region, indicating that a change also occurred in the organic and inorganic fractions of the TDP pool during transport through the estuarine system.

In the MRGB estuarine system, several studies were conducted during the 1980s to elucidate nutrient flux relationships between the marsh and estuary. Durand (1988) found that the marsh appears to export DON and DIP and import nitrate and ammonium. Durand (1988) also determined that all nitrate and ammonium delivered to the marsh was taken up by the marsh, which suggests the marsh vegetation may be nitrogen limited. From these studies, Durand (1988) was unable to determine whether the source of DON could be attributed to subsurface marsh water or biological activity. Durand (1988) attributed marsh DIP export to S. alterniflora actively pumping phosphorus from the sediments to the overlying water.

Variation in MRGB estuary nutrient distributions and fluxes is likely attributed to a combination of river and tidal hydrodynamics, watershed inputs, and estuarine biological production. During this study, a DOM seasonal cycle was observed to varying degrees in each region of the estuary. Throughout the growing season, DOC, DON, and DOP concentrations increased from spring to fall, with maximum concentrations occurring in summer/early fall (Figure 2). The largest DOM export from the lower estuary to the coastal bay inlet also occurred during summer (Table 6). Potential sources and processes controlling the seasonal cycle and export of DOM include terrestrial inputs from vegetation and soil, marsh-estuarine nutrient cycling, and macroalgae and phytoplankton production.

C : N : P stoichiometry can be frequently used to determine DOM source composition (Goni and Thomas, 2000; Mannino and Harvey, 2000). Generally, DOM with high C : N and C : P ratios originate from terrestrial and marsh sources, whereas DOM with lower C : N and C : P ratios originate from macroalgae and plankton sources (Table 7). High DOC : DON (annual mean 24–48) and DOC : DOP (annual mean 1316–2011) ratios found within the estuary suggest that a portion of DOM originates from both terrestrial and marsh sources, with terrestrial sources being most significant in the upper estuary region where the highest C : N ratios were observed. In contrast, the lower estuary and coastal bay inlet tend to have lower C : N and C : P ratios in the summer/fall (Table 4), which suggests that plankton (phytoplankton and bacteria) and macroalgal production may contribute a significant portion of DOM to these regions during the growing season.

Terrestrial and marsh derived DOM inputs may be significant in regulating the observed DOM seasonal cycle within this estuarine system, with the most dominant contributions occurring in the upper estuary and the least dominant in the lower estuary. In this study, seasonal watershed DOM inputs were found to occur, with the largest DOC and DOP input occurring in spring and the largest DON input occurring in spring/summer (Table 5). In addition, the following correlations and distributions indicate that upper estuarine DOM is probably derived from both terrestrial and marsh sources: a significant positive correlation between river discharge and upper estuary DOC concentration (r = 0.32, p <0.05), a significant negative correlation between salinity and DOC (r = −0.78, p < 0.05), and an observed DON decrease with increasing salinity during spring/summer. Several other studies have also shown a positive correlation between river discharge and DOC concentration and suggest an upstream DOC source (Cauwet and Meybeck, 1987; Moore, 1989; Schlesinger and Melack, 1981) in addition to an increase in DOC leaching from soil and plant litter during periods of high discharge (Goni and Gardner, 2003). DOC export from freshwater marshes within the Chesapeake Bay (Virginia) estuarine system (Neubauer, Miller, and Anderson, 2000) and North Inlet Estuary (South Carolina) salt marshes (Goni and Gardner, 2003) have also been implicated as potentially significant DOC sources. Past MRGB marsh-estuary nutrient flux studies have also shown that the marsh functions as a DON source (Durand, 1988).

The additional influence of macroalgae and phytoplankton production is likely significant in controlling the temporal distributions of DOM in the lower estuary, as well as the potential import of phytoplankton-produced DOM from the coastal bay inlet during tidal exchange. In the lower estuary, the release of DOM to the water column from benthic mac-roalgae during growth and decay may serve as an additional seasonal DOM source. Dominant forms of macroalgae in the MRGB estuary include Ulva lactuca, Enteromorpha spp., and Fucus sp. (Kennish, 2004). Macroalgae are prevalent at stations 5 and 6 throughout the growing season and then begin to senesce in early fall (personal observation). In Hog Island Bay (Virginia), a shallow, coastal lagoon system, U. lactuca was found to release DON to the water column during both active growth and senescence; the largest DON and DIN sediment flux also occurred after a macroalgal bloom decline (Tyler, McGlathery, and Anderson, 2001). In the MRGB estuary, maximum DON and ammonium concentrations occurred in September of both sample years. It is possible that the decomposition of macroalgae in fall may account for a portion of the high DON and DIN water column concentrations observed in September and the large DIN flux to the nearshore coastal region in the fall (Table 6).

A portion of the DOM source, especially in the lower estuary and coastal bay inlet regions, is likely attributed to the seasonal production of phytoplankton DOM. Studies have shown the production of DOM during the growing season occurs through several processes such as phytoplankton exudation (Baines and Pace, 1991; Sondergaard et al., 2000), cell lysis (Agusti et al., 1998), and release from grazing (Bronk et al., 1998; Bronk and Ward, 1999). Similar temporal patterns of DOM have also been observed in coastal systems where phytoplankton production was implicated in the seasonal cycle (Amniot, Sayed, and Kerouel, 1990; Bronk et al., 1998; Butler et al., 1979; Wong et al., 2002). In this study, coincident with the increase and peak in DOM, maximum chlorophyll a concentrations in all regions of the estuary were found to occur during summer/early fall with episodes of high and low chlorophyll a concentrations during spring. In addition, at station 6 in the lower estuary, a significant correlation between chlorophyll a and DOP (r = 0.57, p < 0.05) was determined. Significant correlations between DON and ammonium (r ≥ 0.50, p < 0.05) and DOP and DIP (r ≥ 0.71, p < 0.05) were also determined in the lower estuary, which suggests that a tight-coupling exists between production and consumption of organic matter in this region.

Potential biogeochemical processes controlling MRGB estuarine DOM and inorganic nutrient sources were described in detail in the discussion above. The section below provides a brief discussion on potential biogeochemical processes controlling MRGB nutrient sinks. Removal processes of DOC include flocculation (Forsgren and Jansson, 1992; Fox, 1983; Sholkovitz, 1976; Sholkovitz, Boyle, and Price, 1978), photochemical degradation (Mopper et al., 1991), and heterotrophic consumption (Raymond and Bauer, 2000; Wikner, Cuadros, and Jansson, 1999). In the MRGB estuary, the most significant DOC removal process is most likely attributed to flocculation of the humic acid component (Forsgren and Jansson, 1992; Fox, 1983). In the mid-estuary, approximately 13% and 12% of DOC was lost during winter and spring, respectively. The Mullica River watershed drains upland areas consisting of cedar swamps and sphagnum and cranberry bogs, which results in dark brown river water with high humic and fulvic acid content (Kennish, 2004). The percent humic acid component in the Mullica River DOC pool is currently unknown; however, the humic acid fraction typically comprises 3% to 11% of the DOC pool (Sholkovitz, 1976). The Mullica River DOC pool may have an even higher humic acid content that is closer to the “black-water” rivers in southeastern coastal plain estuarine systems, which contain approximately 75% humic acid content (Beck, Reuter, and Perdue, 1974).

Heterotrophic and autotrophic utilization may account for a portion of the observed DON sink (Berman and Chava, 1999; Seitzinger, Sanders, and Styles, 2002; Stepanauskas and Leonardson, 1999). The MRGB marshes appear to function as significant DIN sinks (Durand, 1988). Estimates of (nfix-denit) (Table 6) also indicate denitrification may be more significant as a potential DIN sink in the upper estuary than in the mid-estuary and lower estuary. A significant portion of DIP removal in the upper estuary and mid-estuary may be attributed to abiotic processes such as sorption onto suspended sediments (Fox, 1993; Froelich, 1988) rather than being solely accounted for by the production and consumption of particulate organic matter. Fang (2000) found indications of DIP removal onto suspended particulate matter in the salinity range of 5–10 ppt in the Tanshui Estuary (Northern Taiwan). In the Delaware Estuary, Lebo and Sharp (1992) also found evidence of abiotic DIP removal in the river section of the estuary. In the MRGB estuary, mean turbidity levels range from 5 NTU to 32 NTU, with higher turbidity levels in the river than in the bay (Kennish, 2004). Studies during 1961–1962 and 1978–1979 in the MRGB estuary suggest there is a strong relationship between turbidity and phosphorus concentration in the upper estuary, with maximum DIP concentrations in areas of maximum turbidity (Durand, 1988).