The impact of oyster aquaculture on the estuarine carbonate system

Ocean and coastal acidification (OCA) affects shellfish and other calcifying organisms due to the dissolution pressure it exerts on their exoskeleton and the elevated metabolic cost of calcification (Melzner et al., 2011; Waldbusser and Salisbury, 2014; Waldbusser et al., 2015; Gledhill et al., 2015). While there is a significant body of literature demonstrating this sensitivity and discussion around potential impacts on wild and cultured shellfish from OCA, most of the focus has been on OCA resulting from various processes, particularly global carbon emissions (Feely et al., 2004; Doney et al., 2009), local freshwater runoff (Salisbury et al., 2008; Duarte et al., 2013), ecosystem metabolism (Feely et al., 2010; Duarte et al., 2013; Wallace et al., 2014; Pacella et al., 2018; Rheuban et al., 2019), temperature changes (Salisbury and Jönsson, 2018), and upwelling (Feely et al., 2008; Hauri et al., 2009; Hauri et al., 2013; Reum et al., 2016). However, oysters are powerful ecosystem engineers that can transform their environment (Dame, 2011), including in aquaculture settings where populations are large and dense relative to the water body volume (Dumbauld et al., 2009; Coen et al., 2011). The beneficial effects of oyster aquaculture on water clarity and suspended particulate removal (Cranford et al., 2011), nutrient cycling (Kellogg et al., 2013; Kellogg et al., 2014; Lunstrum et al., 2018; Ray et al., 2019), and carbon sequestration/removal during harvest (Ahmed et al., 2017; Fodrie et al., 2017) are well known, and oysters are considered one of the most sustainable ways to provide animal protein (Ray et al., 2019; Naylor et al., 2021). Given these environmental benefits and a high demand for oysters, global and U.S. oyster aquaculture has risen rapidly over the last several decades (Naylor et al., 2021; NOAA Fisheries, 2021).

Extensive research shows that calcification and respiration in natural populations release CO₂ into the ocean and atmosphere (Frankignoulle et al., 1995; Gattuso et al., 1995; Gattuso et al., 1999; Zeebe and Wolf-Gladrow, 2001). However, only a few studies have focused on the potential for commercial-scale shellfish aquaculture to affect the surrounding seawater carbonate system (Mistri and Munari, 2012; Munari et al., 2013; Filgueira et al., 2015; Morris and Humphreys, 2019; Li et al., 2021), particularly in the context of shellfish aquaculture lowering surrounding carbonate saturations states (Ω) and limiting its own further expansion through decreased growth rates. Additionally, a sought-after characteristic of oysters on the half shell, a large proportion of the oyster market, is a strong shell, one that is easily shucked and does not break or chip (Mizuta and Wikfors, 2019). A reduction in the shell thickness or hardness from decreased Ω would be detrimental to the quality of the product. Therefore, considering the effect of these ecosystem engineers on estuarine carbonate chemistry to ensure environmentally sustainable and profitable expansion of oyster aquaculture is important.

Dissolved inorganic carbon (DIC) and total alkalinity (TA) are two of the most commonly measured parameters of the carbonate system in marine waters. DIC is the sum of concentrations of all inorganic carbon species and includes CO₂, carbonic acid (H₂CO₃), bicarbonate ion (⁠$HCO3–$⁠), and carbonate ion (⁠$CO32–$⁠). Changes in the concentration of these species can affect the DIC of the system. Photosynthesis, respiration, calcification, and dissolution are large drivers of DIC variability, particularly in coastal areas. TA is a measure of the buffering capacity of a water body defined by Dickson (1981) as:

where the ellipses signify weak organic acids and bases that may contribute to TA, particularly in coastal areas. TA is thus significantly affected by calcification and dissolution, as well as processes that impact nutrient concentrations.

Oysters use calcium and carbonate to build their shells via either aragonite as larvae or calcite as juveniles and adults; both are mineral forms of calcium carbonate (CaCO₃). The Ω of either calcite or aragonite indicates how easily that mineral form can be precipitated through physical processes or by organisms.

where [Ca²⁺] and [$CO32–$] represent ion concentrations of calcium and carbonate, and $Ksp′$ is the apparent solubility product of either aragonite or calcite. Oysters directly reduce the Ω of the surrounding water via two pathways: calcification and respiration (Zeebe and Wolf-Gladrow, 2001). Shellfish build and maintain organic tissue by ingesting food (e.g., phytoplankton and detritus) and then convert some of the particulate organic carbon in the food into soft tissue, excreting a fraction of it as organic carbon waste and respiring the rest to CO₂. The respired CO₂ can react with surrounding water to form H₂CO₃ that lowers pH by dissociating into $HCO3–$ and a free proton. Protons bond with $CO32–$⁠, reducing their concentration and lowering the saturation state of calcium carbonate (Doney et al., 2009). This process increases the DIC/TA ratio because the CO₂ addition increases DIC but not TA. The buffering capacity of the system reaches a minimum when the DIC/TA ratio ≈ 1 (Wang et al., 2013).

Calcification changes TA with a decrease of 2 moles of TA per mole of shell CaCO₃ formed (Lejart et al., 2012). Skeletal calcium carbonate formation removes $CO32–$ and therefore TA and DIC. Calcification also converts existing DIC toward a more CO₂ rich solution, lowering pH, via a reaction such as:

(Frankignoulle et al., 1995; Gattuso et al., 1995; Gattuso et al., 1999). Thus, both respiration and calcification produce CO₂ that lowers Ω in surrounding waters. Unlike the open ocean, many estuaries are already supersaturated with CO₂ due to elevated organic matter remineralization (Wallace et al., 2014; Cai et al., 2017), and any further addition of CO₂ increases the flux from the surface water to the overlying air, which can help reduce local acidification. The resultant Ω will depend on the ability of the system to release CO₂ produced by oyster growth into the atmosphere in addition to the impact of other processes such as photosynthesis, dissolution, mixing, and nutrient cycling (Mistri and Munari, 2012; Morris and Humphreys, 2019). Unless all of the CO₂ produced by oysters is removed, the Ω will decrease relative to the initial condition (Ware et al., 1992), creating a less favorable environment for future calcification.

Scaling these organismal effects to industry level within a growing area is complex but important given the increasing role of aquaculture to supply shellfish to market, especially where cultured bivalves are highly concentrated in productive areas. Prior studies have focused on how biotic and abiotic factors affect metabolic rates and used these relationships to scale to ecosystem level energy transfer (Enquist et al., 2003); here, we simplify and focus on the aquaculture-estuarine exchange of inorganic carbon at the industry level. Specifically, the ability of oysters to change the Ω suggests that their own growth should be considered in assessing the inorganic carbon dynamics of a water body that supports oyster aquaculture.

In estuaries that fringe the Gulf of Maine, the oyster aquaculture industry has expanded rapidly. In the last 15 years, oyster aquaculture in Maine alone has grown from 2 million to nearly 14 million oysters harvested per year (Figure 1). A recent analysis predicted that the robust market for eastern oysters (Crassostrea virginica) in the Northeast United States would lead to an expansion in the annual demand from 7 million pieces (roughly half the present harvest level in Maine) to 26 million pieces in 2030 (The Hale Group, Ltd and Gulf of Maine Research Institute, 2016). The footprint of this industry has also expanded. In 2007, there were 44 active limited purpose aquaculture licenses in Maine, and in 2020 there were 769 (Maine Department of Marine Resources, 2021). Clearly the industry is growing, though not homogeneously; most growth has occurred in a few estuaries. From 2011 to 2017, the harvest from the Damariscotta River Estuary in Maine increased by 1.2 million pieces per year (Figure 1). Since 2003, when Maine Department of Marine Resources began compiling harvest data from lease and license holders, the upper Damariscotta River Estuary has consistently been the major hub for oyster aquaculture in northern New England (67% of total oyster production in Maine in 2019; Maine Department of Marine Resources, 2020). Therefore, understanding the dynamics between the industry and local estuarine conditions to enable sustainable expansion to meet the expected increase in demand is vital to the industry, local economies, and stakeholders in these communities.

In this study, we explored the potential magnitude of impact that net oyster calcification and subsequent harvest can have on estuarine carbonate chemistry to quantify whether oyster shell growth should be included as a term in carrying capacity models for shellfish. Our goal was to understand how current oyster growth might limit further industry expansion, and how future climate conditions may moderate or exacerbate the industry effect. We focused primarily on Ω_CA because juvenile and adult oysters predominantly use calcite to build their shells. Although we model oyster respiration, one could assume that respiration of suspended organic matter would occur in the upper Damariscotta River Estuary at similar rates regardless of whether the shellfish aquaculture industry was present or not. Therefore, our exploratory model was primarily concerned with calculating the carbon impact of calcareous versus noncalcareous heterotrophic communities in the hydrographic conditions of the shellfish growing area.

We used the Damariscotta River Estuary oyster growing area as an example estuary because of its rapid growth in oyster aquaculture. We modeled varying industry sizes, estuarine residence times, and future OCA conditions to determine the current potential effects and how these may change in the future with industry expansion and climate change. We assumed that oyster growth was fueled by external sources of calcium, DIC, and particulate organic matter. Oysters transformed these nutrients into CaCO₃, oyster particulate organic matter (i.e., tissue), and CO₂. The CaCO₃ and oyster tissue were permanently removed from the model domain during harvest. The waste CO₂, a byproduct of CaCO₃ precipitation and metabolism, either remained in the water or was removed via oceanic exchange or air-sea exchange. CO₂ that remained in the water reduced local Ω, potentially hindering future shellfish growth.

The growing area (Figure 2) is highly productive because of the relatively long residence time and shallow depth (average 5.5 m) that incubates Gulf of Maine water, increasing its temperatures during the growing season (6.5°C above water at the mouth of the estuary during peak temperatures) and retaining nutrients and particles such as phytoplankton (Mayer et al., 1996; Lieberthal et al., 2019; Newell et al., 2021). The Gulf of Maine represents the principal source of water for the ocean-dominated Damariscotta River Estuary because the largest freshwater source only contributes 1–3 m³ s^–1 on average during the growing season (McAlice, 1977; and unpublished data from Damariscotta Mills Dam, Kruger Energy, 2015–2018; Figure 2). A Land-Ocean Biogeochemical Observatory buoy (LOBO; Jannasch et al., 2008), positioned within the growing area, collected hourly temperature, practical salinity, and dissolved oxygen data throughout the season. Additionally, we measured physiochemical properties (TA, DIC, practical salinity, and temperature at 1-m depth) of the growing area and Gulf of Maine source water by sampling at the LOBO buoy and mouth of the Damariscotta River Estuary where it meets the Gulf of Maine, respectively (Figure 2), monthly to biweekly from June to September 2018 (6 samplings total). A profile of the whole water column in both locations measured temperature, practical salinity, dissolved oxygen, pH, and chlorophyll a during each sampling; it confirmed that the estuary was well-mixed tidally throughout the year.

Salinity, temperature, pH, DO, and chlorophyll a were measured using a YSI 6920 multiparameter sonde that was calibrated the day prior to sampling using YSI solutions or recommended calibration practices. Triplicate TA and DIC samples were collected using a Niskin sampler from 1-m depth and transferred to 60-ml borosilicate glass bottles using clean flexible tubing. The tube was placed at the bottom of the bottle and the bottle was inverted to allow water to coat all sides to avoid forming bubbles. The bottle was then turned upright, and water was allowed to overflow the bottle for 10 s, or approximately 250 ml. Less than 1% headspace was left in the bottles, which were kept on ice until return to the laboratory refrigerator before analysis. Samples were analyzed within 24 h of collection.

TA was analyzed using potentiometric titration in an open cell titration, following the methods described in Guide to Best Practices for Ocean CO₂ (Dickson et al., 2007). Briefly, samples were kept in a water bath at 25°C for 10 min prior to analysis, and a weighed sample was then placed in a jacketed beaker maintained at 25°C. A calibrated liquid junction pH electrode was submerged in the sample with a magnetic stir bar that gently mixed the sample continuously; 0.1N hydrochloric acid was added to the sample at a rate slower than 0.05 ml s^–1 until the pH reached 2.9. The titration data were then used with the package SeaCarb in R (Gattuso et al., 2020) to calculate the TA of the system. The average standard deviation across all triplicate samples for TA was ± 9.7 μmol kg^–1. DIC was analyzed using a Shimadzu TOC-4200 analyzer. A known volume of room temperature seawater was drawn into a chamber where it was acidified with 8.5% phosphoric acid solution and the off-gassed CO₂ was extracted using a pure nitrogen gas stream. CO₂ was then measured using a nondispersive infrared radiation detector. The full methods can be found in O’Sullivan and Millero (1998). The average standard deviation across all triplicate samples for DIC was ± 6.0 μmol kg^–1. Certified reference materials from the A Dickson Lab at Scripps Institution of Oceanography were used to standardize both TA and DIC samples.

All oyster aquaculture in Maine occurs in subtidal waters. In the growing area, some companies start the oysters in surface rafts and then spread animals onto the bottom for the second year of growth, while the majority of companies use surface culture for both years until harvest (C Davis of Pemaquid Oyster Co, personal communication, June 6, 2019). The oyster growing season is generally when the water temperature is above 10°C, most often spanning early May to late October. The growth of oysters over both first and second growing seasons was estimated as the average change in shell height over the season (Figure 3) with input from industry partners (C Davis of Pemaquid Oyster Co and M White of Mook Sea Farm, personal communication, October 13, 2020). These two companies have been farming in the growing area for over 30 years and provided an interannual average change in oyster shell height for both first- and second-year oysters that they grow. The two estimates were averaged to produce the change in shell height by year used in the model. Calcium carbonate uptake (i.e., shell deposition—dissolution) was estimated by applying shell height to shell dry weight relationships developed in the growing area (Figure 3). We measured seed oyster shell height (2–20 mm) and dry shell weight in 2017, and Adams et al. (2019) measured oyster shell height (40–100 mm) and dry shell weight in 2017. Here, shell height is defined as the distance between the bottom-most section of the hinge to the middle of the top edge (Figure 3).

We estimated the dry shell weight of oysters between 2 and 20 mm using the equation derived from our growth study:

We estimated dry shell weight of oysters between 21 and 90 mm using the shell heights and dry weights measured by Adams et al. (2019):

We assumed dry shell weight was 100% CaCO₃ and that other components of shell (e.g., protein matrix, periostracum) were negligible (Marie et al., 2011). Because these oysters are harvested and then largely exported to other watersheds, this net growth represents net calcification within the growing area.

Harvest throughout the growing season removes second-year oysters as they reach different size classes, reducing the total number of oysters left in the system to calcify and respire. To incorporate the impact of harvest on carbonate removal, we assumed that 25% of oysters are harvested at “cocktail” size (65–75 mm), 50% of oysters are removed at “select” size (75–85 mm), and the remaining 25% are harvested as “jumbo” oysters (>85 mm; C Davis of Pemaquid Oyster Co, personal communication, October 8, 2020).

We developed a one box model (Figure 3) of the impact of oyster net calcification and gross respiration on carbonate chemistry change using Damariscotta River Estuary input values of oyster shell height, temperature, practical salinity, TA, DIC, atmospheric pCO₂, and wind (Table 1). We used the seasonally averaged values of TA, DIC, and practical salinity measured at the mouth of the estuary (Figure 2) as the initial and ocean boundary water chemistry conditions for the model. The freshwater contribution to the estuary was not included due to its extremely low input during the growing season. We therefore kept salinity constant throughout the model runs to avoid errors that would arise from using TA and DIC measured at the mouth of the river with slightly higher salinity than in the growing area. Daily averaged temperature data from the LOBO buoy were used to capture seasonal changes in radiative heating and cooling.

Model variables included TA, DIC, industry size, water residence time in the growing area, in situ pCO₂, Ω_CA, and Ω_AR (Table 2). In situ pCO₂ was used to calculate air-sea exchange and adjust remaining DIC concentrations. The model used a daily time step and estimated daily growth rates for oysters in the upper Damariscotta River Estuary to calculate changes in TA and DIC while accounting for air-sea and oceanic exchanges. We focused the model on the potential impact from oyster growth on carbonate chemistry. We did not include other processes affecting the carbonate system, such as primary production, benthic processes, dissolution, and nutrient dynamics, but their potential contributions are discussed.

For each mole of $CO32–$ removed for calcification, the model removed 2 moles of TA and 1 mole of DIC. Lejart et al. (2012) found that molar CO₂ production from respiration was 3 times higher than the net molar calcium carbonate uptake for adult Pacific oysters, Crassostrea gigas. Assuming this ratio for Crassostrea virginica, water column DIC was increased by 3 moles from respiration for each mole of CaCO₃ derived from shell growth. To scale these estimates to commercial aquaculture scenarios, we modeled between 0 and 50 million oysters in the system to represent a spectrum of no aquaculture to a potential future industry size in the growing area (Table 2). When oysters were present, they were modeled as continuously submerged (i.e., no desiccation regime) for their entire life span in the system. The 2018 industry size was estimated as 16 million oysters by assuming that the 2018 oyster landings from the Damariscotta River (approximately 8 million oysters; Maine Department of Marine Resources, 2020) represent half of cultured oysters in the growing area, given 2 years of growing time required for seed to reach market size (C Davis of Pemaquid Oyster Co, personal communication, January 24, 2018). A potential industry size of 50 million oysters was chosen for the year 2050 by extrapolating the current growth rate (approximately 1 million more oysters in the growing area each year) to the year 2050. Although other factors such as physical space, food supply, and competing uses may limit this growth (Johnson et al., 2019), we modeled this possible industry size to aid our exploration of the history of aquaculture in this system, as well as market and policy implications for shellfish aquaculture development in the Gulf of Maine (The Hale Group, Ltd and Gulf of Maine Research Institute, 2016).

We used the model to isolate the impact of individual processes by modeling the changes in Ω resulting from alterations of one process at a time. For example, oyster net calcification effects on TA, DIC, and Ω were modeled by excluding respiration and oceanic exchange and holding temperature and salinity constant.

To explore the role that OCA may have in changing Ω and the interaction with oyster growth, we examined the combined effects of net calcification and gross respiration on Ω_CA with changes in temperature and DIC due anthropogenic increases in CO₂ projected for the year 2050, using the “business as usual,” RCP 8.5 emissions scenario (IPCC, 2014). Gulf of Maine DIC is predicted to increase an average of 1.25 µmol kg^–1 yr^–1 from increasing atmospheric CO₂ (Siedlecki et al., 2021). Siedlecki et al. (2021) utilized a coupled physical-biogeochemical model developed by Gehlen et al. (2015) and a physical model (Regional Ocean Modeling System; Shchepetkin and McWilliams, 2005) paired with an empirical model developed by McGarry et al. (2021) to estimate changes in DIC, TA, and Ω in the Gulf of Maine in 2050 using RCP 8.5. We therefore adjusted the initial DIC concentration from the current concentration of 1926 µmol kg^–1 to the predicted value of 1966 µmol kg^–1. We also increased the atmospheric pCO₂ to 550 ppm in accordance with estimates from RCP 8.5 (IPCC, 2014). The temperature of the Gulf of Maine has been increasing rapidly this century (Friedland et al., 2020a; Friedland et al., 2020b; Gonçalves Neto et al., 2021) and is expected to continue warming. Brickman et al. (2021) also used the same downscaled global models as Siedlecki et al. (2021) to estimate that the Gulf of Maine will warm by an additional 0.5°C by 2050 compared to average 2018 temperatures using the RCP 8.5 emissions scenario. Thus, we added 0.5°C to the daily 2018 temperatures for the 2050 model scenarios.

We estimated the 2018 residence time for the growing area using a realistic 3-dimensional hydrodynamic model we developed for the mid-coast of Maine, including the Damariscotta River Estuary (see Text S1 for full description). We modeled a range of residence times (i.e., oceanic exchange; Figure 3, Table 2) from 6 to 20 days to explore model sensitivity and reflect that how residence time is calculated (e.g., freshwater replacement time, tidal prism method, isohaline decomposition analysis) can affect our understanding of how Ω varies. Additionally, climate change is predicted to reduce estuarine residence times through increased precipitation and therefore flushing (Scavia et al., 2002) and increased sea level rise and tidal changes (Du et al., 2018). The model was run for each combination of residence time and industry size (Table 2). The CO2SYS program (van Heuven et al., 2011) was used to estimate the daily Ω using the variable temperature, TA, DIC, and average seasonal salinity. We used the K₁ and K₂ dissociation constants derived by Mehrbach et al. (1973) and refit by Dickson and Millero (1987), the KSO_4– disassociation constants by Dickson et al. (1990), and total borate constant of Uppström (1974). Thus, the model uses 3 key inputs to estimate the impact of oyster growth and physical processes on Ω (Figure 4): (1) oyster shell daily weight gain (g CaCO₃ day^–1) of both first- and second-year oysters combined, (2) daily average water temperature (°C), and (3) wind speed (m s^–1) measured at the nearby Wiscasset Airport (Menne et al., 2019).

Finally, we modeled a scenario where we excluded the impact of oyster respiration on Ω. This scenario assumes the CO₂ produced through respiration was utilized for photosynthesis within the growing area or removed through another avenue not modeled here rather than directly released into the water column. All other model forcing functions (i.e., air-sea and oceanic exchange formulations, and seasonal temperature) and scenarios (i.e., industry sizes, residence times, and climate conditions) were as described above.

Several assumptions were necessary to simplify model setup and analysis. First, we omitted any freshwater entering the growing area because the largest freshwater source only contributes 1–3 m³ s^–1 on average to the growing area. Variable salinity should be included for other systems where freshwater is a larger contributor. Second, we omitted changes in TA due to ammonia excretion by oysters. Data from this estuary has shown that there has not been an increase in the standing stock of ammonia over the last 25 years (Mayer et al., 1996, and unpublished data from K Thornton, University of Maine, 2017–2019) despite a more than 10-fold increase in the number of oysters in the system (Maine Department of Marine Resources, 2020). Although excretion of ammonia by oysters may lead to an increase in TA (Lejart et al., 2012) on short time scales (hours to days), the ammonia is removed from the system quickly either through primary production or oxidation to yield an assumed net-zero effect on TA. Third, dissolution of dead shells can also be a source of TA to an estuary from natural oyster reefs (Waldbusser et al., 2013; Shen et al., 2019). However, we did not incorporate dissolution because we based our industry sizes on harvest data; thus, all modeled oysters were removed from the system through harvest and exported and did not remain in the growing area where their dissolution could affect Ω. No shells are returned to the estuary after harvest per state law. Loss of oysters (i.e., settling to the bottom to remain in the estuary) during grow-out and harvest is thought to be very low but might provide neutralization of acidification over the decades; however, we have no data to include this process quantitatively in a model of this type. Dissolution should be included in models of systems with natural oyster reefs, when dead shells are returned to the water body, or calculations in aquaculture systems where the amount of product loss is known or well estimated.

Finally, oysters can use a fraction of the CO₂ produced through respiration as the carbon source for calcification rather than simply releasing the CO₂ to the water column. Previous studies have shown that, on average, 10% of bivalve shell mass is derived from respired CO₂ (McConnaughey et al., 1997; Lorrain et al., 2004; Gillikin et al., 2006), although Gillikin et al. (2007) found a significantly higher fraction (5%–37%) in hard clams (Mercenaria mercenaria) from North Carolina. Given the uncertainties in the amount of respired carbon used for calcification in eastern oysters, and the estimated small fraction used by other bivalves, we did not incorporate the impact from metabolic carbon usage on DIC, but we see value in improving this model in the future by incorporating a more detailed carbon tracking process for oysters as requisite details become available.

The $CO32–$ removed per oyster was calculated using the average shell weight gain each day, converted from grams to micromoles. We used the following equation to estimate total $CO32–$ uptake:

where CO_3uptake is the amount of $CO32–$ removed from the water for shell growth (µmol kg^–1 day^–1), shell_wg is the amount of shell gained each day by both year-1 and year-2 oysters averaged together (µmol CaCO₃ day^–1 oyster^–1), estimated using Equations 4 and 5. Density_GA is the average water density in the growing area calculated from temperature and salinity (kg-water l-water^–1; Moataz, 2021), volume_GA is the average volume of the growing area (l-water; McAlice, 1977), and industry size is the number of oysters in the growing area each day for both first- and second-year oysters.

We calculated TA each day using:

where TA_(t) is the total alkalinity after net calcification, F is the fraction of water replaced each day in the growing area due to tidal exchange (F = t/Residence Time), TA_(t–1) is the TA from the previous day, CO_{3 uptake calc} is the amount of $CO32–$ removed from the growing area each day by oyster calcification, and TA_OE is the average TA at the mouth of the Damariscotta River Estuary, held constant at 2091 µmol kg^–1. All concentrations are in µmol C kg^–1.

DIC was calculated similarly to TA but with the addition of oyster respiration and air-sea exchange. DIC was reduced by 1 mole for each mole of calcification and was increased by 3 moles from respiration for each mole of net calcification. Therefore,

where DIC_(t) is the concentration of DIC at the end of the day after all processes are included, F is the fraction of water replaced each day due to oceanic exchange, DIC_(t–1) is the DIC concentration from the previous day, CO_{3 uptake calc} is the daily loss of DIC from net calcification, CO_{3 uptake resp} is the daily addition of DIC from oyster respiration, DIC_ASE is the daily amount of DIC exchanged with the atmosphere due to gas transfer, and DIC_OE is the average DIC at the mouth of the Damariscotta River Estuary (1926 µmol kg^–1). All variables are in µmol C kg^–1.

We used TA and DIC calculated from daily oyster shell growth and respiration to calculate the partial pressure of carbon dioxide (pCO₂) in the water using the program CO2SYS (van Heuven et al., 2011) and the same constants as above. We used the “tracers only” gas transfer velocity (k₆₀₀) equation (Raymond and Cole, 2001) and converted to a saltwater formulation (k₆₆₀) using the calculated Schmidt number and coefficients from Wanninkhof (2014). We chose the Raymond and Cole (2001) gas transfer velocity equation because it was formulated for estuaries with wind speeds between 0 and 7 m s^–1 which are typical in the growing area. The gas transfer velocity was calculated as:

where x is the wind speed (m s^–1) measured at 10 m above the ground and Sc is the Schmidt number calculated each day. The air-sea flux was calculated using:

where CO_{2 sol} is the solubility of CO₂ calculated from temperature, salinity, and pressure, pCO_{2 w} is the calculated in situ pCO₂, and pCO_{2 atm} is the atmospheric CO₂ concentration.

Using the calculated partial pressure of carbon dioxide in the water, water temperature, practical salinity, wind, and an average atmospheric pCO₂ concentration (409 and 550 ppm for 2018 and 2050, respectively), we estimated the amount of CO₂ that transfers to or from the atmosphere each day and adjusted the DIC pool to reflect this flux. The 2018 atmospheric pCO₂ concentration was measured on a buoy in the western Gulf of Maine (NOAA Pacific Marine Environmental Lab/University of New Hampshire Coastal CO₂ Buoy D).

First-year oysters were estimated to grow from 2 to 45 mm and second-year oysters from 46 to 90 mm between early May and late October. Daily oyster shell weight gain assuming equal numbers of first- and second-year oysters combined was 0.0101 g CaCO₃ oyster^–1 day^–1 in the beginning of May, and increased to a maximum of 0.25 g CaCO₃ oyster^–1 day^–1 in early September (Figure 4a). Temperature varied on a seasonal cycle, beginning at 10.9°C in early May, peaking at 24°C in early August, and ending at 10.6°C in late October (Figure 4b). Daily average wind speed was low (0–2 m s^–1) but representative of the relatively protected upper estuary (Figure 4c). The measured DIC and TA concentrations in the center of the growing area for the beginning, middle, and end of the season were 1963, 1960, and 1798 µmol kg^–1, respectively, for DIC and 2166, 2084, and 2129 µmol kg^–1, respectively, for TA. These measured values were corrected to a salinity of 31.36 to match the model salinity, but they incorporate impacts from other estuarine processes not modeled here. The modeled DIC and TA for the beginning, middle, and end of the season in the growing area were 1933, 1885, and 1929 µmol kg^–1, respectively, for DIC and 2091, 2080, and 2080 µmol kg^–1, respectively, for TA.

The modeled DIC decreased from 1929 µmol kg^–1 to 1888 µmol kg^–1 in the first half of the growing season regardless of oyster presence (Figure 5a), largely due to degassing driven by increased temperature. During the second half of the growing season, the modeled DIC with both 0 and 16 million oysters increased toward the initial concentrations, reaching 1939 and 1933 µmol kg^–1, respectively, by the end of October. Between mid-August and the end of October, 16 million oysters lowered the DIC by 3.4 µmol kg^–1 on average (Figure 5a). The TA remains constant throughout the zero oyster scenario because oyster growth is the only process removing TA in the model. In the 16 million oyster scenario, the TA decreases slowly between early May and mid-June (<5 µmol kg^–1 lost by mid-June), and then decreases more rapidly until early August when it levels off at 2080 µmol kg^–1 for the rest of the season (Figure 5a). The Ω_CA increased by 0.6 units between early May and early August reaching a maximum of 3.63 in the zero oyster scenario, whereas the Ω_AR increased by 0.5 units during the same time period, with a maximum of 2.37. The inclusion of oysters lowered Ω_CA and Ω_AR by 0.11 and 0.07 units, respectively, on average between July and the end of October, with a maximum difference of 0.18 and 0.11 units, respectively, in early August. While the trends in DIC, TA, and Ω remained the same for the 2050 OCA model, the average DIC for the growing season increased and Ω decreased (Figure 5c and d). The maximum change in DIC over the season increased from 45.1 µmol kg^–1 in 2018 to 48.9 µmol kg^–1 in 2050. The TA was not altered by the 2050 conditions. The average Ω_CA and Ω_AR for the growing season decreased by 0.55 and 0.35 units, respectively. However, the maximum change in both Ω_CA and Ω_AR over the growing season did not vary between 2018 and 2050.

The net changes in these carbonate parameters result from accumulated impacts of daily biological and physical processes. Oyster respiration and net calcification had minimal impact (<0.05 µmol kg^–1 day^–1) on DIC during slow oyster growth before July, but these impacts increased throughout the season, with respiration increasing ▵DIC twice as much as calcification decreased it (2 and 1 µmol kg^–1 day^–1, respectively; Figure 6a). Air-sea exchange of CO₂ due to oyster growth reduced ▵DIC by a similar magnitude as oyster respiration increased it, relative to the initial conditions. However, the air-sea exchange of CO₂ lagged behind production via respiration by 1 week (Figure 6a). Oceanic exchange raised DIC by 1.6 µmol kg^–1 day^–1 on average, mirroring losses from air-sea exchange. TA was altered by two counterbalancing processes, net calcification and oceanic exchange. Net calcification was slow before July with minimal impact on TA (<0.33 µmol kg^–1 day^–1). However, between July and late October, TA decreased by 0.90 µmol kg^–1 day^–1 on average. Oceanic exchange increased TA by 0.82 µmol kg^–1 day^–1 on average, which restored TA back toward the boundary conditions at the current industry size (Figure 6b).

The daily temperature impact on Ω_CA was net positive (0.0017 units day^–1 on average) over the first half of the season, and net negative (–0.0017 units day^–1 on average) over the second half of the season (Figure 6c), mirroring the seasonal change in temperature. Non-oyster driven air-sea exchange initially decreased Ω_CA by up to 0.06 units day^–1 as the water was undersaturated in the beginning of May and absorbs CO₂ from the atmosphere. However, after mid-May, air-sea exchange increased Ω_CA by up to 0.04 units day^–1 until late October when air-sea exchange lowered Ω_CA by up to 0.06 units day^–1 again. The air-sea exchange pattern is driven mostly by seasonal temperature changes (Figure 6c). Oceanic exchange counterbalanced change in Ω_CA from oyster growth and air-sea exchange throughout the season; ▵Ω_CA values of up to –0.03 units day^–1 from oceanic exchange occurred in August when DIC was being lost to the atmosphere at a higher rate than TA and DIC were being lost due to calcification. Net oyster calcification and gross respiration have minimal impact prior to July but decreased Ω_CA in the latter half of the season by up to 0.01 and 0.026 units day^–1, respectively (Figure 6d). Air-sea exchange due to oyster growth was negligible for the first part of the growing season but became more substantial (increasing Ω_CA by 0.025 units day^–1) toward the end of the growing season.

Increasing industry size and 2050 climate conditions (assuming a RCP 8.5 scenario) both significantly lowered Ω, particularly during the latter half of the growing season (Figure 7a and b). Modeling Ω using only the physical processes (zero oysters) and the 2050 climate conditions decreased Ω_CA and Ω_AR by 0.57 and 0.37 units, respectively, on average over the growing season compared to the 2018 climate conditions with zero oysters. Increasing the industry size to 50 million oysters caused net calcification and gross respiration to decrease Ω_CA and Ω_AR by 0.16 and 0.10 units on average under 2018 climate conditions. Alternatively, keeping the industry size constant at 16 million oysters, but under modeled 2050 carbonate chemistry conditions, lowered Ω_CA and Ω_AR by an average of 0.55 and 0.35 units, respectively, over the growing season with a maximum reduction of 0.57 and 0.36 units, respectively, in early August. The combination of 2050 carbonate chemistry conditions and 50 million oysters lowered Ω_CA and Ω_AR by 0.67 and 0.43 units, respectively, on average for the latter half of the growing season when compared to 2018 climate conditions and industry size.

Variation in water residence time had a smaller effect on Ω than increasing industry size and 2050 climate conditions. In the first third of the growing season under 2018 conditions, residence time had almost no effect on Ω, but from July to September, decreased residence time resulted in lower Ω. Then, for the month of October, decreased residence time increased Ω (Figure 8a). Residence time had minimal impact on Ω when the industry size increased to 50 million oysters and modeled with 2018 climate conditions, with lower residence time slightly lowering or raising Ω at different times in the season (Figure 8b). With 2018 industry size and 2050 climate conditions, shorter residence time led to decreased Ω throughout the season (Figure 8c). The model combination of 2050 climate conditions and 50 million oysters was similar to the 2050 scenario with the 2018 industry size (Figure 8d).

Oyster growth affects Ω through net calcification and respiration and was sensitive to industry size and OCA conditions (Figure 9). In 2018 scenarios (Figure 9a and c), oyster growth lowered Ω relative to the zero oyster scenario regardless of how many oysters were included in the simulation. Calcification of 16 million oysters lowered Ω_CA and Ω_AR by 0.03 and 0.02 units, respectively, on average, whereas calcification and respiration combined lowered Ω_CA and Ω_AR by 0.08 and 0.05 units, respectively, on average over the growing season. Calcification of 50 million oysters reduced Ω by nearly the same magnitude as calcification and respiration of 16 million oysters. The same trends were observed for scenarios modeled with 2050 climate conditions, but the relative difference between impacts of the various scenarios decreased as overall Ω decreased (Figure 9b and d).

Ocean acidification research has consistently shown that most calcifying organisms will likely be negatively impacted by anthropogenic alterations in carbonate chemistry (Miller et al., 2009; Gazeau et al., 2013; Waldbusser and Salisbury, 2014; Gledhill et al., 2015; Waldbusser et al., 2015; Ekstrom et al., 2015). However, studies rarely consider aquaculture settings where calcifying organisms can enhance acidification of their own surroundings (Mistri and Munari, 2012; Munari et al., 2013; Morris and Humphreys, 2019; Li et al., 2021), with implications for the sustainable carrying capacity. This modeling study shows how these organisms potentially influence estuarine carbonate chemistry in addition to the atmospherically mediated changes. This quantitative framework highlighted that the counterbalancing effects of air-sea exchange and higher temperatures on Ω will make measuring the impact on Ω by oysters in situ difficult. Furthermore, there will be unpredictable interactive impacts of other processes that affect Ω, such as photosynthesis and respiration of other organisms, nutrient fluxes in the water column or from the sediments, and large precipitation events.

Our model estimates that, currently, oysters can decrease both Ω_CA and Ω_AR in the Damariscotta’s growing area by approximately 3% each over the growing season. Although these changes are small, they may be important in systems with high DIC/TA ratios, such as the Damariscotta River Estuary, where ambient Ω_CA and Ω_AR remain close to 3 and 2, respectively, during the growing season. The same amount of oyster growth will drive this system closer to dissolution (Ω < 1) than in a more southern estuary with lower DIC/TA ratios (Morris and Humphreys, 2019). Any further reductions in TA and Ω require that oysters use more energy to create new shell (Melzner et al., 2011). Although juvenile and adult bivalves continue calcifying when Ω_AR is undersaturated, their calcification rate decreases (Gazeau et al., 2007; Miller et al., 2009; Ries et al., 2009; Waldbusser et al., 2010; Melzner et al., 2011). Schwaner et al. (2020) showed that larval and juvenile bivalves raised under high pCO₂, low pH treatments (pCO₂ approximately 1166 ppm, pH = 7.57) had increased susceptibility to disease and increased mortality, and Lord and Whitlatch (2014) showed that the thickness of eastern oyster shells is strongly tied to temperature (i.e., lower temperature increases pCO₂ and reduces Ω). If shell thickness and/or hardness are reduced (Beniash et al., 2010; Meng et al., 2018), bivalves are expected to be more susceptible to predators, such as oyster drills, and may also be less valuable to the half shell market (Mizuta and Wikfors, 2019). First-year “seed” oysters are placed on the farm in the beginning of the growing season when Ω is low due to cold temperatures. The combination of low temperatures and reductions in Ω from oyster growth may lead to lower biomineralization rates by the juvenile oysters leading to increased mortality, particularly at higher industry sizes or future OCA conditions where Ω_AR approaches 1. Thus, understanding how cultured oysters and physical processes alter carbonate chemistry in tandem is important to inform sustainable growth of the industry and maintain a high-quality product.

Physical controls on Ω in this system are as important as oyster growth at the current industry size but can increase and decrease Ω and can be hard to predict. Higher temperatures from climate change will increase Ω due to its thermodynamic properties; however, increased atmospheric and in situ CO₂ concentrations from climate change will decrease Ω. Consistent with other models (Siedlecki et al., 2021), our model predicts that increasing in situ pCO₂ to levels projected for the Gulf of Maine will have a larger impact than increased temperature on Ω, leading to an overall decrease in Ω in the growing area. Our model predicts similar change in Ω (19% with zero oysters present) to models used in Siedlecki et al. (2021), which predicted on average a 22% reduction in Ω during the growing season.

The inclusion of residence time, or oceanic exchange, in our model (Figure 3) was important, as it replenished TA lost to calcification under all scenarios. However, depending on the DIC, the oceanic exchange term either added more DIC or diluted the DIC pool. For scenarios modeled with 16 million oysters, between June and September, the DIC/TA ratio of the incoming ocean water was higher than that in the growing area (Figure 8a and c) because elevated temperatures in the growing area increased CO₂ flux to the atmosphere. However, during May and October in the scenarios with 16 million oysters, or during the second half of the growing season with 50 million oysters and 2018 climate conditions, the incoming water had a lower DIC/TA ratio than the growing area, and reduced the overall DIC, increasing Ω. Thus, the impact of residence time was dependent on the seasonal temperature change and associated air-sea exchange, as well as the increase in oyster growth over the season.

Changes in residence time due to climate change may alter the impacts on Ω in several ways. If freshwater flow were included in the model, increased precipitation to the estuary and watershed would decrease the residence time (Scavia et al., 2002), but also likely lower the TA due to the higher DIC/TA ratio of freshwater, leading to lower Ω. A decrease in residence time due to sea-level rise (Du et al., 2018) would lessen the impact from oysters via increased volume exchange. However, it would also decrease Ω by lowering temperature and thus loss of CO₂ to the atmosphere, especially because the incoming water has a higher DIC/TA ratio than water in the growing area. Given the varied, often competing direct (mixing) and indirect (temperature-induced changes in air-water exchange) effects of altered residence time on Ω, more dynamic model simulations are desirable in the future.

Shellfish are just one part of complicated biogeochemical cycles in estuaries that can affect carbonate chemistry. The ratio of local primary production to calcification in various calcifying systems can control in situ pCO₂ (Suzuki, 1998), and therefore acidification. The impact of photosynthesis and respiration in this system does impact carbonate chemistry on short time scales (days–weeks, unpublished data), but averaged over the growing season, the growing area is slightly heterotrophic (–3.7 mmol O₂ m^–2 day^–1; Figure S2), which may indicate that photosynthesis is not able to capture all of the excess DIC released from oysters. As such, the current impact from oysters in this system is likely somewhere between the calcification-only and calcification-plus-respiration scenarios depicted in Figure 9. However, if the balance between photosynthesis and respiration in the growing area changes due to warming, altered watershed inputs, enhanced oyster grazing, or nutrient limitation, the impact of oyster respiration on Ω may change.

Carrying capacity, from a production perspective, is commonly defined as the maximum number of individuals that a given environment can sustain (Byron et al., 2011; Filgueira et al., 2015). Ecological carrying capacity seeks to balance production with minimal harm to the environment, to preserve pelagic and benthic habitats, and to maintain ecosystem functions (Byron et al., 2011; McKindsey, 2012). Ecological carrying capacity models lag behind production-based shellfish models, and most to date have focused on phytoplankton supply (Byron et al., 2011; Filgueira et al., 2014) or the impacts of biodeposition on benthic habitats (Grant and Filgueira, 2011). The effect of aquaculture on ambient carbonate chemistry has not previously been considered when establishing carrying capacity criteria (McKindsey, 2012; Filgueira et al., 2015; Smaal and van Duren, 2019); however, our model results suggest that oyster aquaculture could impose limits on production or sustainable industry growth due to autoacidification, especially in low alkalinity ecosystems and very high industry sizes. In addition to the domain of our model, future carbonate carrying capacity models may need to incorporate sediment sources of acidity and alkalinity as increased biodeposition alters estuarine respiratory pathways (Thomas et al., 2009; Hu and Cai, 2011). To strike a balance between oyster production and ecosystem integrity, considering the impact of the shellfish aquaculture industry on the carbonate system may be important when developing future carrying capacity criteria, but almost certainly only in areas with Ω currently close to dissolution levels.

Our idealized model simulations suggest that cultured oysters have the ability to alter their surrounding carbonate chemistry through calcification and respiration, where the number of oysters is important in controlling the extent of Ω reduction through TA removal and DIC addition. Physical processes are important in controlling Ω, but the effects of air-sea gas and oceanic exchange depend strongly on the internal dynamics of DIC and TA production and consumption. The combination of future OCA conditions and future industry size will be important determinants of carrying capacity estimates. As shellfish aquaculture continues to expand worldwide, the need to consider not only how local and regional ocean acidification may affect growth of the species but also how shellfish growth itself, when at high densities, may impact the local carbonate budget is clear. An increasing number of biogeochemical models include some variables of the carbonate system, including calcification and dissolution (Shen et al., 2019; Gomez et al., 2020); however, depending on the resolution of the model, they may not capture localized impacts of shellfish aquaculture, and these local conditions are often important when predicting future growth rates and determining an ecological carrying capacity. We suggest that biogeochemical models of areas with emerging or sustained shellfish aquaculture industries include local growth rates and number of shellfish as part of the carbonate cycling variables to increase process-level understanding and, ultimately, predictive capacity.

Data from the LOBO buoy are accessible at: http://maine.loboviz.com/loboviz/.

The supplemental files for this article can be found as follows:

Text S1. Figure S1–2. Docx

The authors thank Cheyenne Adams and Alwyn Ecker for their contribution of oyster shell height data, Kathleen Thornton for her guidance in laboratory methods, and Jon Lewis for his inquiry which inspired this project. They also thank the 3 reviewers who gave encouraging and constructive feedback to this manuscript. This is the University of Maryland Center for Environmental Science Publication No. 6058 and CBL Technical Series CBL 2022-018.

This activity was partially supported by the National Science Foundation award #IIA-1355457 to Maine EPSCoR and National Sea Grant (NA18OAR4170330) at the University of Maine (DC Brady, LM Mayer, and CM Liberti), by a NOAA Regional Vulnerability Assessments for Ocean Acidification award (DC Brady, LM Mayer), and by the United States National Oceanographic and Atmospheric Administration Ocean Acidification Program (NOAA-OAP; award # NA15NOS4780184 and NA18NOS4780179) to the University of Maryland Center for Environmental Science (JM Testa).

The authors have declared that no competing interests exist.

Contributed to conception and design: CML, MWG, DCB, LMM, JMT, WL.

Contributed to acquisition of data: CML, MWG.

Contributed to analysis and interpretation of data: CML, MWG, DCB, LMM, WL.

Drafted and/or revised the article: CML, MWG, DCB, LMM, JMT, WL.

Approved the submitted version for publication: CML, MWG, DCB, LMM, JMT, WL.

How to cite this article: Liberti, CM, Gray, MW, Mayer, LM, Testa, JM, Liu, W, Brady, DC. 2022. The impact of oyster aquaculture on the estuarine carbonate system. Elementa: Science of the Anthropocene 10(1). DOI: https://doi.org/10.1525/elementa.2020.00057

Domain Editor-in-Chief: Jody W. Deming, University of Washington, Seattle, WA, USA

Guest Editor: Andrew Pershing, Climate Central, Princeton, NJ, USA

Knowledge Domain: Ocean Science

Part of an Elementa Special Feature: Gulf of Maine 2050: Visioning Regional Resilience and Sustainability