Articles | Volume 17, issue 14
Research article
22 Jul 2020
Research article |  | 22 Jul 2020

Relative impacts of global changes and regional watershed changes on the inorganic carbon balance of the Chesapeake Bay

Pierre St-Laurent, Marjorie A. M. Friedrichs, Raymond G. Najjar, Elizabeth H. Shadwick, Hanqin Tian, and Yuanzhi Yao

The Chesapeake Bay is a large coastal-plain estuary that has experienced considerable anthropogenic change over the past century. At the regional scale, land-use change has doubled the nutrient input from rivers and led to an increase in riverine carbon and alkalinity. The bay has also experienced global changes, including the rise of atmospheric temperature and CO2. Here we seek to understand the relative impact of these changes on the inorganic carbon balance of the bay between the early 1900s and the early 2000s. We use a linked land–estuarine–ocean modeling system that includes both inorganic and organic carbon and nitrogen cycling. Sensitivity experiments are performed to isolate the effect of changes in (1) atmospheric CO2, (2) temperature, (3) riverine nitrogen loading and (4) riverine carbon and alkalinity loading. Specifically, we find that over the past century global changes have increased ingassing by roughly the same amount (∼30 Gg-C yr−1) as has the increased riverine loadings. While the former is due primarily to increases in atmospheric CO2, the latter results from increased net ecosystem production that enhances ingassing. Interestingly, these increases in ingassing are partially mitigated by increased temperatures and increased riverine carbon and alkalinity inputs, both of which enhance outgassing. Overall, the bay has evolved over the century to take up more atmospheric CO2 and produce more organic carbon. These results suggest that over the past century, changes in riverine nutrient loads have played an important role in altering coastal carbon budgets, but that ongoing global changes have also substantially affected coastal carbonate chemistry.

1 Introduction

The well-documented rise in atmospheric CO2 concentrations is one of the most ubiquitous changes in global biogeochemical cycling over the past century (e.g., Keeling et al.2003). Although the ocean's biological pump maintains atmospheric CO2 significantly lower than it would otherwise be, the uptake of anthropogenic CO2 by the ocean is governed largely by chemical and physical processes. These processes include the diffusion of CO2 across the air–sea interface, the dissolution of CO2 and its dissociation into bicarbonate and carbonate ions, and the transport of anthropogenic dissolved inorganic carbon into the ocean interior by vertical mixing and subduction. Thus, early estimates of the uptake of anthropogenic CO2 by the ocean did not explicitly include marine biological processes (Oeschger et al.1975). However, if biological processes change during the uptake of anthropogenic CO2 into the ocean, then they can alter that uptake. Such changes could occur in at least three ways: (1) the CO2 invasion itself, which could influence photosynthesis and calcification (see Riebesell et al.2007, 2000); (2) climate change, which could influence biogeochemistry via warming and changes in mixing and advection (e.g., Sarmiento et al.1998); and (3) the delivery of nutrients and carbon via river runoff and atmospheric nitrogen deposition (Da et al.2018; Duce et al.2008; Ver et al.1999; Walsh et al.1981). Coastal regions, especially estuaries, have unique susceptibility to changes due to their proximity to anthropogenic nutrient and carbon sources and therefore may be particularly important for understanding how biological processes may influence the uptake of anthropogenic CO2 by the ocean.

Different perspectives on the role of the coastal ocean in the uptake of anthropogenic CO2 have been proposed over the past decades (see Cai2011, for a review). For example, Walsh et al. (1981) argued that the input of anthropogenic nitrogen to the ocean by rivers has stimulated primary production and enhanced the ocean's uptake of atmospheric CO2. On the other hand, Ver et al. (1999) found that increases in the riverine input of organic matter to the ocean have had a larger and counteracting effect by stimulating heterotrophy. Such disagreements are to be expected given the great heterogeneity of coastal waters and the differences in their dominant processes.

Process-based biogeochemical models afford the opportunity to isolate the various ways in which the exchange of CO2 between the atmosphere and coastal waters has changed during the industrial period. Such models represent many of the important forcing mechanisms, such as the essentially global changes of increasing temperature and atmospheric CO2, as well as regional shifts in the delivery of freshwater, nutrients, carbon and alkalinity by rivers. Despite the considerable advancement of estuarine biogeochemical models in recent years (Ganju et al.2016), the relative impact of these global and regional changes on carbon cycling in coastal waters is not always clear. Here, we examine these changes and quantify them in the context of the Chesapeake Bay.

The Chesapeake Bay is a coastal-plain estuary and the largest estuary in the continental United States. Its watershed provides ∼80 km3 yr−1 of freshwater with nearly half of this input coming from one river positioned at the northern end of the bay (the Susquehanna River; Fig. 1). At its southern end, the bay is in direct contact with the shelf water of the Mid-Atlantic Bight (Fig. 1). This configuration leads to a meridional gradient of salinity but also of dissolved inorganic carbon (DIC) and total alkalinity (TA) (e.g., Shen et al.2019a; Friedman et al.2020). The gradient is apparent throughout the year, although the seasonal discharge of the rivers (maximum around March–April) modulates the salinity, DIC and TA (e.g., Shadwick et al.2019b), especially in the northern part of the bay (see the observations in Brodeur et al.2019).

Figure 1Map of the study area with the key tributaries labeled. The gray shading represents the model grid cells (see Da et al.2018, for a map of the full model domain). Red circles represent the locations of riverine inflow in the model (10 rivers total). The largest rivers (Susquehanna, Potomac, York and James rivers) have their discharge distributed over two neighboring model grid cells (i.e., two overlapping red circles in the figure). Blue circles represent locations where riverine alkalinity and DIC data are available. Each location is identified with an eight-digit number (see Methods section).

To better understand the evolution of the bay over the last century, we perform a process-oriented study based on a numerical model of the Chesapeake Bay. The study includes a set of numerical experiments quantifying the sensitivity of the inorganic carbon budget to the global and regional changes described above. The paper is structured as follows. The modeling system and the numerical experiments are described in the next section. The results from the control experiment (years 2000–2014) are then presented, compared to observations and contrasted with sensitivity experiments representative of the period 1900–1914. Finally, the results of the study are discussed in the context of the existing literature and of the ongoing global and regional changes impacting the Chesapeake Bay region.

Table 1List of experiments conducted in the study (see Methods section for their description) with their differences from the control experiment highlighted in bold.

Download Print Version | Download XLSX

2 Methods

The study uses a numerical model of the Chesapeake Bay (Feng et al.2015; Irby et al.2018; Da et al.2018, with modifications described below) and includes a total of six numerical experiments (Table 1). The first experiment (control experiment) represents contemporary conditions with realistic forcings for a period of 15 years (2000–2014). Then, four sensitivity experiments are used to isolate the effect of specific parameters on the inorganic carbon balance: atmospheric CO2 concentration, temperature, and riverine inputs of nitrogen, carbon and alkalinity. In each of those four experiments, the parameter of interest is modified to represent the period 1900–1914 while keeping all other components of the model the same as in the control experiment (Table 1). The relative importance of global and regional changes is assessed by combining the experiments as follows: 1900CO2+1900T (global) and 1900N+1900C (regional; see Table 1). The last of the six experiments includes the four perturbations at once to evaluate potential synergies. All sensitivity experiments are preceded by a 1-year period during which the model solution adjusts itself to the modification. This adjustment period is not part of the 15-year-long experiments (it precedes them) so that all the experiments represent the bay in a stationary state (trends ≈0).

2.1 Control experiment (2000–2014)

2.1.1 Estuarine model

The numerical experiments are based on an implementation of the Regional Ocean Modeling System (ROMS, Shchepetkin and McWilliams2005) for the Chesapeake Bay (ChesROMS-ECB; see Da et al.2018; Feng et al.2015). The model domain includes the bay and a portion of the continental shelf (Fig. 1) with a curvilinear discretization on the horizontal (resolution O(1 km) in the bay) and 20 topography-following levels on the vertical (Xu and Hood2006). The model domain assumes permanent coastlines and thus no flooding of land areas nor drying of shoals.

The ROMS physical kernel is coupled to a biogeochemical module (Estuarine Carbon Biogeochemistry, ECB) at every baroclinic time step (60 s) using a positive-definite advection scheme (Smolarkiewicz and Margolin1998). The module represents the nitrogen and carbon cycles of the lower trophic levels (Druon et al.2010) with additional processes specific to estuarine systems (see Da et al.2018; Feng et al.2015). The ECB module includes 17 state variables: nitrate (NO3-), ammonium (NH4+), oxygen, inorganic suspended solids (ISSs), dissolved inorganic carbon (DIC), total alkalinity (TA), phytoplankton, chlorophyll, zooplankton, small and large nitrogen and carbon detritus, and separate semilabile and refractory dissolved organic carbon and nitrogen components (DOC and DON). Hereafter we refer to the sum of nitrate and ammonium as dissolved inorganic nitrogen (DIN). Note that ECB does not represent the oxidation of hydrogen sulfide (see Cai et al.2017). The equation for each state variable is documented in the Supplement (Tables S3–S6).

A number of modifications are made to the ECB module described in Da et al. (2018). Specifically, the parameters controlling the growth and fate of phytoplankton are modified to better represent the observed seasonal cycle of the bay. First, the initial slope of the photosynthesis–irradiance curve is set to 0.04 (W m−2 d)−1, similar to the “spring phytoplankton group” of Cerco and Noel (2004). For T>20C, the maximum phytoplankton-specific growth rate is set to 0.6 exp (0.078 T) d−1, where T is the water temperature in degrees Celsius and the coefficient 0.078 C−1 is from Lomas et al. (2002). A constant rate of 2.15 d−1 is assumed when T<20C (as in Feng et al.2015; Da et al.2018) to reflect the observed temperature independence in this range (Lomas et al.2002). The phytoplankton mortality rate is also decreased to 0.05 d−1 and the aggregation rate is increased to 0.008 (mmol-N m−3 d)−1 to better represent the nonzero phytoplankton concentrations observed at depth during the winter period. Finally, a minimum value of 0.6 m−1 is enforced for the coefficient of diffuse attenuation to represent the effect of ISS resuspension in the lower part of the bay. All these model parameters are documented in the Supplement.

2.1.2 Atmospheric forcing for the control experiment (2000–2014)

The model is forced with the atmospheric forcings (North American Regional Reanalysis, NARR; Mesinger et al.2006) described in Da et al. (2018). In addition, we assume that atmospheric CO2 concentrations vary slowly over the period 2000–2014 with a mixing ratio represented by a quadratic polynomial:

(1) mixing ratio = 371.19 + 1.86 t - t 0 + 0.0125 t - t 0 2 ,

where t is the time in years and t0=2001. The coefficients are based on a fit to historical global values from the period 1950–2011 assembled by Miller et al. (2014) (see Fig. 2). Seasonal variations in atmospheric CO2 concentrations are not considered given our focus on long-term changes.

Figure 2CO2 mixing ratio used in the numerical experiments (green; see Methods section). The historical values (blue) are from Miller et al. (2014). The dashed line represents the constant value assumed in the experiments of 1900–1914.


At the air–sea interface, the model calculates CO2 fluxes using

(2) F = k w α p CO 2 a - p CO 2 w ,

where F is the flux (mmol-C m−2 d−1), kw is the transfer velocity for CO2 (m d−1) (Wanninkhof1992, his Eq. 3), α is the CO2 solubility in seawater (mmol-C m−3µatm−1;  Weiss1974), pCO2a is the atmospheric partial pressure of CO2, and pCO2w is the partial pressure of CO2 at the water surface. Note that F is defined as positive for ingassing; we use this convention because the carbon budget of the bay is being assessed and all carbon sources are treated as positive. An algorithm adapted from Zeebe and Wolf-Gladrow (2001) is applied to compute pCO2w (as in Fennel et al.2008) using modeled surface temperature, salinity, DIC and TA at each model time step (60 s). The algorithm uses the dissociation constants from Mehrbach et al. (1973) as fitted by Millero (1995).

2.1.3 Oceanic forcing for the control experiment (2000–2014)

Oceanic conditions are prescribed at the model open boundary positioned on the continental shelf of the Mid-Atlantic Bight. Temperature, salinity, oxygen, and dissolved nitrogen (organic and inorganic) are derived using a combination of climatology, in situ (i.e., observational) data and satellite data, as described in Da et al. (2018). For TA and DIC, data from 12 cruises conducted between 2005 and 2006 in the vicinity of the bay's mouth (Filippino et al.2009, 2011) were used to derive the following relationships with salinity (S):


where TA is in milliequivalents per cubic meter, DIC is in millimoles of carbon per cubic meter, N is the number of measurements and R2 is the coefficient of determination. These relationships are combined with the seasonal climatology used for salinity to prescribe TA and DIC at the model open boundary. The pH at the oceanic model boundary, calculated from these TA and DIC values, varies seasonally and spatially within the range 7.75<pH<8.05 with an average value pH=7.89 (total scale). This range is consistent with the measurements in Wang et al. (2013) (their Fig. 8b, transect “MA”, pH7.9±0.1 where ± represents 1 standard deviation). Note that the same oceanic conditions are used in the 1900–1914 and 2000–2014 experiments since we are primarily interested in historical changes that occurred inside the bay and at its surface (e.g., atmospheric CO2). The potential impact of the historical change in DIC on the continental shelf (i.e., the anthropogenic DIC) is thus not represented here, but it should be considered in future studies.

2.1.4 Riverine forcing for the control experiment (2000–2014)

At the land–estuary interface, the model is linked to the Dynamic Land Ecosystem Model (DLEM; Yang et al.2015b, a; Tian et al.2015) as in Feng et al. (2015). The version of DLEM used here has a resolution of 4 km and provides daily fluxes for 1900–2015 for the entire watershed of the bay. For this study these fluxes (defined at the coastlines of the bay) are aggregated into 10 river sources positioned along the bay (Fig. 1) and include freshwater, NO3-, NH4+, DON, DOC, and particulate organic nitrogen and carbon (PON and POC). Riverine fluxes of ISS are provided by the Chesapeake Bay watershed model (Shenk and Linker2013).

Riverine fluxes of DIC and TA are calculated from the freshwater discharge of DLEM coupled with our best estimates of riverine concentrations. Numerous studies have shown that riverine TA and DIC exhibit interannual and seasonal variability (e.g., Raymond and Oh2009). However, the observational coverage varies considerably from one river to another and thus these are described individually below (in order of decreasing freshwater discharge).

Figure 3Long-term changes in the concentration of total alkalinity (TA) and dissolved inorganic carbon (DIC) in the Susquehanna River. (a) Comparison between TA time series from two locations in the river (see Methods section). The interannual variability of the time series is emphasized with a 1-year moving average. The dotted line represents the idealized linear trend used in the control experiment. (b) Comparison between time series of TA and DIC. (c) Comparison between the seasonality of TA and the idealized seasonal cycle used in model simulations.


The Susquehanna River is the river with the most extensive observational record. Two time series of TA that together span a period of 58 years (1960–2017) are compared in Fig. 3a. The blue time series is from Raymond and Oh (2009, site 01540500; see Fig. 1) and the green time series is derived from the United States Geological Survey data (USGS site 01578310; see Najjar et al.2020). The figure shows actual concentrations with no statistical treatment other than a 1-year moving average to emphasize long-term changes. Both time series suggest a long-term increase of ∼9 meq m−3 yr−1 between 1960 and 2017 which has been attributed to decreasing acid inputs following the decline in coal mining activity (Raymond and Oh2009). Given our focus on long-term changes in the bay, we use the linear trend of Fig. 3a in the 2000–2014 experiment and neglect the remaining year-to-year variations visible in the time series. A separate analysis of DIC (which was calculated from TA, pH and temperature using PHREEQC Parkhurst and Appelo1999; see Najjar et al.2020, for details on the analysis) suggests that TA and DIC followed similar trajectories over these decades (Fig. 3b); therefore we also assume linearly increasing DIC in the 2000–2014 experiment (∼10 mmol-C m−3 yr−1). The remaining year-to-year variability apparent in Fig. 3b is considered beyond the scope of the study and not represented in the model experiments. Finally, a seasonal cycle in TA and DIC concentrations is isolated by applying a 1-month moving average over the original time series and then subtracting the 1-year moving average of Fig. 3a, b. This seasonal cycle is fitted to a sinusoid with a 1-year period (Fig. 3c) to represent the low concentrations around March and the high concentrations around September associated with the relative contribution of surface runoff and groundwater (Najjar et al.2020). This idealized seasonality is superimposed on the linear trend described (see Table 2) so that correlations on seasonal timescales between concentrations and freshwater discharge are represented in the model.

Table 2River concentrations of total alkalinitya (TA) and dissolved inorganic carbona (DIC) for the 10 rivers of the model (see Methods section).

a Concentrations are parameterized as DIC(t)=DIC^+aDICcos5π/4-ωt, where DIC^ is a long-term average, ω=2π/(365d), and t is days since year 0 (proleptic calendar). b Value for the year 2007 (the concentrations of the Susquehanna River include a long-term trend during 2000–2014; see Methods section). n/a indicates that no seasonality is prescribed.

Download Print Version | Download XLSX

The Potomac and James rivers have the next highest freshwater discharge. As for the Susquehanna River, USGS data (Najjar et al.2020) are used to parameterize a seasonal cycle for TA and DIC concentrations (Table 2). We do not include a long-term trend for those rivers as the temporal coverage is more limited than for the Susquehanna. A long-term arithmetic mean of the TA and DIC concentrations is calculated instead (based on the years 1975–2005 for the Potomac and 1975–1995 for the James) and superimposed on the seasonal cycle (Table 2). Annual mean concentrations are also calculated from the USGS time series for the Patuxent (based on the years 1985–1999), Rappahannock (1968–1994) and York rivers (1990–1998; Table 2) but no attempt was made at parameterizing their seasonal cycle given their smaller discharge and influence on the bay. The concentration for the York River is calculated from its two tributaries (the Mattaponi and Pamunkey rivers) weighted according to their mean freshwater discharge.

No time series of TA or DIC were available for the remaining four rivers (Elk, Chester, Choptank and Nanticoke rivers) and thus the zero salinity intercept (785 meq m−3) of an alkalinity–salinity relationship derived for the eastern shore of the Chesapeake Bay is used (Najjar et al.2020). This value is assumed constant in time (Table 2). For the DIC concentrations of the same rivers, no relationship is available, and thus we assume that the ratio TA : DIC is 1:1 at a salinity of zero (e.g., Fig. 2 of Friedman et al.2020).

With these assumptions, the bay's mean riverine loading over the period 2000–2014 is 89 Geq yr−1 for TA and 1169 Gg-C yr−1 for DIC (Table 3). The Susquehanna River itself contributes to 45 % of TA, 45 % of DIC and 47 % of freshwater discharge during this period.

Table 3Mean riverine loadings over the two periods of interest (see Methods section)a.

a The values are for the 10 rivers of the model (combined). b The experiments assume the same riverine freshwater discharge in both periods (see Methods).

Download Print Version | Download XLSX

2.1.5 Data available for the evaluation of the control experiment

The hydrodynamic and nitrogen components of the model have been extensively evaluated in previous publications (e.g., Irby et al.2016; Da et al.2018) based on the data from the monitoring program of the Chesapeake Bay (USEPA2012). In order to evaluate the inorganic carbon system of our control experiment, a dataset of TA, DIC and pCO2w is used (Shadwick et al.2019a; Friedman et al.2020). The data were collected along the main stem of the bay (37 to 39.5 N) over the period June 2016 to June 2018 and cover all four seasons. The dataset includes a total of 204 data points of surface TA, DIC and pCO2w. In order to compare these data with the control experiment (years 2000–2014), a seasonal climatology was assembled for the months of December to February, March to May, June to August, and September to November. This combination is chosen so that all seasons include a comparable number of data points.

2.2 Sensitivity experiments

2.2.1 Sensitivity to atmospheric CO2 concentration (1900CO2)

This sensitivity experiment is designed to isolate the impact of atmospheric CO2 from all other drivers of change. It is identical to the control experiment except that it uses atmospheric CO2 concentrations that are representative of the early 1900s. More precisely, the experiment assumes a constant mixing ratio of 300 ppm (representative of the values during 1900–1914 in Miller et al.2014) which is ∼90 ppm lower than during 2000–2014 (Fig. 2). Note that a change in atmospheric CO2 alone does not affect the primary production and respiration in the model, but it does affect the DIC and the carbon budget.

2.2.2 Sensitivity to temperature (1900T)

This sensitivity experiment isolates the impact of temperature change by using water temperatures that are always 1.5 C lower than the control experiment throughout the water column. Few long-term records of water temperature date back to the early 1900s, and thus the 1.5 C value should be viewed as an approximation. The 1.5 C corresponds approximately to the difference in water temperature between the years 1990–2005 (∼16C) and the years 1940–1950 (∼14.5C) at the mouth of the Patuxent River (Najjar et al.2010, their Fig. 3). The uniformity of the change in the vertical is justified by the shallow depths of the bay and supported by the study of Preston (2004). Note also that in this experiment the change in temperature only affects the biogeochemical fields; the uniform change considered here would have only a minor impact on the physical circulation of an estuary like the Chesapeake Bay, and thus we simply use the same physical fields as in the control experiment.

With the temperature-dependent formulations used in the model, the historical increase of 1.5C represents a +11 % increase in the maximum phytoplankton growth rate, maximum grazing rate and remineralization rate. Note that phytoplankton production is also limited by nutrients and light and these two can mitigate the increase expected from temperature alone. Respiration depends on the amount of organic matter present in the water column and is thus expected to mirror changes in production.

2.2.3 Sensitivity to riverine inputs of nitrogen (1900N)

This sensitivity experiment isolates the impact of change in nitrogen loading by using riverine concentrations of nitrogen that are modified to represent conditions of the early 1900s. DLEM riverine concentrations of NO3-, NH4+, PON and DON for the period 1900–1914 are used for this purpose, following the same protocol as for the control experiment. However, the river freshwater discharge remains the same as in the control experiment (2000–2014) so that the physical fields of the model (notably the currents and stratification) are unaffected, and thus the differences between the two runs are solely due to the riverine concentrations of nitrogen. Note that the mean freshwater discharge of the bay's watershed is comparable during the periods 1900–1914 and 2000–2014 (78 and 86 km3 yr−1, respectively, DLEM).

The DLEM results suggest that the main change that occurred between 1900–1914 and 2000–2014 is a large increase in riverine nitrate concentrations (Fig. 4) that occurred primarily in the late 1960s and early 1970s. This increase is well documented (see Harding et al.2016, and the observations in Fig. 4) and is generally attributed to a large increase in nitrogen fertilizer usage after World War II. The bay's DIN loading increased by 33 % between the 1900s and the 2000s (DLEM, Table 3).

Figure 4Long-term changes in the nitrate loading of the Susquehanna River. Observations are from Harding et al. (2016) (their Fig. 7) while model values are from DLEM (see Methods section).


2.2.4 Sensitivity to riverine inputs of carbon and alkalinity (1900C)

This sensitivity experiment uses riverine fluxes of carbon and alkalinity that are modified to represent conditions of 1900–1914. The modifications to TA and DIC are limited to the Susquehanna River as it accounts for approximately half of the freshwater discharge to the bay and it is the only river with >50 years of observations (Fig. 3). For simplicity, we assume a constant annually averaged TA and DIC throughout the period 1900–1914 that is equal to the estimates for the year 1960 from the linear trend (Table 2). The idealized seasonal cycle of TA and DIC remains the same as in the control experiment.

With these assumptions, the bay's riverine loading increases by 25 % (TA) and 22 % (DIC) between the early 1900s and early 2000s (Table 3). Note that such increases in DIC and TA inputs, taken alone, would not affect the primary production and respiration of the model, but they certainly affect the model budgets of inorganic carbon. In the case of total organic carbon (TOC), the bay's riverine loading increases by 11 % from the 1900s to the 2000s (DLEM, Table 3).

2.2.5 Combined effect of atmospheric CO2, temperature, riverine N, C and TA (1900all)

This numerical experiment simultaneously applies all four perturbations described above (atmospheric CO2, temperature, and riverine inputs of nitrogen, carbon and alkalinity), allowing us to test the additivity of the changes caused by the individual perturbations. The sum of the four individual perturbation experiments is unlikely to match experiment 1900all exactly since the perturbations are not acting independently. As one example, the air–sea CO2 flux (F) depends nonlinearly on T, DIC and TA through pCO2w.

Figure 5Overview of the spatial and seasonal variability of the inorganic carbon system: (a) surface dissolved inorganic carbon, (b) surface total alkalinity, (c) partial pressure of CO2 at water surface and (d) air–sea CO2 flux. The shading represents a seasonal climatology from the control experiment (years 2000–2014). The circles are derived from observations (see Methods section). DJF is December to February, MAM is March to May, JJA is June to August and SON is September to November.


2.3 Carbon budgets assembled from the model results

A carbon budget for the Chesapeake Bay (including the tributaries and integrated from surface to bottom) is calculated at every model time step and then averaged over the simulation periods (15 years). The equations of the budget are (e.g., Wakelin et al.2012)


where the terms on the left-hand side of Eqs. (5)–(6) represent changes in DIC and TOC inventory over time. The first two terms on the right-hand side of Eqs. (5)–(6) represent the input from rivers and the net horizontal flux (export) across the bay's mouth (positive seaward). NEP is the net ecosystem production and represents the difference between production and respiration (Eq. 7). Airseaflux is the net air–sea CO2 flux over the domain and it is defined as positive if this term represents a source of carbon to the bay (net ingassing). Burial represents the fraction of the bottom TOC flux that is permanently buried (i.e., not resuspended nor respired).

In the model, the inorganic and organic budgets (Eqs. 56) are closed and there is no residual term. We report budget values rounded to the nearest integer (Gg-C yr−1) as this corresponds to the order of magnitude of the smallest term in the time-averaged budget. Note, however, that the terms have year-to-year variations that far exceed 1 Gg-C yr−1. This variance is quantified by the standard deviation of the annually averaged budget terms and is indicated with the symbol ±. The standard deviation is rounded to the nearest 10 Gg-C yr−1 in the text and in the tables.

3 Results

3.1 Control experiment (2000–2014)

3.1.1 Overview of the inorganic carbon system

The inorganic carbon system of the model exhibits important regional differences (Fig. 5). The tributaries of the bay, including the major Susquehanna River in the north (Fig. 1), are associated with relatively low DIC and TA, and they produce a gradient increasing seaward along the main stem of the bay (Fig. 5a, b; see also Brodeur et al.2019, and Friedman et al.2020). The lower DIC and TA concentrations are particularly apparent in the northern half of the bay (Fig. 5a, b) where the Susquehanna River delivers ∼45 % of the freshwater discharge to the bay. One exception to the low tributary DIC and TA concentrations is the Potomac River where concentrations are higher than all other rivers (Table 2) and approach those of shelf water in the fall season.

Figure 6Seasonal climatology of the net primary production (NPP) over the period 2002–2011. The blue curves are the modeled results and the green curves are from the empirical model of Son et al. (2014) (their Fig. 7). The upper, middle and lower bay regions are defined as in Son et al. (2014) (their Fig. 1). The “+” symbols represent the annual mean value of the curves. The model results are from the control experiment (Table 1).


The low TA of the river water is accompanied by relatively high surface pCO2w inside the tributaries and downstream of the Susquehanna River (Fig. 5c). Away from the tributaries, surface pCO2w values are generally close to atmospheric levels (∼385µatm in 2000–2014). This spatial distribution of surface pCO2w drives strong outgassing within the tributaries and in the northern half of the bay and either ingassing or near-neutral conditions in the southern half (Fig. 5d).

The inorganic carbon system also exhibits large seasonal variations. The signature of seasonal river inputs is most apparent by comparing the period March–May (after the spring freshet results in low surface salinity) to the period September–November (after the low freshwater inputs of the summer season). In March–May, the fresh riverine water contributes to low DIC and TA concentrations in the estuary whereas September–November show higher DIC and TA (Fig. 5a, b). Similarly, pCO2w oscillates seasonally with lower values in March–May and the highest values in September–November (Fig. 5c). Biological production and temperature contribute to this seasonality of pCO2w, with warming temperatures increasing pCO2w and increasing production mitigating this (e.g., Friedman et al.2020). The seasonality of pCO2w is reflected in the air–sea CO2 fluxes and results in strong outgassing or near-neutral conditions in September–November, while March–May are characterized by weaker outgassing (in the northern bay) and even strong ingassing in the southern bay (Fig. 5d). In the period June–August, pCO2w is relatively low (close to atmospheric concentrations), resulting in near-neutral air–sea fluxes (Fig. 5c, d).

3.1.2 Evaluation of the modeled inorganic carbon system

As noted earlier, the hydrodynamics and nitrogen cycle of the linked DLEM–ChesROMS–ECB modeling system have been well evaluated (Feng et al.2015; Irby et al.2016; see also the Supplement for an assessment of the main model variables in Table S1, Figs. S1–S2). Therefore, here the model evaluation is focused on the carbon cycling component of the model (Fig. 5). Regional and seasonal variabilities highlighted in the previous section are generally supported by the observations of Friedman et al. (2020) (see also Brodeur et al.2019). These include the north–south gradient in properties, the seasonally varying influence of the rivers and the seasonality of pCO2w away from the tributaries. Some discrepancies are apparent, however. At the mouth of the bay, the model generally underestimates surface DIC concentrations (Fig. 5a; mean DIC bias is −215 mmol-C m−3 at the mouth of the bay). This bias is particularly apparent in March–May and extends into the southern half of the bay (Fig. 5a; bay-averaged DIC bias in March–May is −155 mmol-C m−3).

The bias in DIC directly affects surface pCO2w, which is similarly biased low in the southern bay in March–May (Fig. 5c; bay-averaged bias in March–May is −130µatm). We note, however, that this bias is less and less apparent away from the bay's mouth and that the modeled pCO2w in the main stem bay agrees well with the data points in the vicinity of the Potomac River (Fig. 5c). In the northern half of the bay, observed pCO2w sometimes exhibits noisy patterns (particularly in September–November; Fig. 5c) that are not reproduced by the model. The potential causes of these differences will be discussed later (see Discussion section).

The spatiotemporal variability of the inorganic carbon system can be strongly influenced by biological production within the bay. This important component of the model is evaluated for the period 2002–2011 using results from an empirical satellite productivity model calibrated with in situ observations (Son et al.2014). The empirical model provides a seasonal climatology of net primary production (NPP) for three subregions of the bay's main stem (Fig. 6). Results from both models exhibit a strong seasonal cycle with peak NPP between the months of May and July (consistent with in situ data in Fig. 4 of Harding et al.2002) and a similar magnitude of NPP during the summer. They also agree on the differences between the three regions, with summer NPP being highest in the upper bay and lowest in the lower bay. The primary difference between the two sets of model results is that ChesROMS–ECB generates higher production in the winter months in the lower bay (Fig. 6).

Table 4Carbon budget termsa (Gg-C yr−1) for the control experiment (2000–2014) and deviationd of the sensitivity experiments from the control experiment (see Table 1 and Eqs. 57)

a The values are averaged over the period of the simulation and rounded to the nearest integer. b The symbol ± indicates the interannual variability (see Methods). c “Sum” is the sum of the deviations associated with experiments 1900CO2, 1900T, 1900N and 1900C. d A version of this table with absolute values (rather than deviations) is available in the Supplement.

Download Print Version | Download XLSX

3.1.3 Combined inorganic and organic carbon budget

A carbon budget (Eqs. 56) for the Chesapeake Bay domain (including the tributaries but excluding the continental shelf) is calculated over the simulation period of the control experiment (years 2000–2014, Table 4). If we first consider the total carbon (the sum of inorganic and organic carbon), the budget shows a near balance between riverine carbon inputs and advective output at the mouth of the bay (“Export”; see Table 4). The difference between the two (“Riv-exp.”, 196 Gg-C yr−1, Table 4) is equivalent to ∼12 % of the annual riverine carbon input and is largely balanced by burial within the bay (221±20 Gg-C yr−1). In comparison with these terms, the carbon inventory of the bay shows a very small positive trend over the 15 years of the simulation but substantial year-to-year variability (+8±60 Gg-C yr−1, Table 4).

Figure 7Summary of the changes in the inorganic carbon system for the six model experiments (Table 1). “Rivers minus export” combines the riverine DIC input and the export of DIC at the bay's mouth (Table 4). NEP is net ecological production (see Methods section).


Despite the near balance between the riverine carbon input and the carbon export at the bay's mouth, considerable biogeochemical transformations take place within the budget domain. Production and respiration are each equivalent to ∼270 % of the annual riverine carbon input (Table 4). The difference between production and respiration (NEP) is, however, an order of magnitude smaller (+259±60 Gg-C yr−1, positive indicating net autotrophy; Fig. 7b). NEP is thus equivalent to ∼15 % of the riverine carbon input and is comparable in magnitude to burial (Table 4). The standard deviation of NEP is relatively small (<20 % of the standard deviation associated with production or respiration) as years of high production are also years of high respiration (not shown).

As noted in Fig. 5, the air–sea CO2 flux exhibits ingassing or near-neutral fluxes in the southern half of the bay and strong outgassing within the tributaries and in the northern half of the bay (i.e., downstream of the Susquehanna River). The net air–sea CO2 flux of the bay, defined as ingassing minus outgassing, is very close to zero (+34 Gg-C yr−1, i.e., slightly ingassing; Table 4 and Fig. 7b). The sign of the net flux is thus sensitive to environmental changes and fluctuates substantially from one year to another (negative for 4 years and positive for 11 years). The standard deviation over 2000–2014 of the net air–sea flux is ±90 Gg-C yr−1.

3.2 Sensitivity experiment results: changes in the carbon budget

3.2.1 Experiment 1900CO2 versus control experiment

The increase in atmospheric CO2 concentrations (experiment 1900CO2 versus control experiment) only affects the inorganic component of the carbon budget (Table 4, Fig. 7c). The production, respiration, burial and export of organic carbon in experiment 1900CO2 are thus identical to the control experiment. The historical change in atmospheric CO2 is large enough to reverse the sign of the net air–sea flux from -20±90 Gg-C yr−1 (slightly outgassing in 1900CO2) to +34±90 Gg-C yr−1 (slightly ingassing in the control experiment; Table 4). The increase in the net air–sea CO2 flux is accompanied by an increase in DIC export of similar magnitude (+51 Gg-C yr−1, Table 4, Fig. 7c). Note that this increase in export reflects higher DIC concentrations within the bay and not a change in the physical circulation of the bay. As in the control experiment, the trends in inorganic and organic carbon are very small over the 15 years of the experiment (Table 4).

3.2.2 Experiment 1900T versus control experiment

The increase in water temperature (experiment 1900T versus control experiment) mostly affects the production and respiration, with increases of +252 and +265 Gg-C yr−1, respectively (Table 4) and only a small resulting change in NEP (−13 Gg-C yr−1; Fig. 7d). The net air–sea CO2 flux over the domain changes from +57 to +34 Gg-C yr−1; i.e., the change in temperature brings the bay closer to being neutral (Fig. 7d). Note that the increase in temperature affects surface pCO2w and contributes to the change in air–sea CO2 flux and to a 13 Gg-C yr−1 reduction in the DIC export (Fig. 7d, Table 4). The other components of the budget (burial and export of organic carbon) are mostly unchanged by the warming.

3.2.3 Experiment 1900N versus control experiment

The increase in riverine inputs of nitrogen (experiment 1900N versus control experiment) has a strong impact on production and respiration (Table 4). These two terms are increased by +492 and +391 Gg-C yr−1 (respectively), resulting in a NEP increase of +101 Gg-C yr−1 (Fig. 7e, Table 4). This change in NEP affects the organic component of the budget, increasing both burial and TOC export at the bay's mouth by similar amounts (+55 and +46 Gg-C yr−1, respectively). The net air–sea CO2 flux shows the largest change of all the sensitivity experiments, changing from −36 (slightly outgassing) to +34 Gg-C yr−1 (slightly ingassing). This change is consistent with the increase in NEP and with a decrease in surface DIC concentrations. Finally, the increase in riverine inputs of nitrogen produces a decrease in DIC export of 31 Gg-C yr−1 (Fig. 7e, Table 4).

3.2.4 Experiment 1900C versus control experiment

The increase in riverine inputs of carbon and alkalinity (experiment 1900C versus control experiment) leads to an increase in the respiration term (+30 Gg-C yr−1, Table 4). The latter is solely a result of the increased riverine TOC loading as TA and DIC alone would not affect respiration. Since the production is nearly unchanged, overall the NEP exhibits a decrease of 28 Gg-C yr−1 (i.e., the bay is becoming less autotrophic). The net air–sea CO2 flux decreases from +76 to +34 Gg-C yr−1, meaning that the change in the riverine carbon and TA brings the bay closer to being neutral (Table 4, Fig. 7f). The change in TA and DIC contributes to this change in air–sea CO2 flux by their impact on surface pCO2w. Assuming an annually averaged water temperature of 15C, conservative mixing between the properties of the Susquehanna River (Table 2) and an oceanic end-member defined by S=33 psu and Eqs. (3)–(4), we estimate an increase of up to ∼9 % in surface pCO2w between the 1900s and the 2000s. Finally, the increase in riverine inputs of DIC is accompanied by a similar increase in DIC export, leading to a small net effect on the horizontal DIC transport (+19 Gg-C yr−1; Table 4, Fig. 7f).

3.2.5 Experiment 1900all versus control experiment

When all four changes are made simultaneously (i.e., experiment 1900all with simultaneously increased atmospheric CO2, temperature, nitrogen loading, DIC and TA loadings; Fig. 7a), the results differ from the control experiment primarily in terms of air–sea CO2 flux and NEP. The air–sea CO2 flux switches from a small net source to the atmosphere (−8 Gg-C yr−1) to a net sink (+34 Gg-C yr−1), and NEP becomes increasingly autotrophic (209 to 259 Gg-C yr−1; Fig. 7a). The results generated in experiment 1900all are also similar to what is obtained by adding the results of the four sensitivity experiments described above (i.e., compare the last two columns in Table 4), suggesting a substantial linearity between these four experiments. Some differences do exist, however. Specifically, when the four experiments are run simultaneously, the changes in the net air–sea CO2 flux, burial and NEP are all slightly smaller than when the results of the four experiments are simply added together.

3.3 Relative importance of global and regional changes

The relative importance of global and regional changes is assessed by combining the experiments as follows: 1900CO2+1900T (global) and 1900N+1900C (regional). An important result is that this grouping narrows the gap between the early 1900s and early 2000s by combining changes of opposite signs (Fig. 7c–f). For example, the change in air–sea CO2 flux from rising atmospheric CO2 concentrations is partially mitigated by the increased outgassing from rising temperatures. The net effect of global changes is to decrease the net horizontal flux of DIC (“Rivers – export”) by 38 Gg-C yr−1, increase the net air–sea CO2 flux by 31 Gg-C yr−1 and decrease the NEP by 12 Gg-C yr−1 (Fig. 7, Table 4).

In the case of regional drivers, the large increases in ingassing and NEP from increased DIN loadings are partially mitigated by the changes in riverine inputs of carbon and alkalinity (Fig. 7e–f). The net effect of regional changes is to increase the net horizontal flux of DIC (Rivers – export) by 50 Gg-C yr−1, increase the net air–sea CO2 flux by 27 Gg-C yr−1 and increase the NEP by 73 Gg-C yr−1 (Table 4). Global and regional drivers are thus of similar importance when assessing changes in the inorganic carbon balance, with the exception of NEP which is primarily affected by regional drivers. Note that global and regional drivers both push the bay toward net ingassing but they influence the horizontal DIC flux and NEP in opposite ways.

4 Discussion

4.1 Uncertainties and comparison with other studies

The inorganic carbon budget of the control experiment (Table 4) can be compared to that of Shen et al. (2019b). The two model-derived budgets share the same key features, namely, a positive net horizontal flux (Rivers – export) of DIC (234 and 157 Gg-C yr−1, respectively), a positive NEP (259 and 165 Gg-C yr−1, respectively) and a comparatively small net air–sea CO2 flux (34 and −50 Gg-C yr−1, respectively). The main discrepancy is the riverine DIC loading (1169 in this study and 821 Gg-C yr−1 in Shen et al.2019b). The cause of this difference is unclear as the two studies assume similar DIC concentrations for the largest river (Susquehanna River, not shown). Potential explanations include differences in the years examined (2000–2014 versus 1986–2015), differences in the riverine freshwater discharge used in the simulations or differences in the DIC concentrations assumed for the smaller rivers. Because most of the DIC loading from the rivers is exported to the coastal ocean, these differences are unlikely to cause major discrepancies in the 1900s versus 2000s changes reported in this study.

The budget of the control experiment (Table 4) can also be compared to that of Kemp et al. (1997). They estimate riverine loadings of 55.8 Gg-N yr−1 (DIN), 39.5 Gg-N yr−1 (TON, total organic nitrogen) and 261.8 Gg-C yr−1 (TOC). While their TON loading is similar to that in our budget, their DIN and TOC loadings are 42 % and 48 % lower (respectively) than in our budget. The difference in DIN loading must originate from the smaller tributaries of the bay since the DIN loading of the Susquehanna is nearly identical between the two studies. In the case of the TOC loading, the discrepancy remains whether we focus on the Susquehanna or the watershed. The TOC export is also 48 % lower in Kemp et al. (1997) than in the present study (with the caveat that the two budgets represent different years). The other components of the budget in Kemp et al. (1997) are not directly comparable to our study as they are specific to the main stem of the bay.

Although in general the model results represent recent climate-quality data in the Chesapeake Bay quite well (Fig. 5), it is worth discussing the origin and impact of the model biases and how these may be reduced in future work. For example, the low DIC bias at the mouth of the bay (Fig. 5a) most likely originates from uncertainties in the DIC concentrations prescribed at the model's oceanic boundary, which are derived from limited measurements (Methods section). The low DIC bias leads, in turn, to a low bias in pCO2w in March–August in the southern half of the bay (Fig. 5c). The observed pCO2w values suggest a relatively weak outgassing in this region and time of the year, while the model exhibits a weak ingassing (Fig. 5d). The bias is, however, unlikely to have a major impact on net air–sea CO2 flux of the model as it appears to be geographically confined to the southern part of the bay (Fig. 5c). In future implementations of the model, more climate-quality data will be used to improve this outer boundary condition issue.

Differences in modeled and observed pCO2w are also apparent in the northern half of the bay (Fig. 5c). A possible explanation for these differences is a temporal mismatch between the observations and the model results (which are from different years; see Methods). Such a mismatch in years can cause substantial differences in the water properties of this area as the freshwater discharge of the Susquehanna River varies substantially between years and often controls the alongshore gradients (Zhang et al.2006).

Historical changes that were not considered in the present study include alkalinity sinks within tributaries such as the Potomac River (Najjar et al.2020) due to biogeochemical processes not accounted for here. Other historical changes not considered in the present study include the warming, salification and acidification of continental shelf waters (Wallace et al.2018; Saba et al.2015). Future studies should consider the role that these oceanic changes have played over the past century. Finally, it is worth pointing out that the important topic of coastal acidification (Cai et al.2011, 2017) was not examined in this study but that it should be a focal point of future studies.

4.2 Changes in Chesapeake Bay carbonate chemistry over the past century

There have been considerable changes to the inorganic carbonate system of the Chesapeake Bay over the past century (Table 4). Causes include both global factors, including increases in atmospheric CO2 and increases in temperature, and more regional factors within the watershed, including increases in nitrogen and alkalinity loadings. The results from this study demonstrate that together these changes have only slightly altered the net advective flux of DIC into the bay: the difference between DIC river inputs and export to the coastal ocean has changed by only 6 % over the past century (Fig. 7a, b; Table 4). In contrast the changes in NEP and air–sea flux have been considerably larger. The bay has become 19 % more autotrophic over this time period (Fig. 7a, b; Table 4), and the bay's net air–sea flux has switched from being a small net source of CO2 to the atmosphere, to a sink of CO2 from the atmosphere. In the sections below, the causes of these overall changes, identified via the sensitivity experiments described above, are discussed individually, including both global changes (atmospheric CO2 and temperature) and regional watershed changes (riverine nitrogen, carbon and TA). In each case, there are mitigating factors that cause the changes to be lower than otherwise expected.

4.2.1 Global changes and their impact on Chesapeake Bay carbonate chemistry

Between the early 1900s (1900–1914) and the early 2000s (2000–2014) atmospheric CO2 concentration has increased by roughly 100 ppm (Etheridge et al.1996; Keeling et al.2003). As expected, the impact of this single change on the inorganic carbon budget of the Chesapeake Bay is significant, resulting in the transformation of the bay from an average net source of CO2 to the atmosphere (outgassing: -20±80 Gg-C yr−1) to a net sink (ingassing: +34±90 Gg-C yr−1). It is important to note that the standard deviations associated with these interannual means in air–sea CO2 flux represent interannual variability and are significantly larger than the estimated long-term change. Thus, although the increase in atmospheric CO2 is clearly increasing ingassing on average, there are still large year-to-year differences that may cause certain years in the early 1900s to be net sinks of atmospheric CO2 and certain years in the early 2000s to be net sources of atmospheric CO2. This interannual variability makes it difficult to determine the average direction of the net air–sea CO2 flux over the estuary unless long time series of climate-quality observations are available.

In addition to increases in atmospheric CO2, atmospheric and estuarine temperatures have also been rising over the past century (Ding and Elmore2015; Muhling et al.2018; Irby et al.2018). The increased ingassing due to elevated atmospheric CO2 is partially mitigated (by roughly 50 %; Fig. 7c, d) via these increasing temperatures, which enhance pCO2w (because of solubility but also more respiration; Fig. 7d). As a result, the change in air–sea CO2 flux due to the global changes of the past century (+31 Gg-C yr−1) is only half as large as it would be without the concomitant increase in temperature. The increase in water temperature also leads to a 5 % decrease in net ecosystem production through enhanced respiration of organic matter, consistent with heterotrophic processes being more sensitive to temperature than production (Lomas et al.2002).

4.2.2 Regional watershed changes and their impacts on Chesapeake Bay carbonate chemistry

The increase in riverine DIN loading associated with urbanization and increased fertilizer usage has caused changes in the inorganic carbon budget over the last century that are nearly equal to those induced by global changes (Fig. 7). Specifically, increased nitrogen loading has caused NEP to increase substantially (Fig. 7e, +39 %), and this, in turn, leads to lower surface DIC (-10 mmol-C m−3) and surface pCO2w (-25µatm) in the southern half of the bay. In response to these changes, the net air–sea CO2 flux into the bay has increased considerably (Fig. 7e, +70 Gg-C yr−1), which is an even larger increase than that due to the higher atmospheric CO2 (+54 Gg-C yr−1). Another consequence of the enhanced NEP and lower DIC concentrations is a reduction in DIC export to the shelf (Fig. 7e).

The historical increase in riverine carbon (TOC, DIC) and TA loadings has had relatively minor impacts on the inorganic carbon budget, compared to those due to increased nitrogen loading discussed above (compare Fig. 7e and f). Although significantly increased DIC loading to the bay is assumed (Table 2), and although the DIC concentrations of the bay are increased substantially, much of the extra riverine DIC is simply exported to the coastal ocean (94 %, Table 4). In terms of TOC, only 38 % of the increase is exported to the coastal ocean. The remaining increase in TOC serves to increase respiration (decrease in NEP), partially offsetting the increase in production that resulted from the increased nitrogen loading discussed above. Regarding the air–sea CO2 flux, the net effect of the increased respiration and increased riverine DIC and TA loadings is a relatively small increase in pCO2w (approximately +6 % on average over seasons and over the bay between experiment 1900C and the control experiment) that brings the net air–sea CO2 closer to being neutral. This ultimately serves to largely counteract the increased ingassing resulting from the increased nitrogen loading. Thus when considered together, the increases in nitrogen and carbon loading over the past century have resulted in the Chesapeake Bay becoming a greater sink for atmospheric CO2 (by 27 Gg-C yr−1), which is similar in magnitude to the increased sink due to global changes (+31 Gg-C yr−1).

5 Summary and concluding remarks

Sensitivity experiments were performed to isolate the effect of changes in (1) atmospheric CO2, (2) temperature, (3) riverine nitrogen loading, and (4) riverine carbon and alkalinity loading, on the inorganic carbon balance of the Chesapeake Bay between the early 1900s and early 2000s. Limited information is available for the early 1900s and thus these experiments are meant to highlight the aforementioned changes rather than to model actual early 1900s conditions. Both regional and global changes have enhanced the bay's sink for atmospheric CO2 by similar amounts. The increased riverine nitrogen load, a regional change, increased production, which resulted in the bay having a 19 % higher (more autotrophic) NEP. Overall, the results of the study help clarify the impact that local management efforts (past or future) can have on the bay's inorganic carbon balance and the limits of these efforts in the context of ongoing global changes. The temporal and spatial scope of this study also highlights the usefulness of modeling studies and how difficult it is to answer questions on these spatial and temporal scales from observations alone.

The comparison between the early 1900s and early 2000s suggests that the ongoing increase in atmospheric CO2 concentrations overshadows the temperature-driven increase in pCO2w and outgassing. In other words, the bay's trend toward more uptake of atmospheric CO2 will likely continue in the decades to come. This is in contrast with regional changes in riverine loadings (mostly DIN, TA and DIC) which were particularly large in the past century and are not expected to continue in the future. Management efforts in the bay's watershed, notably the implementation of a total maximum daily load (USEPA2010; Irby and Friedrichs2019), are expected to stabilize or reduce the nutrient inputs to the bay over the next several decades. Similarly, Raymond and Oh (2009) (their Fig. 1) suggest that coal production has stabilized in the past decades, and thus one would expect the Susquehanna's alkalinity and DIC concentrations to also stabilize. Overall, these results suggest that although changes in riverine nutrient inputs have played an important role in altering coastal carbon budgets over the past century, in the future ongoing global changes may have an even greater affect on coastal carbonate chemistry.

Code and data availability

Documentation and source code for the numerical model used in this study (ROMS) are publically available at (last access: 23 February 2018) (ROMS/TOMS Group2018). Additional documentation about the biogeochemical equations and parameters is provided in the Supplement of this publication. The model results used in the paper are permanently archived on an online public repository: (St-Laurent and Friedrichs2020).


The supplement related to this article is available online at:

Author contributions

EHS, MAMF and RGN acquired the financial support and supervised the project leading to this publication. YY and HT contributed model inputs and forcings necessary for the model simulations. MAMF and PS designed the numerical experiments and PS carried them out. PS analyzed the model results and created the figures and tables while all co-authors contributed to their interpretation. PS prepared the manuscript with contributions from all co-authors. PS was responsible for data curation of the model outputs.

Competing interests

The authors declare that they have no conflict of interest.


We thank the referees for their careful reading of the manuscript and for providing helpful and thoughtful comments. We acknowledge the contribution of Edward G. Stets, who generously provided us with time series of alkalinity and DIC for the main tributaries of the Chesapeake Bay and provided insightful comments on the manuscript. This paper is the result of research funded by the National Oceanic and Atmospheric Administration's Ocean Acidification Program under award NA18OAR0170430 to the Virginia Institute of Marine Science. This research was also supported by the National Science Foundation (collaborative grants OCE-1537013, 1536996), NASA IDS (NNX14AF93G) and NASA Carbon Cycle Science to WETCARB (NNX14AM37G). This work used the Extreme Science and Engineering Discovery Environment (XSEDE) supercomputer Comet at SDSC through allocation OCE-160013, which is supported by National Science Foundation grant number ACI-1548562. The authors also acknowledge William & Mary Research Computing (, last access: July 2020) for providing computational resources and/or technical support that have contributed to the results reported within this paper. This is contribution no. 3918 of the Virginia Institute of Marine Science, William & Mary.

Financial support

This research has been supported by the National Science Foundation, Division of Ocean Sciences (grant nos. 1537013 and 1536996), the National Aeronautics and Space Administration (grant nos. NNX14AF93G and NNX14AM37G) and the National Oceanic and Atmospheric Administration's Ocean Acidification Program (grant no. NA18OAR0170430).

Review statement

This paper was edited by Stefano Ciavatta and reviewed by two anonymous referees.


Brodeur, J. R., Chen, B., Su, J., Xu, Y. Y., Hussain, N., Scaboo, K. M., Zhang, Y., Testa, J. M., and Cai, W. J.: Chesapeake Bay inorganic carbon: Spatial distribution and seasonal variability, Front. Mar. Sci., 6, 99,, 2019. a, b, c

Cai, W. J.: Estuarine and coastal ocean carbon paradox: CO2 sinks or sites of terrestrial carbon incineration?, Annu. Rev. Mar. Sci., 3, 123–145,, 2011. a

Cai, W. J., Hu, X., Huang, W. J., Murrell, M. C., Lehrter, J. C., Lohrenz, S. E., Chou, W. C., Zhai, W., Hollibaugh, J. T., Wang, Y., Zhao, P., Guo, X., Gundersen, K., Dai, M., and Gong, G. C.: Acidification of subsurface coastal waters enhanced by eutrophication, Nat. Geosci., 4, 766–770,, 2011. a

Cai, W. J., Huang, W. J., Luther III, G. W., Pierrot, D., Li, M., Testa, J., Xue, M., Joesoef, A., Mann, R., Brodeur, J., Xu, Y. Y., Chen, B., Hussain, N., Waldbusser, G. G., Cornwell, J., and Kemp, W. M.: Redox reactions and weak buffering capacity lead to acidification in the Chesapeake Bay, Nat. Commun., 8, 1–12,, 2017. a, b

Cerco, C. F. and Noel, M. R.: Process-based primary production modeling in Chesapeake Bay, Mar. Ecol. Prog. Ser., 282, 45–58, 2004. a

Da, F., Friedrichs, M. A. M., and St-Laurent, P.: Impacts of atmospheric nitrogen deposition and coastal nitrogen fluxes on oxygen concentrations in Chesapeake Bay, J. Geophys. Res.-Oceans, 123, 5004–5025,, 2018. a, b, c, d, e, f, g, h, i, j

Ding, H. and Elmore, A. J.: Spatio-temporal patterns in water surface temperature from Landsat time series data in the Chesapeake Bay, U.S.A., Remote Sens. Environ., 168, 335–348,, 2015. a

Druon, J., Mannino, A., Signorini, S., McClain, C., Friedrichs, M. A. M., Wilkin, J., and Fennel, K.: Modeling the dynamics and export of dissolved organic matter in the Northeastern US continental shelf, Estuar. Coast. Shelf Sci., 88, 488–507,, 2010. a

Duce, R. A., LaRoche, J., Altieri, K., Arrigo, K. R., Baker, A. R., Capone, D. G., Cornell, S., Dentener, F., Galloway, J., Ganeshram, R. S., Geider, R. J., Jickells, T., Kuypers, M. M., Langlois, R., Liss, P. S., Liu, S. M., Middelburg, J. J., Moore, C. M., Nickovic, S., Oschlies, A., Pedersen, T., Prospero, J., Schlitzer, R., Seitzinger, S., Sorensen, L. L., Uematsu, M., Ulloa, O., Voss, M., Ward, B., and Zamora, L.: Impacts of Atmospheric Anthropogenic Nitrogen on the Open Ocean, Science, 320, 893–897,, 2008. a

Etheridge, D., Steele, L., Langenfelds, R., Francey, R., Barnola, J., and Morgan, V.: Natural and anthropogenic changes in atmospheric CO2 over the last 1000 years from air in Antarctic ice and firn, J. Geophys. Res., 101, 4115–4128,, 1996. a

Feng, Y., Friedrichs, M. A. M., Wilkin, J., Tian, H., Yang, Q., Hofmann, E. E., Wiggert, J. D., and Hood, R. R.: Chesapeake Bay nitrogen fluxes derived from a land-estuarine ocean biogeochemical modeling system: Model description, evaluation, and nitrogen budgets, J. Geophys. Res.-Biogeo., 120, 1666–1695,, 2015. a, b, c, d, e, f

Fennel, K., Wilkin, J., Previdi, M., and Najjar, R.: Denitrification effects on air-sea CO2 flux in the coastal ocean: Simulations for the northwest North Atlantic, Geophys. Res. Lett., 35, L24608,, 2008. a

Filippino, K. C., Bernhardt, P. W., and Mulholland, M. R.: Chesapeake Bay Plume Morphology and the Effects on Nutrient Dynamics and Primary Productivity in the Coastal Zone, Estuar. Coast., 32, 410–424,, 2009. a

Filippino, K. C., Mulholland, M. R., and Bernhardt, P. W.: Nitrogen uptake and primary productivity rates in the Mid-Atlantic Bight (MAB), Estuar. Coast. Shelf Sci., 91, 13–23,, 2011. a

Friedman, J. R., Shadwick, E. H., Friedrichs, M. A. M., Najjar, R. G., De Meo, O. A., Da, F., and Smith, J. L.: Seasonal variability of the CO2 system in a large coastal plain estuary, J. Geophys. Res.-Oceans, 125, e2019JC015609,, 2020. a, b, c, d, e, f

Ganju, N. K., Brush, M. J., Rashleigh, B., Aretxabaleta, A. L., Del Barrio, P., Grear, J. S., Harris, L. A., Lake, S. J., McCardell, G., and O'Donnell, J.: Progress and challenges in coupled hydrodynamic-ecological estuarine modeling, Estuar. Coast., 39, 311–332,, 2016. a

Harding, L. W., Mallonee, M. E., and Perry, E. S.: Toward a predictive understanding of primary productivity in a temperate, partially stratified estuary, Estuarine, Coast. Shelf Sci., 55, 437–463,, 2002. a

Harding, L. W., Gallegos, C. L., Perry, E. S., Miller, W. D., Adolf, J. E., Mallonee, M. E., and Paerl, H. W.: Long-Term Trends of Nutrients and Phytoplankton in Chesapeake Bay, Estuar. Coast., 39, 664–681,, 2016. a, b

Irby, I. D. and Friedrichs, M. A. M.: Evaluating Confidence in the Impact of Regulatory Nutrient Reduction on Chesapeake Bay Water Quality, Estuar. Coast., 42, 16–32,, 2019. a

Irby, I. D., Friedrichs, M. A. M., Friedrichs, C. T., Bever, A. J., Hood, R. R., Lanerolle, L. W. J., Li, M., Linker, L., Scully, M. E., Sellner, K., Shen, J., Testa, J., Wang, H., Wang, P., and Xia, M.: Challenges associated with modeling low-oxygen waters in Chesapeake Bay: a multiple model comparison, Biogeosciences, 13, 2011–2028,, 2016. a, b

Irby, I. D., Friedrichs, M. A. M., Da, F., and Hinson, K. E.: The competing impacts of climate change and nutrient reductions on dissolved oxygen in Chesapeake Bay, Biogeosciences, 15, 2649–2668,, 2018. a, b

Keeling, C. D., Piper, S. C., Whorf, T. P., and Keeling, R. F.: Evolution of natural and anthropogenic fluxes of atmospheric CO2 from 1957 to 2003, Tellus B, 63, 1–22,, 2003. a, b

Kemp, W. M., Smith, E. M., Marvin-DiPasquale, M. M., and Boynton, W. R.: Organic carbon balance and net ecosystem metabolism in Chesapeake Bay, Mar. Ecol. Prog. Ser., 150, 229–248, 1997. a, b, c

Lomas, M. W., Gilbert, P. M., Shiah, F. K., and Smith, E. M.: Microbial processes and temperature in Chesapeake Bay: Current relationships and potential impacts of regional warming, Glob. Change Biol., 8, 51–70,, 2002. a, b, c

Mehrbach, C., Culberson, C. H., Hawley, J. E., and Pytkowicz, R. M.: Measurement of the apparent dissociation constants of carbonic acid in seawater at atmospheric pressure, Limnol. Oceanogr., 18, 897–907, 1973. a

Mesinger, F., DiMego, G., Kalnay, E., Mitchell, K., Shafran, P. C., Ebisuzaki, W., Jovic, D., Woollen, J., Rogers, E., Berbery, E. H., Ek, M. B., Fan, Y., Grumbine, R., Higgins, W., Li, H., Lin, Y., Manikin, G., Parrish, D., and Shi, W.: North American Regional Reanalysis: A long-term, consistent, high-resolution climate dataset for the North American domain, as a major improvement upon the earlier global reanalysis datasets in both resolution and accuracy, B. Am. Meteorol. Soc., 87, 343–360, 2006. a

Miller, R. L., Schmidt, G. A., Nazarenko, L. S., Tausnev, N., Bauer, S. E., Genio, A. D. D., Kelley, M., Lo, K. K., Ruedy, R., Shindell, D. T., Aleinov, I., Bauer, M., Bleck, R., Canuto, V., Chen, Y. H., Cheng, Y., Clune, T. L., Faluvegi, G., Hansen, J. E., Healy, R. J., Kiang, N. Y., Koch, D., Lacis, A. A., LeGrande, A. N., Lerner, J., Menon, S., Oinas, V., PerezGarcia-Pando, C., Perlwitz, J. P., Puma, M. J., Rind, D., Romanou, A., Russell, G. L., Sato, M., Sun, S., Tsigaridis, K., Unger, N., Voulgarakis, A., Yao, M. S., and Zhang, J.: CMIP5 historical simulations (1850-2012) with GISS ModelE2, J. Adv. Model. Earth Sy., 6, 441–477,, 2014. a, b, c

Millero, F. J.: Thermodynamics of the carbon dioxide system in the oceans, Geochim. Cosmochim. Ac., 59, 661–677, 1995. a

Muhling, B. A., Gaitan, C. F., Stock, C. A., Saba, V. S., Tommasi, D., and Dixon, K. W.: Potential Salinity and Temperature Futures for the Chesapeake Bay Using a Statistical Downscaling Spatial Disaggregation Framework, Estuar. Coast., 41, 349–372,, 2018. a

Najjar, R. G., Pyke, C. R., Adams, M. B., Breitburg, D., Hershner, C., Kemp, M., Howarth, R., Mulholland, M. R., Paolisso, M., Secor, D., Sellner, K., Wardrop, D., and Wood, R.: Potential climate-change impacts on the Chesapeake Bay, Estuarine, Coast. Shelf Sci., 86, 1–20,, 2010. a

Najjar, R. G., Herrmann, M., Cintron Del Valle, S., Friedman, J. R., Friedrichs, M. A. M., Harris, L. A., Shadwick, E. H., Stets, E. G., and Woodland, R. J.: Alkalinity in tidal tributaries of the Chesapeake Bay, J. Geophys. Res.-Oceans, 125, e2019JC015597,, 2020. a, b, c, d, e, f

Oeschger, H., Siegenthaler, U., Schotterer, U., and Gugelmann, A.: A box diffusion model to study the carbon dioxide exchange in nature, Tellus, 27, 168–192,, 1975. a

Parkhurst, D. L. and Appelo, C. A. J.: User's guide to PHREEQC–A computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations, Tech. rep., U.S. Geological Survey, 1999. a

Preston, B. L.: Observed Winter Warming of the Chesapeake Bay Estuary (1949–2002): Implications for Ecosystem Management, Environ. Manage., 34, 125–139,, 2004. a

Raymond, P. A. and Oh, N. H.: Long term changes of chemical weathering products in rivers heavily impacted from acid mine drainage: Insights on the impact of coal mining on regional and global carbon and sulfur budgets, Earth Planet. Sc. Lett., 284, 50–56,, 2009. a, b, c, d

Riebesell, U., Zondervan, I., Rost, B., Tortell, P. D., Zeebe, R. E., and Morel, F. M.: Reduced calcification of marine plankton in response to increased atmospheric CO2, Nature, 407, 311–313,, 2000. a

Riebesell, U., Schulz, K. G., Bellerby, R. G. J., Botros, M., Fritsche, P., Meyerhofer, M., Neill, C., Nondal, G., Oschlies, A., and Wohlers, J.: Enhanced biological carbon consumption in a high CO2 ocean, Nature, 450, 545–548,, 2007. a

ROMS/TOMS Group: Source code for the Regional Ocean Modeling System (ROMS, revision 898), available at:, last access: 23 February 2018. a

Saba, V. S., Griffies, S. M., Anderson, W. G., Winton, M., Alexander, M. A., Delworth, T. L., Hare, J. A., Harrison, M. J., Rosati, A., Vecchi, G. A., and Zhang, R.: Enhanced warming of the Northwest Atlantic Ocean under climate change, J. Geophys. Res.-Oceans, 121, 118–132,, 2015. a

Sarmiento, J. L., Hughes, T. M. C., Stouffer, R. J., and Manabe, S.: Simulated response of the ocean carbon cycle to anthropogenic climate warming, Nature, 393, 245–249,, 1998. a

Shadwick, E. H., De Meo, O. A., and Friedman, J. R.: Discrete CO2 System Measurements in the Chesapeake Bay Mainstem between 2016 and 2018, Dataset archived at W&M ScholarWorks,, 2019a. a

Shadwick, E. H., Friedrichs, M. A. M., Najjar, R. G., De Meo, O. A., Friedman, J. R., Da, F., and Reay, W. G.: High-frequency CO2-system variability over the winter-to-spring transition in a coastal plain estuary, J. Geophys. Res.-Oceans, 124, 7626–7642,, 2019b. a

Shchepetkin, A. F. and McWilliams, J. C.: The Regional Oceanic Modeling System (ROMS): A split-explicit, free-surface, topography-following-coordinate oceanic model, Ocean Model., 9, 347–404,, 2005. a

Shen, C., Testa, J. M., Li, M., Cai, W. J., Waldbusser, G. G., Ni, W., Kemp, W. M., Cornwell, J., Chen, B., Brodeur, J., and Su, J.: Controls on carbonate system dynamics in a coastal plain estuary: A modeling study, J. Geophys. Res.-Biogeo., 124, 61–78,, 2019a. a

Shen, C., Testa, J. M., Ni, W., Cai, W. J., Li, M., and Kemp, W. M.: Ecosystem metabolism and carbon balance in Chesapeake Bay: A 30-year analysis using a coupled hydrodynamic-biogeochemical model, J. Geophys. Res.-Oceans, 124, 6141–6153,, 2019b. a, b

Shenk, G. W. and Linker, L. C.: Development and application of the 2010 Chesapeake Bay watershed total maximum daily load model, J. Am. Water Resour. Assoc., 49, 1042–1056,, 2013. a

Smolarkiewicz, P. K. and Margolin, L. G.: MPDATA: A finite difference solver for geophysical flows, J. Comput. Phys., 140, 459–480, 1998. a

Son, S., Wang, M., and Harding, L. W.: Satellite-measured net primary production in the Chesapeake Bay, Remote Sens. Environ., 144, 109–119,, 2014. a, b, c

St-Laurent, P. and Friedrichs, M. A. M.: Associated dataset: Relative impacts of global changes and regional watershed changes on the inorganic carbon balance of the Chesapeake Bay, W&M ScholarWorks,, 2020. a

Tian, H., Yang, Q., Najjar, R. G., Ren, W., Friedrichs, M. A. M., Hopkinson, C. S., and Pan, S.: Anthropogenic and climatic influences on carbon fluxes from eastern North America to the Atlantic Ocean: A process-based modeling study, J. Geophys. Res.-Biogeo., 120, 757–772,, 2015. a

USEPA: Chesapeake Bay total maximum daily load for nitrogen, phosphorus, and sediment, Tech. rep., U.S. Environmental Protection Agency, Chesapeake Bay Program Office, 410 Severn Avenue, Annapolis, Maryland 21403, available at: (last access: 19 July 2020), 2010. a

USEPA: Guide to using Chesapeake Bay Program water quality monitoring data, Tech. rep., U.S. Environmental Protection Agency, Chesapeake Bay Program Office, publication 903-R-12-001 CBP/TRS 304-12, 410 Severn Avenue, Annapolis, Maryland 21403, available at: (last access: 30 March 2018), 2012. a

Ver, L. M. B., Mackenzie, F. T., and Lerman, A.: Carbon cycle in the coastal zone: Effects of global perturbations and change in the past three centuries, Chem. Geol., 159, 283–304,, 1999. a, b

Wakelin, S. L., Holt, J. T., Blackford, J. C., Allen, J. I., Butenshon, M., and Artioli, Y.: Modeling the carbon fluxes of the northwest European continental shelf: Validation and budgets, J. Geophys. Res.-Oceans, 117, C05020,, 2012. a

Wallace, E. J., Looney, L. B., and Donglai, G.: Multi-Decadal Trends and Variability in Temperature and Salinity in the Mid-Atlantic Bight, Georges Bank, and Gulf of Maine, J. Mar. Res., 76, 163–215,, 2018. a

Walsh, J. J., Rowe, G. T., Iverson, R. L., and McRoy, C. P.: Biological export of shelf carbon is a sink of the global CO2 cycle, Nature, 291, 196–201,, 1981. a, b

Wang, Z. A., Wanninkhof, R., Cai, W. J., Byrne, R. H., Hu, X., Peng, T. H., and Huang, W. J.: The marine inorganic carbon system along the Gulf of Mexico and Atlantic coasts of the United States: Insights from a transregional coastal carbon study, Limnol. Oceanogr., 58, 325–342,, 2013. a

Wanninkhof, R.: Relationship between wind speed and gas exchange over the ocean, J. Geophys. Res., 97, 7373–7382,, 1992. a

Weiss, R. F.: Carbon dioxide in water and seawater: The solubility of a non-ideal gas, Mar. Chem., 2, 203–215,, 1974. a

Xu, J. and Hood, R. R.: Modeling biogeochemical cycles in Chesapeake Bay with a coupled physical-biogeochemical model, Estuar. Coast. Shelf Sci., 69, 19–46,, 2006. a

Yang, Q., Tian, H., Friedrichs, M. A. M., Hopkinson, C. S., Lu, C., and Najjar, R. G.: Increased nitrogen export from eastern North America to the Atlantic Ocean due to climatic and anthropogenic changes during 1901–2008, J. Geophys. Res.-Biogeo., 120, 1046–1068,, 2015a. a

Yang, Q., Tian, H., Friedrichs, M. A. M., Liu, M., Li, X., and Yang, J.: Hydrological responses to climate and land-use changes along the North American east coast: A 110-year historical reconstruction, J. Am. Water Resour. Assoc., 51, 47–67, 2015b. a

Zeebe, R. E. and Wolf-Gladrow, D.: CO2 in seawater: Equilibrium, kinetics, isotopes, vol. 65, Elsevier Ocenographic Series, 2001.  a

Zhang, X., Roman, M., Kimmel, D., McGilliard, C., and Boicourt, W.: Spatial variability in plankton biomass and hydrographic variables along an axial transect in Chesapeake Bay, J. Geophys. Res.-Oceans, 111, C05S11,, 2006. a

Short summary
Over the past century, estuaries have experienced global (atmospheric CO2 concentrations and temperature) and regional changes (river inputs, land use), but their relative impact remains poorly known. In the Chesapeake Bay, we find that global and regional changes have worked together to enhance how much atmospheric CO2 is taken up by the estuary. The increased uptake is roughly equally due to the global and regional changes, providing crucial perspective for managers of the bay's watershed.
Final-revised paper