Projecting ocean acidification impacts for the Gulf of Maine to 2050: New tools and expectations

Ocean acidification (OA) is increasing predictably in the global ocean as rising levels of atmospheric carbon dioxide (CO₂) are absorbed by the ocean air–sea interface, resulting in higher oceanic inorganic carbon concentrations (Gattuso et al., 2015; Le Quéré et al., 2018; Friedlingstein et al., 2019). Over time, ocean uptake of CO₂ perturbs the entire carbonate system, with a concomitant reduction in surface pH, carbonate ion (CO₃^2–), and carbonate saturations states (here, we focus on aragonite saturation state or Ω_a; Bates et al., 2014). Since the beginning of the 19th century, the pH of the world’s surface ocean has decreased by 0.1 units (Caldeira and Wickett, 2003; Orr et al., 2005; Doney et al., 2009), and further pH reductions on the order of 0.2–0.3 units are predicted by 2100 under Intergovernmental Panel on Climate Change (IPCC) AR2 scenarios (Feely et al., 2009) due to the projected increase in atmospheric CO₂.

The carbonate system in coastal waters is impacted not only by changes in atmospheric CO₂ concentrations but also by varying fluxes of dissolved inorganic carbon (DIC), total alkalinity (TA), and nutrients derived from local and remote processes. In the Northwest Atlantic, local processes include low alkalinity river discharge (Salisbury et al., 2008), atmospheric deposition of acidic and alkaline compounds (Doney et al., 2007), sedimentary processes (Fennel et al., 2008), and coastal eutrophication (Wallace et al., 2014; Rheuban et al., 2019). Therefore, the Gulf of Maine (GOM) is a region of confluence for many processes that further affect the carbonate system. Two main regions within the GOM currently experience suboptimal carbonate conditions for many marine organisms (i.e., Ω_a < 1.5)—the nearshore and subsurface deep shelf. In nearshore regions, close to the coastline and less than 50 m in water depth, for example, five large rivers discharge into the GOM and mix with more saline waters in part due to the large tides. In addition, acidic and alkaline compounds from many of the coal-fired power plants upwind of the GOM are deposited in the surface waters (Doney et al., 2007) and watershed, and there is a south-to-north decreasing gradient in both human population density and eutrophication. The interannual to decadal variability of large-scale, more remote circulation patterns, like the Gulf Stream, Labrador Current, and Atlantic meridional overturning circulation, delivers warmer, saltier waters from the south or cooler fresher waters from the north, all of which contribute to temporal and spatial variability in the carbonate chemistry of the GOM. Seasonal variation in the delivery of locally produced and imported organic materials to subsurface deeper shelf offshore regions (beyond 50-m isobaths) contributes to the variability in biologically suboptimal conditions, which may be further influenced by sedimentary processes (Wang et al., 2017). Regional studies over the last decade have identified that these drivers of acidification in the GOM (summarized in the 2015 Maine OA Study Commission; Salisbury et al., 2008; Wang et al., 2013; Strong et al., 2014; Gledhill et al., 2015; State of Maine Legislature, 2015; Salisbury et al., 2018), the relative contributions of each driver, and in particular the mixing of water masses, remain poorly understood.

Sutton et al. (2016) performed an analysis of GOM data using the Pacific Marine Environmental Laboratory (PMEL)/University of New Hampshire (UNH) Coastal CO₂ Buoy D (43.02°N, 70.54°E, “Buoy D”) among other moorings and, after adopting assumptions about past atmospheric CO₂ values, concluded that Ω_a had never dropped below 1.6 in the GOM during preindustrial times. That study also showed that during the period 2010–2013, coastal GOM waters at Buoy D regularly and persistently dipped below the preindustrial minima of Ω_a and pH. This historical view and the observed impacts on mollusk larvae (e.g., Salisbury et al., 2008; Waldbusser and Salisbury, 2014; Ekstrom et al., 2015; Waldbusser et al., 2015) informed our choice of the critical Ω_a threshold of 1.5 in this study. The surface waters of the GOM at Buoy D recently experienced conditions below that threshold 9% of the year using the Buoy D record spanning 2010–2020, with peak exposure to low Ω_a in the winter months (December–March; Figure 1). Nearshore, at the University of New Hampshire Coastal Marine Lab (UNH CML) station (43.07°N, 70.71°E), conditions existed below the Ω_a threshold of 1.5 in the surface waters (Figure 1) 41% of the time over that same time period. The nearshore CML location experiences more freshwater influence than Buoy D, which is about 18 km offshore (Figure 2), and freshwater has a low Ω_a in this area due to the low alkalinity content of local rivers (Hunt et al., 2011). Despite the existence of some oyster growers upstream, the area around the CML location was not an area of prominent shellfish harvest over the observation period presented here (2010–2019).

The history of direct high-quality carbonate chemistry observations in the GOM is relatively short, and each observational program provided for limited spatial and mostly surface water coverage. The deployment of “Buoy D” equipped with an MAPCO2 system (Sutton et al., 2019) was one of the earliest observational efforts in the region. In addition to this system, in 2004, UNH researchers began monthly-to-seasonal transects in the Western GOM, from shore to Wilkinson Basin (WB), making discrete measurements of TA and DIC at various depths, as well as partial pressure of CO₂ (pCO₂) in surface water. Regular, GOM-wide, dedicated research cruises conducting near-synoptic inorganic carbon sampling of the full water column along select transects did not begin until 2007. Since then, these surveys have been repeated every 3–4 years and are supplemented by annual seasonal surveys performed as part of the NOAA National Marine Fisheries Service Ecosystem Monitoring surveys and by individual studies (e.g., Wang et al., 2017). Today, the observational network has expanded and includes many subsurface observations (Figure 2) and community science-based sampling events like “Shell Day” (Rheuban et al., 2021; Figure 2). Such direct observations have indicated that the least buffered waters on the east coast of the United States reside on the Northwest Atlantic shelves, which may make them exceptionally susceptible to acidification (Wang et al., 2013).

In an effort to extend the observational record back several decades, Salisbury and Jönsson (2018) used surface temperature and salinity outputs of a general circulation model applied to empirical algorithms trained from the Surface Ocean CO_2(aq) Atlas (Bakker et al., 2016) and informed by satellite chlorophyll to derive the estimates of past OA conditions in the GOM. The analysis of a 34-year reconstructed record (1981–2014) revealed decades when the carbonate system was perturbed to an extent greater than expected solely from atmospheric uptake of CO₂. In fact, recent temperature and salinity changes in the GOM (between 2005 and 2014) appear to have partially compensated for the effects of OA over the past decade by causing an increase in Ω_arag of 0.4 in surface waters. The analysis also revealed markedly different behavior between pH and Ω_ar--. The pH was found to decline at a rate of –0.018 decade^–1, driven primarily by OA and was largely in agreement with global surface water rates of decline over the past 30 years of –0.0013 to −0.0025 per year (Bates et al., 2014). The OA effect on Ω_a, on the other hand, is partially obscured by increased temperature and buffering by higher salinity waters entering the GOM (Salisbury and Jönsson, 2018). The differing responses of individual carbon variables, such as pH and Ω_ar, to changes in temperature and salinity can lead to varying seasonal amplification or attenuation in these variables globally, which can also impact trends (Kwiatkowski and Orr, 2018; Jiang et al., 2019). Using this reconstructed time series, Salisbury and Jönsson (2018) estimated the “time of emergence” or the time required for the atmospheric CO₂ signal to emerge from background natural variability. They found that observing a discernible OA signal in pH may require 30 years of sustained measurements and up to a century to observe the OA signal in Ω_arag (Salisbury and Jönsson, 2018). Despite these long emergence times, nearshore and offshore GOM regions influenced by local rainfall, freshwater, and respiration are already experiencing event-scale reductions in Ω_a considered suboptimal for many marine organisms.

In the future, OA is projected to cause global annual economic losses up to US$100 billion due to decreased shellfish production, based on the use of the highest CO₂ emission scenarios for 2100 (in this case scenario A1B) by the IPCC to drive spatially coarse resolution global models (Narita et al., 2012). Although oceanographic, social, and economic variability all dictate that the impacts of OA will not be uniform, New England is increasingly dependent on shellfisheries (Figure S1; Gledhill et al., 2015). Large-scale and regional studies show that the combination of global and local drivers of acidification of the waters off New England makes its shellfisheries, both the predominantly offshore wild harvest fisheries and predominantly nearshore aquacultural production, potentially among the most vulnerable to OA in the United States (Ekstrom et al., 2015, Gledhill et al., 2015), with socioeconomic consequences for the communities that rely on these industries. This vulnerability to OA also has implications for the shellfisheries and shellfish-dependent communities of Atlantic Canada. Regional economic and vulnerability analyses, focused on the scallop industry and scallop landings in Massachusetts, found that OA is projected to threaten tens of thousands of jobs and cause hundreds of millions of dollars in losses during the 21st century (Cooley and Doney, 2009; Cooley et al., 2015; Ekstrom et al., 2015; Rheuban et al., 2018).

Global projections with spatial resolution of a degree or more, relied upon by the IPCC, cannot resolve coastal processes with complex feedbacks known to be important to the regional evolution of the OA signal. Some high-resolution global projections have been developed which perform well in some coastal settings (Saba et al., 2016; Alexander et al., 2020) but lack regional biogeochemical processes which can amplify or dampen these global changes, particularly in coastal shelf regions like the GOM.

Here, we use a combination of regional high-resolution simulations that include coastal processes to project OA conditions for the GOM in 2050 and consider what this projection portends for its ecosystems with a particular focus on the nearshore and subsurface environments. First, we review the science specific to acidification of coastal waters (nearshore to shelf break) extending from the Canadian Scotian Shelf southward to the northern reaches of the Mid-Atlantic Bight. Next, we review what is known about the sensitivity of the regional ecosystem inhabitants with a focus on those of socioeconomic value to the region. We then describe the methods employed to compute trends and provide a suite of regional projections. Finally, we discuss the results of the projections and their implications for GOM ecosystems.

The GOM is a productive temperate continental shelf sea bounded by Cape Cod to the west, Georges Bank to the south, and Nova Scotia to the northeast. It is well known for its large semidiurnal tides and their resulting impact on mixing and also for the high commercial value of its fish and shellfish landings. It is separated from the open northwest Atlantic by both Georges and Browns banks. Considerable control on seasonal to interannual circulation patterns is exerted via shelf–sea exchange through the narrow northeast channel (NEC), which separates Georges and Browns banks and by the fresher coastal source waters from along the Scotian Shelf (Townsend, 1991; Pringle, 2006). The GOM has an average tidal range > 3.0 m and experiences large seasonal amplitudes in surface salinity (Geyer et al., 2004), net primary productivity (O’Reilly et al., 1987), sea surface temperature (SST) and pCO₂ (Vandemark et al., 2011). The key circulation feature influencing the region is the Maine Coastal Current, which flows counterclockwise and delivers freshwater and constituents from the northeast GOM to the southwest GOM (Pettigrew et al., 2005).

Physical processes that deliver and mix the nutrient pools supporting the GOM ecosystem (e.g., strong tides, wind-driven mixing, large annual ranges in temperature and salinity) also generate thermodynamic variability in carbonate system parameters at diurnal through annual scales. Further, such processes contribute to an intense annual cycle of CO₂, whereby disequilibrium with the atmosphere is partially balanced by an air–sea flux of CO₂ (Shadwick et al., 2011; Vandemark et al., 2011; Wang et al., 2017). The degree to which physical variability can modify OA in the GOM is substantial (Salisbury and Jönsson, 2018). In a climatological study using data from 1950 to 2013, the GOM recorded an intraannual SST range of 5 °C–15.5 °C and salinity range of 31.5–32.5 (Richaud et al., 2016). The annual temperature range alone elicits a significant change in the carbonate system. For example, assuming the approximate mean GOM salinity (32.2), the mean TA (2,184 µmol kg^–1) and an atmospherically equilibrated seawater surface pCO₂ (400 uatm), temperature variability causes an annual increase of 0.013 in pH_TOT (pH total scale) and 1.06 in Ω_a.

Within the GOM, the period spanning 2003–present (the period for which moored data exist via the Northeast Regional Association of Coastal Ocean Observing Systems) has been characterized by temporally coherent oceanographic phases that last on the order of 1–3 years. These phases result from periods of primarily warm slope-water intrusions into deeper layers of the GOM, bringing warm and salty high-nutrient waters, followed by episodes of increased shelf-water fluxes from the Scotian Shelf that bring colder and fresher low-nutrient waters to the GOM (e.g., Townsend et al., 2014, 2015). Beginning in about 2008, the GOM entered an extended warm and salty phase, which was an important contributor to the anomalously warm surface waters reported by Pershing et al. (2015). However, in the last 7 years, since 2013, greater short-term variability has been observed, with some of the coldest (and freshest) and warmest (and saltiest) episodes in the more than 15-year moored record. The timing of these water mass fluxes and associated surface temperature changes maintained by atmospheric heat flux (Chen et al., 2014, 2018) have combined to produce short-term, rapid warming of surface waters in the Gulf (on the order of 0.15 °C y^–1) as reported by Pershing et al. (2015). Such water mass fluctuations may be related to changes in the Labrador Current (Townsend et al., 2010, 2015; Mountain, 2012; Claret et al., 2018) and a northerly shift of the Gulf Stream as described in Saba et al. (2016; also see Grodsky et al., 2017). Although the majority of the warming trend over the period 1982–2018 can be attributed to natural variability, a third of the trend has been attributed to anthropogenic warming (Chen et al., 2020).

Seasonal variability in the carbonate system of the water column in the GOM is affected mostly by and synchronized largely with the seasonal progression of stratification overturn, primary production, respiration–remineralization, and water mixing (Wang et al., 2017). Different processes may alternate to exert the main control on the variability at certain times of the year (Wang et al., 2017). In the winter, saturation states at the surface are at their lowest (Figure 1) and the water column is well mixed. Subsequently, spring blooms draw down carbon at the surface subsequently raising the Ω_a, while at the same time, the freshet arrives bringing stratification and low Ω_a waters to the surface. Over the summer, surface waters warm and stratification continues to isolate the bottom waters, while circulation changes act to retain local organic material in subsurface waters where remineralization of that material draws down saturation states. The fall brings increased storms, precipitation, mixing, and a decline in saturation states near the surface as seasonal trends cool the surface and reduce productivity. Overall, photosynthesis–respiration is the main driver controlling the seasonal variability of Ω_a (Wang et al., 2017). Less is known about subsurface seasonal patterns. Interannual variability in the system is primarily governed by the magnitude of ocean productivity and by the relative supply of differing water sources from regional currents. Waters sourced from the Labrador Current and Gulf of St. Lawrence have lower carbonate ion concentrations (and thus lower Ω_a) versus those supplied by the higher salinity Gulf Stream. Within the GOM, these waters are subsequently modified by remineralization and respiration as well as riverine inputs (Wanninkhof et al., 2015; Wang et al., 2017).

Many marine species are likely to be affected by projected OA conditions in 2050. Evidence supporting the effects of OA on commercially important marine bivalves is clear (Gazeau et al., 2013; Waldbusser and Salisbury, 2014), and more information continuously becomes available for other groups such as crustaceans and fish (e.g., Dixson et al., 2010; Keppel et al., 2012; Long et al., 2016a; Rodriguez-Dominguez et al., 2018). These data indicate that acute or chronic responses are observed between Ω_a of 1.0 and 2.0. The response is species-dependent, but Ω_a = 1.5 is a conservative estimate of when physiological effects have been observed and as such will be referred to as suboptimal conditions for marine organisms (i.e., Ω_a < 1.5). The physiological effects of OA on life stages and populations of particular species are typically better understood than the ecological consequences at ecosystem levels. Similarly, the short-term effects of Ω_a, pH, and buffering capacity on marine organisms have been better characterized than their long-term effects. Thus, a unified understanding of long-term combined effects of OA with accompanying changes in temperature, circulation patterns, freshwater inputs, and productivity, among other factors that are also projected to change in the region, is still lacking.

In the GOM, there has been progress in understanding the effects of carbonate chemistry in key fisheries and aquaculture industries, such as those that target the American lobster, Eastern oyster, and Atlantic sea scallop (reviewed by Gledhill et al., 2015). However, this progress falls short when considering the number of shell-forming species with high and very high climate vulnerability on the Northeast U.S. continental shelf (Hare et al., 2016), their complex life cycles, and the interactions of environmental factors. Table 1 of Gledhill et al. (2015) summarizes the acidification studies of seven commercially important species of the New England/Nova Scotia region, including the American lobster, four species of bivalves, the summer flounder, and the Atlantic longfin squid. Studies on approximately 140 additional species living in or related to species living in the region are also available (Table S1 in Gledhill et al., 2015). Here, we included more than 60 studies, conducted with species from the GOM since 2015 and focused on commercially important species (Tables S1 and S2). Below, we discuss the highlights from these recent studies in the context of prior results discussed in previous efforts (e.g., Gledhill et al., 2015; Hare et al. 2016).

Fisheries in the GOM rendered over US$1.2 billion in 2017 in the United States alone, with 74% of the total landings composed of crustaceans and bivalves. These two groups have been identified as highly susceptible to OA and temperature increases, and both have complicated life histories with their larval stages primarily pelagic and adult stages mostly benthic. American lobster is currently the highest value fishery in the United States and Canada and the top species landed by weight and revenue in the region, with US$551 million harvested from the U.S. GOM in 2017 (NOAA Fisheries, 2019), 1.42 billion CDN in Atlantic Canada in 2018, with New Brunswick and Nova Scotia contributing a landed value of 1.06 billion CDN, of which US$730 million were harvested from Nova Scotia (Department of Fisheries and Oceans, 2018, Zonal Interchange File [database], Ottawa, accessed August 4, 2020). Although American lobsters have been observed in waters as deep as 800 m, the majority of landings occur from waters shallower than 80 m in the GOM (Tanaka and Chen, 2016; Oppenheim et al., 2019). Warming effects in lobsters have been studied extensively (e.g., Rheuban et al., 2017; Jaini et al., 2018; Boavida-Portugal et al., 2018). In Table S1, we included recent studies examining the effects of OA alone or both OA and warming in different life stages, often at experimental CO₂ levels higher than projected for 2050. Larval studies up to Stage IV of development found that warming had a greater adverse effect than increased pCO₂ alone (Waller et al., 2017) and when combined, increased temperature and pCO₂ altered dry mass, carapace length, swimming speed, feeding rates, and survival of lobsters. On the other hand, experiments with older larvae (Stage V) to juveniles found increases in mortality, slower development, and increases in aerobic capacity under increased pCO₂ alone (Menu-Courey et al., 2019). Similarly, McLean et al. (2018) found that increasing pCO₂ resulted in decreases in length and biomass, longer intermolting periods, and more susceptibility to shell disease in juvenile lobsters. Juvenile lobsters also showed immunosuppression and reduced cardiac performance when exposed to increased pCO₂ and acute thermal stress for 60 days (Harrington and Hamlin, 2019). These recent studies further support detrimental physiological responses to increased OA for different life stages of the American lobster, which will be exacerbated by accompanying warming conditions.

In GOM, bivalves accounted for 36% of the total U.S. landings revenues in 2017, with the Atlantic sea scallop yielding US$341 million (28%) followed by the eastern oyster yielding US$35 million (3%). In Canada, the sea scallop fisheries for Nova Scotia and New Brunswick combined yielded US$164 million CDN, while eastern oyster aquaculture for the same regions yielded US$20 million CDN (Statistics Canada. Table 32-10-0107-01 Aquaculture, production and value. Retrieved November 22, 2019). The majority of their adult benthic habitat occupies the nearshore and shallow shelf regions of the GOM, except for Atlantic sea scallops that reside further offshore with a regional hot spot around Georges Bank. Recent integrated assessment models for the Atlantic sea scallop (Cooley et al., 2015; Rheuban et al., 2018) suggest that under different CO₂ emission scenarios and different impact and management conditions, sea scallop biomass could be reduced between 13% and 50% by the end of the century (Table S1). These projections and the high revenue underscore an urgent need for laboratory and field research on how sea scallops respond to OA and increased water temperature, so that the industry and managers can respond accordingly to maintain sustainable stocks. Similarly, the Atlantic surfclam is another commercially prized species, with ample distribution from Canada’s Gulf of St. Lawrence to North Carolina in the United States, that has already shifted northward in its distribution (Timbs et al., 2019), but studies of its response to OA have begun only recently (Pousse et al., 2020; Meseck et al., 2021; Table S1).

In comparison with scallops and surfclams, eastern oysters are better studied, and their response to OA varies across life stages, with early stages typically being more susceptible than adults (Table S1). Other bivalve species such as soft-shell clams and mussels have also shown this pattern (Tables S1 and S2). Although there are more studies on early life stages, a few recent studies of the eastern oyster Crassostrea virginica have examined the short-term effects of OA on adults and found minimal change in length or in gaping behavior (Table S1). Both studies of adults in Table S1 involved experiments of relatively short duration, and calcite in adult shells has less corrosivity than aragonite in larval shells in general. Although gaping behavior in adult oysters did not seem affected, research on other adult bivalves, such as blue mussels, bay scallops, and hard clams, has shown detrimental effects with reduced growth rates and tissue weights when exposed to OA (Table S1). Less frequent field experiments also support adverse OA effects in bivalves as evidenced by reductions in settlement rates when low pH conditions prevail in sediments for Mya arenaria and Nucula spp (Table S1). An increased repertoire of studies, including more species, different life stages, multiple concurrent stressors, field and lab experiments, modeling projections, transgenerational observations (Griffith and Gobler, 2017), and potential for mitigation in multispecies studies (Young and Gobler, 2018), have considerably improved the understanding of potential responses to OA for individual bivalve species in the GOM.

Several studies from the region have also highlighted the complex response of finfish to OA (Tables S1 and S2). Atlantic sea herring (Clupea harengus) is an important commercial and keystone species at the base of the food web in marine ecosystems in the GOM. Studies on the effects of OA on swimming behavior, energetic costs, and survivorship of this species have yielded contrasting results (Table S1). This variability suggests that different populations in the region might vary in their response to OA, which should be explored further (Leo et al., 2018). Similar to results from bivalve species above, studies on the forage fish Atlantic silverside (Menidia menidia) and inland silverside (M. beryllina) have highlighted differences between experiments where OA was a unique factor and experiments that included additional stressors such as high temperature (Table S2).

Experimental studies in the GOM have also illustrated contrasting responses of different physiological parameters within a life stage for a single species. For example, Stiasny et al. (2018) observed that cod larvae were larger under energy limitation in the elevated CO₂ treatment compared to the ambient CO₂ treatment, but larvae under elevated CO₂ had developmental impairments in critical organs, such as the liver. All these studies make clear that our capacity to replicate environmental and biological complexity in experimental studies is limited, yet more sophisticated and longer term studies are critical to enable further understanding of physiological mechanisms, energetic balances in organisms, and population dynamic consequences of OA in marine species. A handful of other studies, including responses through trophic webs (e.g., Jin et al., 2015) and other ecological interactions (e.g., Fay et al., 2017; Schulz et al., 2017; Gilbert, 2020; Raven et al., 2020; Tester et al., 2020; Wilson et al., 2020), have suggested potentially larger but varied ecological and ecosystem consequences of OA in the GOM. Thus, a better integration of our current knowledge, predictive capability, and collaboration will be key to addressing management and societal changes that ameliorate OA consequences in the region and preserve the sustainability of marine resources and the services they provide to the GOM.

In this section, we describe the methods employed to compute regional trends from observations, the evaluation methods we relied on for the regional models of the area, and the approach we used for the entire suite of regional projections produced for this work. A combination of regional high-resolution simulations that include coastal processes was used to project OA conditions for the GOM into 2050. We also discuss the approach for each model and include some discussion of the global models as well.

Several different data sets were used for this work including observations form Buoy D, the CML site, and some regional transects including those within WB. Surface observations from the moored NOAA PMEL/UNH Coastal CO₂ Buoy D, as described in Sutton et al. (2019), were used in this work. The Buoy D data set spans the period 2006–2018. The CML data set spans 2008–2019. Buoy temperature and salinity were measured by a Sea-Bird Electronics SBE-16plus (Sea-bird electronics, Bellevue WA, salinity and temperature accuracies of ±0.0005 Siemens (S) m^–1 and ±0.005 °C, respectively), while those from CML were measured with an Aanderaa 4319B (Xylem, Rye Brook NY, salinity and temperature accuracies of ±0.0018 S m^–1 and ±0.05 °C, respectively). PMEL/UNH buoy CO₂ measurements were from an MAPCO2 system (Sutton et al., 2019), while those from the CML site were made using a showerhead-type equilibrator with headspace air recirculated through a Li-cor LI840 nondispersive infrared detector (details of this system in Hunt et al., 2017, and references therein). Each CO₂ system was calibrated with a CO₂-free gas stream and a CO₂/nitrogen mixture that was cross-calibrated to standard gases from NOAA/Earth System Research Laboratories. The uncertainty of the PMEL/UNH CO₂ measurements was better than ±3 µatm, while that of the CML CO₂ measurements was estimated as ±3 µatm (Vandemark et al., 2011). Calculations were performed using the CO2SYS program (van Heuven et al., 2011) with K1 and K2 constants from Millero (2010), KSO₄ from Dickson (1990a, 1990b), and KB from Uppstrom (1974). Ω_a at both the PMEL/UNH Coastal CO₂ Buoy and UNH CML sites (Figure 2) was calculated from in situ measurements of salinity, water temperature, pCO₂, and TA modeled from salinity and water temperature using an empirical model (McGarry et al., 2021). CML monthly salinity and Ω_a climatologies (Figure 2) were constructed from mean monthly values, which were interpolated to a daily time step (“pchip,” MATLAB, MathWorks Natick MA) and smoothed with a 60-day moving mean (“smooth,” MATLAB).

Monthly-to-seasonal cross-shore transects have been conducted by the University of New Hampshire since 2004. The transect line runs more than 80 km west to east, ending at station WB-7 (bottom depth approximately 250 m) in WB, one of the three deeper basins in the GOM (Figure 2). Discrete samples for lab analyses are collected from Niskin bottles closed at depth intervals through the water column. Ω_a at the WB-7 site was calculated from in situ measurements of salinity and water temperature (from a SBE-16plus) and bottle measurements of DIC and TA. DIC of unfiltered water was measured using an Apollo SciTech AS-C2 analyzer (Newark DE), while TA of unfiltered water was measured via Gran titration using an Apollo SciTech AS-ALK2. DIC and TA method precisions were ±2 µmol kg^–1. Further details of DIC and TA sample collection and lab analysis methods can be found in Hunt et al. (2013).

The water column at WB-7 was divided into three depth strata for trend analysis, based on sampling consistency over the program period, with the aim to isolate the surface, mixed layer, and bottom conditions. Linear trend lines were calculated for each depth stratum by using deseasonalized data based on the methods described in Takahashi et al. (2009). For comparison, the data collected from similar locations as WB-7 by Woods Hole Oceanographic Institution (WHOI) (Wang et al., 2017), NOAA Gulf of Mexico and East Coast Carbon cruises 1&2 (Wang et al., 2013; Wanninkhof et al., 2015), and East Coast Ocean Acidification 1 were also included in the analysis.

The Coupled Model Intercomparison Project (CMIP) considers various representative concentration pathways (RCPs) that are primarily based on human consumption of fossil fuels (Bopp et al., 2013). Four pathways, namely, RCP2.6, RCP4.5, RCP6, and RCP8.5, are identified according to a possible range of radiative forcing values in the year 2100 (2.6, 4.5, 6.0, and 8.5 Wm^–2, respectively). The models used for CMIP5 and CMIP6 coarsely resolve the ocean (approximately 1°) and consequently neglect important coastal processes that can modify global rates of acidification locally in regions like the GOM. However, global models can be used to drive regional higher resolution models designed to resolve coastal processes important for constraining regional rates of acidification and provide important information on how boundary conditions are likely to change. The regional projections used for this work are described in Table 1.

Lavoie et al. (2020) projected ocean conditions using a coupled physical–biogeochemical regional climate model, developed for the Gulf of St. Lawrence and the Scotian shelf called the Gulf of St. Lawrence Biogeochemical Model (GSBM; Gehlen et al., 2015), that was driven by three global projections from CMIP5 for the RCP 8.5 scenario. The three ensemble members include the second generation Canadian Earth System Model (CanESM2), the low-resolution version of the Max-Planck-Institut Earth System Model (MPI-ESM-LR), and the latest Hadley Center Earth System climate model (HadGEM2-ES) and simulate the period 1970–2100. Future boundary conditions were built by adding the simulated trends from the ESMs to the observed climatology for the physical variables (see Long et al., 2016b, for temperature and salinity), as well as for nitrate, DIC, TA, and dissolved oxygen. The trends from the ESMs that were used to build future boundary conditions for the regional model are generally consistent (Lavoie et al., 2019). The regional model includes the GOM, but the results for that region had not been analyzed until now. Because this model includes biogeochemistry, the projections are referred to as the dynamically downscaled biogeochemical projections.

Other dynamically downscaled projections exist for the GOM as described in Brickman et al. (2021), but most do not include biogeochemistry. A set of simulations were performed using the Regional Ocean Modeling System (ROMS; Shchepetkin and McWilliams, 2005), referred to as ROMS throughout this work, each forced by a global climate model. The model domain covers the entire Northwest Atlantic, including GOM shelf, slope, and the Gulf Stream. The grid is configured with horizontal spacing of approximately 7 km and 40 vertical terrain-following levels (Kang and Curchitser, 2013). A control (CTRL) ROMS simulation was performed using data-assimilated fields at the surface and along the side boundaries in the ocean over the period 1976–2005. The ROMS climate change simulations were performed using the delta method. The deltas were computed by subtracting the long-term monthly mean values during the period 1976–2005 from those during 2070–2099 in each of the three general circulation models: NOAA Geophysical Fluid Dynamics Lab (GFDL), Institut Pierre Simon Laplace (IPSL), and HadGEM. Then, the deltas were added to the CTRL providing three independent future simulations using RCP 8.5 (Alexander et al., 2020). Because the ROMS simulations used 2070–2099 for the future period, the results needed to be adjusted to estimate the changes in 2050. We tested several methods to accomplish this goal using the time-varying CMIP5 simulations, which include values in 2050, to evaluate the accuracy of the adjusted values. A method that worked well uses the global average change in radiative forcing from the RCP8.5 scenario. The change in forcing from 1976–2005 to 2050 is 3.1 Wm^–2 and to 2085 is 5.68 Wm^–2, so the grid point values of temperature and salinity were scaled by 0.546 (3.1/5.68) to estimate the changes that occur by 2050. These results are presented in more detail in Brickman et al. (2021).

To include biogeochemical fields in the projected conditions without dynamic biogeochemistry, an empirical model for the carbonate system to the future conditions was applied and CO₂ was added to the DIC pool using the atmospheric concentration in 2050 in the RCP scenario, assuming the anthropogenic signal is well mixed throughout the water column. This approach has the benefit that relative contributions of physical changes, like warming, are easily decipherable. An empirical model for each carbonate system parameter was developed using a least-squares multiple linear regression (MLR) of the carbonate system parameter on hydrographic variable observations collected on the northeast U.S. shelf between 2007 and 2015 (McGarry et al., 2021). Although empirical models including other variables were generated, the models used here rely solely on temperature and salinity but skillfully reproduce observed DIC and TA (McGarry et al., 2021; their Equations I and IV). Ω_a was computed from the empirically derived DIC and TA using CO2SYS (Lewis and Wallace, 1998). Additional model evaluation was performed as part of this work, as described below.

Anthropogenic CO₂ concentrations were added to the DIC computed from the MLR by first calculating the seawater pCO₂ in equilibrium with the atmospheric concentration in 2050 from RCP 8.5 (550 ppm) using the solubility equations from Weiss et al. (1970) for the projected temperature and salinity. The DIC in equilibrium with future atmospheric concentrations was then computed using the equilibrium pCO₂ and MLR-computed TA as inputs in CO2SYS (Lewis and Wallace, 1998). The same method was used to calculate DIC for the control conditions and the difference was added to the projected DIC from the MLR. Note that this method will only be representative of the conditions experienced by the training data for the MLR. Because the dynamics of oxygen and carbon are linked, the inclusion of oxygen in the MLR improves the representation of the carbonate system. However, this variable was not available for regional models used here, so the MLR was based on temperature and salinity only. The results of these efforts are reported below and referred to as dynamically downscaled physical projections with empirical biogeochemistry.

The moored observations described in Salisbury and Jonsson (2018) from Buoy D were used to compare to the base and future projections of the downscaled model. Statistical analyses were used to quantify the comparisons.

Duration is defined as the number of months each year that Ω_a falls below the critical threshold Ω_a,_β (1.5). Intensity represents the integrated annual acidification stress defined as $I=∫t=1t=12$(Ω_a,_β − Ω_a,t)dt, where I is the intensity, Ω_a,_β is the threshold (1.5), and Ω_a,t is the monthly modeled value when Ω_a,t < Ω_a,_β (i.e., values > Ω_a,_βv are disregarded). This acidification intensity score is similar to other threshold-based stress indices (e.g., degree heating week for gauging coral thermal stress) that adopt the assumption that low-level stress experienced over extended periods can develop into a higher level accumulated stress response similar to that experienced in a short-term chronic or acute stress response. The benefit of the intensity score is that it combines both the magnitude by which the threshold is exceeded as well as the duration.

To determine the role of warming and circulation changes in the future projections of OA, a “slug” of anthropogenic CO₂ was added to the modern/base case DIC and compared to the projected fields for the ROMS empirical projections. Thus, it yields a modeled field assuming only an increase in DIC in response to OA independent of other changes.

The UNH repeated transect measurements off the NH coast during the period of 2005–2017 showed that subsurface Ω_a from depths of 18–34 m exhibited trends which closely followed those of surface water shallower than 10 m at the deepest station in the WB of the GOM (Figure 3). Bottom water Ω_a deeper than 210 m at the same station (WB7, Figure 3) was less correlated with the surface distributions. No clear seasonal pattern of Ω_a was observed for either the surface or subsurface offshore waters, which is consistent with the finding by Wang et al. (2017). Inconsistent sampling frequency of these data might have contributed to the lack of observed seasonality or missed capturing additional shifts or events affecting aragonite saturations, particularly as GOM studies have shown that low Ω_a or high CO₂ events in both surface and bottom waters can dissipate within days or weeks (Waldbusser and Salisbury, 2014). Surface and subsurface waters experience occasional low Ω_a conditions at this offshore location below the threshold value (Ω_a = 1.5; red line in Figure 3). The observed Ω_a values near bottom were persistently below this threshold with only a few exceptions during the 12-year period. A few shorter data sets collected by WHOI and NOAA Ocean Acidification Program generally agree with these patterns (Figure 3; Wang et al., 2017). A long-term trend analysis was conducted by subtracting the respective monthly climatological mean Ω_a from individual Ω_a values in Figure 3 based on the method of Takahashi et al. (2009). The result shows a positive trend in ΔΩ_a for all three layers of the water column in WB, with an annual increasing rate of 0.02–0.04 per year (Figure 3). The rate of Ω_a increase at the surface was about twice that of the near-bottom water mass. These trends coincide with (nonsignificant) positive trends in salinity and temperature at the same site (the same colocated data as in Figure 3; Figure S2), as well as the Ω_a rate estimated for the GOM for the period 2004–2015 (0.0152 ± 0.003; Salisbury and Jonsson, 2018). These observations confirm that the recent warming and observed salinity increase over the last decade, likely due to the increase of slope water intrusion (Salisbury and Jonsson, 2018), may play an important role in controlling carbonate chemistry in the GOM.

The moored surface observations shown in Figure 1 from the PMEL/UNH Coastal CO₂ Buoy D, as described in Salisbury and Jonsson (2018), were compared to the downscaled modeled base and future projections of Ω_a (Figure 4). The monthly mean values for the entire time series (2006–2018) are shown in black squares and the standard deviation of the observations each month is provided, with additional evaluation provided in Figure S1. The observations indicate that Ω_a experiences strong seasonality at this location, with higher values at the surface during summer and lower values during winter, which can also be seen in Figure 1. Presently, surface conditions away from the shore rarely experience suboptimal conditions for marine organisms and rarely decline below the saturation horizon (Figures 1 and 4).

Both of the downscaled models used for the projected conditions in this work perform well when compared to the historical observations of Ω_a (GSBM: R² = 0.75; ROMS: R² = 0.86; Figure 4). GSBM underestimates the seasonality of the observations, with slightly higher saturation states in the winter than observed, although that is the observed trend. In the spring, ROMS with empirical biogeochemistry likely underestimates the role of the spring bloom in raising the Ω_a by scavenging DIC from the surface mixed layer, and so the model underestimates the saturation state at that time. As a result, winter and spring projections are more uncertain as they are not as well simulated in the model suite. Despite these minor inconsistencies, the ability of the models to capture the historical climatological conditions suggests that they are robust enough to reasonably project future conditions.

The empirical biogeochemical model can also be applied to the GSBM simulation to understand the relative importance of the physical simulation in GSBM to the biogeochemical feedbacks. The GSBM model control simulation produces a surface average Ω_a for the GOM of 1.68 and 1.48 at the bottom. When the empirical biogeochemistry is applied to the GSBM temperature and salinity fields, the mean surface Ω_a for the control simulation is 1.84, and the bottom is 1.63. Both of the results from the empirical biogeochemical model produce higher saturation states than the dynamic biogeochemical model.

The CMIP5 ensemble average under RCP 8.5 projects a decline in Ω_a of 0.44 at the surface (not shown) for the entire GOM region when the control values (for the period 1976–2005) are subtracted from the average for the future (2005–2055). When the MLR is applied to the temperature and salinity from the CMIP5 future projection, the projected decline is reduced (0.17). In the subsurface, the global models only simulate bottom conditions at the shallowest around 300-m depth and no longer serve the DIC and TA values for the region, so the MLR method is used. Using that depth as representative of the bottom conditions projected for the GOM, the CMIP5 ensemble average for the GOM projects a decline in Ω_a of 0.15.

Surface and bottom Ω_a decline in the future in the GSBM projections for 2050 under RCP 8.5. The decline at the surface is projected to be 0.40 for the MPI-ESM-LR and CanESM2 forced projections and 0.42 for the HadGEM2-ES ensemble member (Table 2). The ensemble members largely agree on the magnitude and direction of the trend. Bottom conditions are largely consistent with the surface, with a slightly weaker decline (0.36–0.37; Table 2). Although the surface is comparable to the global model projected change, the subsurface change is larger than the global models, in part due to the lack of resolution over these depth ranges in the global models. In addition to the future projected changes, the absolute values for Ω_a in the global model are higher than in the downscaled models predicted in both the control (C) and future (F) simulated conditions, that is, C:2.73 and F:2.29 versus C:2.00 and F:1.68.

The decline in Ω_a is driven in the regional model mainly by the increase in DIC in the region, which occurs at all depths and in all three simulations (Figure 5). We note that increasing DIC is mainly attributable to increasing atmospheric CO₂ but may also be imported into the region or result from increasing rates of respiration with warming conditions. The DIC increase is greatest at depth. The decrease in pH in the 150–300 m layer is greater in the GOM than in the Scotian Shelf or GSL (Lavoie et al., 2020) because the DIC increase is greater in the GOM. The model suggests that this difference is caused by an increase in primary production, mainly by diatoms, and greater mesozooplankton mortality in summer due to higher water temperatures, with both leading to an increase in local respiration of organic matter at depth in the region. Stratification is also stronger in the projection which reduces the ventilation of subsurface waters.

Alkalinity generally decreases at the surface and increases in all the regional simulations at depth (150–300 m; Figure 6). The surface trends differ widely among the three ensemble members in magnitude with HadGEM2-ES producing the steepest decline in TA, and due to this disagreement, the trends are highly uncertain at the surface and midwater column (Figure 6A and B). These changes in alkalinity primarily correspond to the changes in salinity in the region with fresher surface waters and saltier deep waters.

Because the results were simulated continuously between 1970 and 2100, the simulations project the timing when key thresholds are exceeded in the region. Like pH, Ω continues to decline under RCP 8.5. For key thresholds, like Ω_a and saturation state for calcite, the thermodynamic threshold of 1 (annual mean) is reached in the benthic environment in all three ensemble members by the year 2060 in the southern GOM for aragonite. By 2081, calcite remains above saturation for the entire domain. Suboptimal Ω_a conditions of less than 1.5 are projected for the surface by 2040–2060 with some uncertainty between the simulations and already present in the benthic environment by 2020 for aragonite and by 2060 for calcite in the southern portion of the GOM (Figure 7). Georges Bank is a region that remains above the threshold for the longest period of time, and one (CanESM2) ensemble member entirely avoids the threshold at the surface. As such, Georges Bank is one of the regions of largest disagreement between the ensemble members. The Scotian Shelf is one of the regions that experiences the earliest threshold crossing in the domain at both the surface and the bottom.

Consistent with the dynamical biogeochemical model results from GSBM, the surface and bottom Ω_a conditions are also projected to decline with the empirical biogeochemistry results (Table 3). The surface Ω_a decline is projected to be 0.12 and 0.14, respectively, from the GFDL and IPSL forced ensemble members and 0.03 for the HadGEM simulation. The bottom Ω_a decline is projected to be 0.11 and 0.15, respectively, for the GFDL and IPSL forced ensemble members and 0.02 for the HadGEM simulation. These changes are smaller than projected by the dynamical biogeochemistry GSBM model and not as similar with respect to depth. However, the empirical biogeochemical CTRL conditions are lower than the GSBM model. When the empirical model is applied to the GSBM temperature and salinity values for the projection, the projected decline in surface Ω_a for the mean change across all three ensemble members in the GOM is 0.37 and, in the bottom, 0.20, suggesting that even in those simulations, the empirical model does not support uniform change vertically in the future. Although the surface is comparable to the global model MLR-based projected change (0.17) and lower than the decline projected by the DIC and TA-based projections (0.44), the subsurface declines less severely for two ensemble members (0.02 and 0.11; Table 3) or consistently for one ensemble member (0.15; Table 3) when compared to the global model projected change (0.15; Table 3). In addition, the absolute values for Ω_a in the global model (C:2.04, F:1.87) are higher than in the downscaled models (C:1.75, F:1.64) predicted in both the control and future simulated conditions.

Under RCP 8.5 in 2050, the ROMS model ensemble projects that persistently suboptimal conditions for marine organisms (i.e., Ω_a <1.5) will exist at the surface for the first 5 months of the year (January–May; Figure 4). Maximum supersaturation will shift from late summer (August) to early fall (September). One of the three ensemble members (HadGEM) does not project as severe a seasonal evolution, with fewer months experiencing suboptimal conditions for marine organisms (i.e., Ω_a < 1.5), because it is the projection with the warmest temperatures which continue to counteract the effects of OA. This shift in the presence of undersaturated waters in the spring and the timing of maximum saturation state later in the year may have biological impacts, as these seasonal transitions often correspond to reproductive and early life stage transitions in many marine organisms, which are typically more vulnerable to OA.

Projected conditions for 2050 are very different from control conditions for RCP 8.5 in terms of DIC for all three ensemble members (Figure 8). However, in all cases, DIC anomalies are highest offshore and lowest near the coast. Although the magnitude varies across ensemble members, the fall produces the sharpest gradients in the anomaly between the coast and offshore and experiences the most variability between ensemble members. In all three ensemble members, using the empirical model alone without the adjustment for the higher atmospheric CO₂ concentrations or the “slug,” DIC concentrations are less than modern conditions (not shown), suggesting that physical changes to the GOM along with biological feedbacks are acting to mollify the increases from the atmospheric change.

All three ROMS simulations experienced lower TA at the surface (Figure 9), consistent with the projected salinity and dynamically downscaled biogeochemical projections described above. Subsurface TA increases in the future, consistent with the GSBM results. Seasonally, the late fall conditions (October–December) experienced the greatest magnitude of change, although the three simulations vary widely in the magnitude of the changes.

Biologically relevant impacts from OA, particularly for shell-forming organisms, can be interpreted from the duration of exposure to conditions below the identified Ω_a threshold of 1.5 (i.e., number of months when Ω_a is < 1.5) as well as the intensity of exposure (magnitude below Ω_a 1.5). Currently, the maximum duration of exposure in the modern baseline model at the surface does not exceed 4 months (“CTRL” in Figure 10) and is concentrated nearshore, with some of the southern GOM region further offshore experiencing approximately a month of exposure. In the future projections for 2050 under RCP 8.5, the maximum duration of nonoptimal exposure at the surface continues to be in the nearshore but is longer (6–8 months) and expands spatially to include the majority of GOM surface waters, with approximately 4 months of exposure. The ensemble members project a range of conditions, with the IPSL projection providing the most severe duration and HadGEM projecting the smallest duration encompassing the smallest subsurface area.

At the bottom, baseline control conditions are more severe than at the surface, and the majority of the GOM experiences below-threshold conditions for nearly half the year with the WB experiencing year-round exposure. In the 2050 projections, the duration of below-threshold conditions increases to year-round over the majority of the region, with the exception of Georges Bank where it increases from 4 months year^–1 in the control to 8 months year^–1 (Figure 10).

Intensity of exposure at the surface is tightly coupled to the nearshore coastline in the modern baseline (“CTRL” in Figure 11). In the future, the intensity of exposure becomes more severe along the coastline and spreads over the GOM with a wide range projected among the ensemble members. The disagreement between the projections stems from the degree of warming simulated in each of the projections, with HadGEM projecting the most intense warming and thus the least severe intensity of Ω_a < 1.5. Georges Bank remains one of the most favorable environments in the GOM region, remaining largely above the threshold. The difference between the “slug” and the future ensemble conditions illustrates the partial compensation that continued warming and circulation changes impart on future OA conditions, which partly limit acidification intensity at the surface and even more so at depth (Figure 11).

Results indicate that warming continues to partly reduce the impacts of OA in the future as demonstrated by the “slug” compared to the projected change for the region (Figures 4, 10, and 11). Under the RCP 8.5 scenario, both intensity and duration of below-threshold conditions would be more severe in the GOM in 2050 without the additional warming and circulation changes (Figure S2). The duration doubles without the amelioration of circulation and warming, and the intensity score (I) is a full unit more severe. The slug map indicates that the warming and circulation act to mollify the OA duration, particularly in the nearshore surface environment and more so than at the bottom (Figure 10).

Regional high-resolution downscaled projections for the GOM of Ω_a conditions inclusive of coastal processes under RCP 8.5 project that the entire region is likely to experience Ω_a < 1.5 by 2050. These projections are in contrast to the CMIP5 projections for the region, particularly in the subsurface environment, which predict a decline in future Ω_a, but the projected values do not fall below the Ω_a = 1.5 threshold over much of the region by 2050. Subsurface waters are not well represented by coarse global models, yet house many important shellfish species in the region. In the discussion that follows, we first put the projected results in the context of prior projections for the region. We compare the models used in this work and suggest which processes are vital to best reduce uncertainties in model projections moving forward. We also consider the implications of the results for the vulnerable ecosystems of the GOM, with a particular focus on the nearshore and subsurface environments and economically important species. Finally, we summarize adaptive societal measures that could be employed in the region to potentially avert or mitigate these future changes.

Globally, the oceans simulated using CMIP5 under RCP 8.5 are projected to be higher in pCO₂ and more acidified (surface pH decrease of –0.33 ± 0.003; Gattuso et al., 2015). Regionally, the direction of the trend predicted here is the same, but there are differences between the global and downscaled results, mostly on the subsurface shelf which is home to many valuable shellfisheries in the region including scallops. Prior projections used to assess OA impacts to the scallop fishery in the region relied on the global model projections from CMIP5 and did not report Ω_a (Rheuban et al., 2018), but the previous iteration of that study by Cooley et al. (2015) reported deep box Ω_a values that were greater than those reported here. Our summary of CMIP5 projections for the region suggests that the benthic habitat where adult shellfish reside will endure more severe change in the future than previously estimated.

Global coarse resolution model projections differ from the regional projections discussed here because the intrusion of anthropogenic CO₂ is not the only mechanism that can reduce Ω_a within coastal surface waters. Local processes like freshwater delivery, eutrophication, productivity, and sediment interactions that drive variability on regional and subseasonal scales can also modify spatial variability in Ω_a (Feely et al., 2008, 2018; Cai et al., 2011; Qi et al., 2017; Siedlecki et al., 2017; Pilcher et al., 2018). Although the downscaled simulations presented here do not simulate the highly variable and more severe nearshore Ω_a conditions depicted by the CML observations in Figure 1, the shelf variability is well captured by the models. The downscaled surface trends are comparable to the global model projected change, but the subsurface change is larger than projected by the global models, suggesting that there are processes in the downscaled simulations amplifying this rate of change within this critical habitat for scallops, lobsters, and other vulnerable benthic dwellers. Such processes could include those previously reported (e.g., water column metabolism and sediment interactions), but additional research is needed to understand the role of these processes on decadal timescales. Additionally, existing projections, including the ones presented here, should be used as conservative estimates for nearshore environments as they already experience suboptimal conditions for marine organisms (Figure 1) more severely than regions offshore in the GOM.

Differences exist between the regional projections using dynamical biogeochemistry and the empirical biogeochemistry, and those differences can be exploited to prioritize observation needs and process-based understanding of potential biogeochemical feedbacks for future work. The empirical model of biogeochemistry assumes the net biological and sedimentary flux contributions at present will be similar in the future. Differences between the empirical and dynamic biogeochemical model projections suggest that the balance of those processes may change and result in biogeochemical feedbacks in the future. Here, the projections differ in their vertical structure for future change: The surface experiences the same amount of change, but the subsurface is projected to decline more in the dynamical biogeochemistry in conjunction with the surface. This effect likely results from a biogeochemical feedback but could also correspond to altered residence times due to circulation changes in the region relating to weakening of the Gulf Stream, as is further described in Alexander et al. (2020). The ROMS projections used here and in Alexander et al. (2020) show an enhanced counterclockwise circulation in the region in the future which could alter the residence time of the water masses and thus their effective age. The dynamical model (GSBM) results indicate an increase in primary production, mainly by diatoms, and a greater mesozooplankton mortality in the summer of 2050 due to the higher water temperature, which both lead to an increase in local respiration of organic matter at depth in the region. Increase in local subsurface respiration would cause a further decline in Ω_a. This kind of feedback is not represented by the empirical biogeochemistry but could be resolved by monitoring the biomass and rates of net primary and secondary production, along with their interactions with changing hydrographic conditions in the region. Targeted research to constrain these uncertainties would increase confidence in the biogeochemical model response and reduce uncertainties in the projections.

The ensemble spread of the downscaled simulations represents a measure of uncertainty in the system for the projected future conditions. The ensemble spread was small for the GSBM projections as the projected conditions largely agree in both the absolute values for future Ω_a and the magnitude of the change from present. However, the absolute values for the surface and bottom projections did vary based on the degree of warming, mostly per ensemble member (Table 2). ROMS ensemble spread is similar for the surface and subsurface, with the HadGEM-forced simulation which has the warmest SST projection producing the smallest change from the control. The HadGEM forced simulation also projects the highest Ω_a in the future (Table 3). The decline in Ω_a from control conditions due to the addition of anthropogenic CO₂ for 2050 under RCP 8.5 is –0.22 to –0.29 across the ROMS ensemble, and the majority of that is compensated by the future warming (+0.16 to +0.25). The most compensation is provided by the warmest ensemble member, HadGEM. Salinity changes in the future act to further reduce Ω_a, but the magnitude of this change is small on an annually averaged basis (–0.01). As a result, reducing the uncertainty surrounding the future warming and rates of warming in the region is essential to reducing the uncertainty on the OA projections.

All the downscaled simulations agree that the region will experience more suboptimal conditions for marine organisms for the majority of the year by 2050 in both the surface and bottom conditions. The degree to which this change occurs is modulated by the degree of projected warming and freshening in the physical simulations. All three simulations of the ROMS model projected large increases in SSTs and changes in bottom water temperature over most of the domain (Alexander et al., 2020; Brickman et al., 2021). The enhanced warming may result from a number of processes discussed in Alexander et al. (2020) including very strong warming of the atmosphere over eastern Canada that is transported over the Atlantic by westerly winds, heating the ocean via changes in the surface fluxes. In addition, strong warming occurs in portions of the North Atlantic, especially to the southeast of Newfoundland, which are subsequently advected into the GOM including through the NEC. In all simulations, the Gulf Stream is projected to be weaker (Alexander et al., 2020). Reducing the uncertainty surrounding rates of warming in the future would also reduce the uncertainty in the OA projected conditions.

Over the past 15 years, waters in the GOM have taken up CO₂ at a rate significantly slower than that observed in the open oceans due to a combination of the extreme warming experienced in the region and an increased presence of well-buffered Gulf Stream water (Salisbury and Jonsson, 2018). The reduced uptake of CO₂ by the region has also altered local acidification rates, which differ from the global rates in the GOM. In Salisbury and Jonsson (2018), the rates, averaged for the entire GOM from 1981 to 2014, showed a decline in Ω_a at the surface of –0.049 (±0.009) per decade, which is slower than global declines of Ω_a and falls outside the range reported over many sites around the world (–0.0066 to –0.0115 per year, 1982–2012; Bates et al., 2014). When the time period is limited to 2004–2015 when observations were in place at Buoy D and warming was greatest in the region, the Ω_a increase at the surface was 0.015 per year, which is consistent with the WB observational trends reported here for the subsurface (0.02–0.04 per year; Figure 3). These trends correspond to a time period when the region became increasingly dependent on vulnerable shellfish fisheries (Gledhill et al., 2015; Figure S1). In the projections presented here, the GSBM model simulates a constant rate of Ω_a decline of –0.008 per year over the upper 150 m between 1980 and 2060. In a companion paper about the future warming for the region using the same simulations discussed here (Brickman et al., 2021), warming continues to be projected and continues to ameliorate the impacts of OA in the region. However, under RCP 8.5, the impacts of OA still emerge in our projections by 2050 over the entire region with implications expected for the vulnerable ecosystems that reside in the GOM.

Developing knowledge on life stage and population-specific impacts of OA on marine organisms is critical to projecting vulnerabilities and developing adaptation plans for marine biodiversity and for sustaining ecosystems and the socioeconomically dependent coastal communities that rely on these resources. Multiple aspects of the carbonate system have been shown to have impacts on marine biota. Reductions in CO₃^2– availability and Ω_a below critical thresholds have been shown to represent a stressor to a variety of shelled marine invertebrates (e.g., Clements et al., 2017, 2018; Rheuban et al., 2018; Young and Gobler, 2018; Pousse et al., 2020) including crabs, lobsters, oysters, scallops, and mussels in the GOM (Tables S1 and S2). The pH can impact physiological processes that rely on the transfer of [H+], including olfactory senses in fish (e.g., Dixson et al., 2015) and crab megalopae (Miller et al., 2016). Finally, CO₂ concentrations in excess can cause hypercapnia (excess CO₂ in blood) in organisms (Mintenbeck et al., 2012; Macek et al., 2016; Pimentel et al., 2016; Feely et al., 2018; Davenport, 2019). Each of these carbonate system variables can cause stress for marine organisms and may be compounded by concurrent stressor combinations (Kroeker et al., 2017) of temperature, salinity, and dietary availability.

Considerable research has provided insight into how organisms from the GOM respond to OA (Tables S1 and S2), but there are still critical gaps for key environmental components, such as carbonate chemistry data at the sediment–water interface where commercial species like lobster, sea scallops, and surfclams live. Current bottom water measurements, which are meters above where the organisms reside, suggest that Ω_a is below critical levels for marine organisms living there. Although water quality does not likely improve closer to the seafloor, measurements at the sediment/water interface would provide data needed to understand the environment where benthic organisms spend the majority of their lives and reduce uncertainty in our projections. In this study, the seasonal variability and spatial differences of Ω_a suggest that the northern reaches (i.e., Scotian Shelf) and near-shore GOM will cross critical biological thresholds 10–20 years before the southern GOM (i.e., Georges Bank) due to the warming in the region (Figure 9). Experiments on spatial and seasonal impacts on marine organisms at multiple life stages are limited and require the synergy of chemical and physical oceanography with biological responses. Nearshore regions are common nursery grounds for both vertebrate and invertebrate marine species (Stevenson et al., 2014; Goode et al., 2019) but are poorly represented in our models. Subsequent projections for such critical habitats and inclusions of new studies on key regional economic species, including the understudied Jonah crab and goosefish, are needed to fully assess vulnerabilities and promote adaptive management plans.

Dependency of coastal communities on shellfish fisheries resources has increased over the past decade in the GOM, which necessitates understanding the dynamics governing the current and future trends for OA in the region. Many of these commercial resources are considered vulnerable to OA. An increase in biological sensitivity and impact studies, summarized in Tables S1 and S2 in this work, has considerably improved the understanding of potential responses of marine species to OA in combination with other climate stressors like temperature in the GOM, particularly for bivalves. Such efforts still only represent a small fraction of knowledge needed on commercially important species and scarcely address complex interactions through food webs and integrated risks of life stages through different habitats. Although surface OA trends historically have been masked by recent warming and increases in inputs from the Gulf Stream, Ω_a below suboptimal thresholds for marine organisms (Ω_a < 1.5) have already been observed in areas of high freshwater influence and subsurface zones.

The combination of regional projections detailed here indicates that Ω_a is projected to decline everywhere in the GOM by 2050 under RCP 8.5 with most pronounced impacts in locations near the coast, subsurface, and those associated with freshening. By 2050, the entire GOM will experience suboptimal conditions of aragonite saturation state (Ω_a < 1.5) for most of the year under RCP 8.5. Models with the most pronounced warming experienced the weakest decline in Ω_a as a consequence of the compensatory effects of largely warming and lesser salinity changes. Achieving reduced uncertainty of the regional warming trends in the future is a high priority to reduce uncertainty in OA projections for the GOM as well. Despite continued masking of the OA signal into the future, projections still indicate increased duration and intensity of suboptimal conditions for marine organisms for Ω_a throughout the GOM under RCP 8.5. These projections are more severe than the global projections for the same domain, emissions scenario, and time horizon but are providing downscaled regional projections. Our projections provide important spatial information for decision making in fisheries and future planning; for example, the Georges Bank is a region that remains above the threshold for the longest period. The Scotian Shelf is one of the regions that experiences the earliest sensitivity threshold, both the surface and the bottom. Future work needs to focus on a better understanding of the influence of variability in the Labrador Current, as well as potential biogeochemical feedbacks (i.e., freshwater alterations and sediment–water interactions) and potential adaptation strategies for land–ocean management for the nearshore regions.

Data in Figure 1 from the University of New Hampshire (UNH) Coastal Marine Laboratory are archived with the Northeast Regional Association of Coastal Ocean Observing Systems (http://www.neracoos.org/erddap/index.html). UNH/Pacific Marine Environmental Laboratory CO2 buoy data shown in Figure 1 are archived at the National Center for Environmental Information (https://www.nodc.noaa.gov/ocads/oceans/Coastal/unh_ts.html). UNH total alkalinity and dissolved inorganic carbon data used in the construction of Figure 4 are archived by the UNH Ocean Process Analysis Laboratory (https://eos.unh.edu/ocean-process-analysis-laboratory/data). WHOI data in Figure 4 are archived by the Biological & Chemical Oceanography Data Management Office (https://www.bco-dmo.org/project/472577). GOMECC-1 data used in Figure 4 are archived by NOAA Atlantic Oceanographic and Meteorlogical Laboratory (https://www.aoml.noaa.gov/ocd/gcc/GOMECC1/data.php); those from GOMECC-2 are archived at National Centers for Environmental Information (NCEI; https://www.nodc.noaa.gov/oads/data/0117971.xml, NCEI Accession Number 0117971); those from ECOA-1 are archived at NCEI (https://www.nodc.noaa.gov/oads/data/0159428.xml, NCEI Accession number 0159428). Decadal commercial landing statistics for New England fisheries species as reported from the National Oceanic Atmospheric Administration Fisheries Service (2020, www.st.nmfs.noaa.gov/). Output from the three Regional Ocean Modeling System physical simulations used for the projections as well as the control run is provided on a website (https://www.psl.noaa.gov/ipcc/roms/) and described further in Alexander et al. (2020).

The supplemental material provided as a Microsoft Word document (DOCX) contains additional model evaluation, model metrics, and control/base maps used to generate anomalies in the text. Finally, two supplemental tables are also included which detail the extensive review of the literature on species sensitivities in the region to OA as an update to the table in Gledhill et al. (2015).

Figure S1. Increasing importance of shellfish to New England fisheries in terms of total landed value.

Figure S2. Trends at WB7 for temperature and salinity.

Figure S3. Model evaluation at Coastal Marine Lab and Buoy D for Ω_a.

Figure S4. Model acidification annual metrics at Buoy D.

Figure S5. Total alkalinity (µmol kg^–1) seasonal CTRL for the Regional Ocean Modeling System model with empirical biogeochemistry.

Figure S6. Dissolved inorganic carbon (µmol kg^–1) seasonal CTRL for the Regional Ocean Modeling System model with empirical biogeochemistry.

Table S1. Responses of six commercially important species of the Gulf of Maine to increased ocean acidification conditions.

Table S2. Studies on organismal response of species from the Gulf of Maine to ocean acidification conditions.

Table S3. Gulf of Maine regional ocean acidification observing capacity description of observations plotted in Figure 2.

We would like to acknowledge Shani Rousseau for the extraction of the GSBM model data and for the generation of Figures 6–8 and David Brickman for fruitful discussions. We would also like to thank the constructive contributions of guest editor Stephenson and our two reviewers: Helen Gurney-Smith and one anonymous reviewer, for their help in improving this article. The views expressed in this article are those of the authors and do not necessarily represent the views or policies of the U.S. Environmental Protection Agency. Any mention of trade names, products, or services does not imply an endorsement by the U.S. government or the U.S. Environmental Protection Agency (EPA). The EPA does not endorse any commercial products, services, or enterprises.

SAS and KM acknowledge support from the NOAA Ocean Acidification Program grant NA19OAR0170351, the National Center for Atmospheric Research, which is a major facility sponsored by the National Science Foundation (NSF) under Cooperative Agreement No. 1852977 which supports the Early Career Innovators program.

CWH, JS, and DV support comes from NOAA grants NA15NOS0120155 and NA16NOS0120023. DV support also comes from NASA grant NNX17AK08G. JS also received support from NSF grant CO1658377.

KA and DL acknowledge support from the DFO’s ACCASP-Ocean Chemistry.

MAA acknowledges support from grants by the NOAA Coastal and Climate Applications (COCA) program NA-15OAR4310133 and NA-17OAR4310268. Support for this work was also provided by the NOAA Integrated Ecosystem Assessment Program.

ZAW acknowledges funding support provided by the Coastal Ocean Institute at Woods Hole Oceanographic Institution, National Science Foundation (OCE-1316040), the Pickman Foundation, and the Tom Haas Fund at the New Hampshire Charitable Foundation.

IM acknowledges the support from the U.S. Environmental Protection Agency Region 1’s Water Division in supporting her role on this work.

CB acknowledges funding support provided by the MIT Sea Grant award NA18OAR4170105.

The authors have declared that no competing interests exist.

Contributed to conception and design: SAS, JS, DKG, RM.

Contributed to acquisition of data: SAS, JS, DKG, CB, SM, KM, CWH, MA, DL, ZAW, JS, IM, DV.

Contributed to analysis and interpretation of data: SAS, JS, DKG, CB, SM, KM, CWH, MA, DL, ZAW, JS, IM.

Drafted and/or revised the article: SAS, JS, DKG, CB, SM, KM, CWH, MA, DL, ZAW, JS, IM, KA, CL, CM, MW, AP, DV, RM.

Approved the submitted version for publication: SAS, JS, DKG, CB, SM, KM, CWH, MA, DL, ZAW, JS, IM, KA, CL, CM, MW, AP, DV, RM.

How to cite this article: Siedlecki, SA, Salisbury, J, Gledhill, DK, Bastidas, C, Meseck, S, McGarry, K, Hunt, CW, Alexander, M, Lavoie, D, Wang, ZA, Scott, J, Brady, DC, Mlsna, I, Azetsu-Scott, K, Liberti, CM, Melrose, DC, White, MM, Pershing, A, Vandemark, D, Townsend, DW, Chen, C, Mook, W, Morrison, R. 2021. Projecting ocean acidification impacts for the Gulf of Maine to 2050: New tools and expectations. Elementa: Science of the Anthropocene 9(1). DOI: https://doi.org/10.1525/elementa.2020.00062

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

Guest Editor: Robert Stephenson, Fisheries and Oceans Canada Andrews Biological Station, Canada

Knowledge Domain: Ocean Science

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