Changing Hydrographic, Biogeochemical, and Acidification Properties in the Gulf of Maine as Measured by the Gulf of Maine North Atlantic Time Series, GNATS, Between 1998 and 2018.

The Gulf of Maine North Atlantic Time Series (GNATS) has been run since 1998, across the Gulf of Maine (GoM), between Maine and Nova Scotia. GNATS goals are to provide ocean color satellite validation and to examine change in this coastal ecosystem. We have sampled hydrographical, biological, chemical, biogeochemical, and bio-optical variables. After 2008, warm water intrusions (likely North Atlantic Slope Water [NASW]) were observed in the eastern GoM at 50-180 m depths. Shallow waters (<50 m) significantly warmed in winter, summer, and fall but cooled during spring. Surface salinity and density of the GoM also significantly increased over the 20 years. Phytoplankton standing stock and primary production showed highly-significant decreases during the period. Concentrations of phosphate increased, silicate decreased, residual nitrate [N*; nitrate-silicate] increased, and the ratio of dissolved inorganic nitrogen:phosphate decreased, suggesting increasing nitrogen limitation. Dissolved organic carbon (DOC) and its optical indices generally increased over two decades, suggesting changes to the DOC cycle. Surface seawater carbonate chemistry showed winter periods where the aragonite saturation (Ωar) dropped below 1.6 gulf-wide due to upward winter mixing of cool, corrosive water. However, associated with increased average GoM temperatures, Ωar has significantly increased. These results reinforce the hypothesis that the observed decrease in surface GoM primary production resulted from a switch from Labrador Sea Water to NASW entering the GoM. A multifactor analysis shows that decreasing GoM primary production is most significantly correlated to decreases in chlorophyll and particulate organic carbon plus increases in N* and temperature.

Introduction

The Gulf of Maine (GoM) is a semi‐enclosed shelf sea off of NE North America and has an average depth of about 150 m with greatest depth of 379 m (Brooks, 1992). It is bordered by coastline to the west, north, and east plus, bordered to the south and east by two shoal areas, Georges Bank and Browns Bank, respectively (Figure 1; note, a list of the acronyms and variables used in this paper can be found in Table 1). The GoM communicates with the NW North Atlantic via two deep channels, the Northeast Channel (sill depth 230 m) and the Southwest Channel (sill depth 70 m) (Brooks, 1992). Deep slope waters periodically flow through the Northeast Channel (Mountain, 2004; Ramp et al., 1985) into the GoM. This represents a deep source of inorganic nitrogen, which is subsequently mixed upwards due to intense tidal mixing and deep convective winter overturn (Hopkins & Garfield, 1979; Townsend, 1991), fueling the next season's primary production. The circulation around the Gulf of Maine is generally counter clockwise with Scotian Shelf (SS) water flowing around the eastern side of the GoM after entering along the southern flank of Nova Scotia. On the western side of the Gulf is the Eastern Maine Coastal Current (EMCC), emanating from the tidally‐resonating, well‐mixed, Bay of Fundy, flowing southwestwards along the Maine coast (Brooks & Townsend, 1989; Townsend et al., 1987). The EMCC plume turns offshore just seaward of where the Penobscot watershed empties into the GoM and here is referred to as the extension of the EMCC (ExtEMCC) (Pettigrew et al., 1998) (Figure 1). This plume is thought to provide 44% of the inorganic nitrogen to support new production in the entire GoM (Townsend et al., 1987). The Western Maine Coastal Current (WMCC) originates south of Penobscot Bay and continues southwestwards toward Massachusetts Bay (Brooks, 1985). There are three deep basins in the Gulf of Maine, all isolated below 200 m depth: Wilkinson Basin (∼275 m depth), Jordan Basin (∼275 m depth), and Georges Basin (∼370 m depth) (Brooks, 1992).

Figure 1

Region of study showing the Gulf of Maine with bathymetry. GNATS transect across the Gulf of Maine appears as a solid line crossing between Portland, ME and Yarmouth, NS. Water masses that are traversed along the transect are also indicated: Western Maine Coastal Current (WMCC), Eastern Maine Coastal Current (EMCC), Extension of the Eastern Maine Coastal Current (Ext EMCC), Jordan Basin (JB), and Scotian Shelf Water (SS). NERACOOS buoy M is designated with “⊗”. Notable bathymetric features (banks and basins) are also designated. Approximate mean residual surface circulation is shown with white arrows. Figure modified from Balch et al. (2004).

Table 1

Definitions for Acronyms and Variables Used in This Work

Term	Definition
a Φ440	Phytoplankton absorption at 440 nm (m−1)
a g412	Dissolved absorption at 412 nm (m−1)
AIRICA	Automated Infra‐Red Inorganic Carbon Analyzer
AMOC	Atlantic Meridional Overturning Circulation
a pg440	Particulate and dissolved absorption from CDOM plus phytoplankton at 440 nm (m−1)
b	Total scattering (m−1)
b bp	Particulate backscattering (m−1)
b b'	Acid‐labile backscattering from PIC (m−1)
BSi	Biogenic silica (μM)
CDOM	Chromophoric dissolved organic matter
DIC	Dissolved inorganic carbon (μequiv Kg−1)
DIN	Dissolved inorganic nitrogen ([nitrate + nitrite]; μM)
DOC	Dissolved organic carbon (mol m−3)
EMCC	Eastern Maine Coastal Current
ExtEMCC	Extension of the Eastern Maine Coastal Current
FDOM	Fluorescent dissolved organic matter (quinine sulfate equivalent units)
GNATS	Gulf of Maine North Atlantic Time Series
GoM	Gulf of Maine
GoSL	Gulf of Saint Lawrence
JB	Jordan Basin
LSW	Labrador Sea Water
MARMAP	MArine Resources Monitoring, Assessment, & Prediction
μPIC	PIC‐specific rate of calcite precipitation (d−1)
μPOC	POC‐specific growth rate (d−1)
MVP	Moving Vessel Profiler
NASW	North Atlantic Slope Water
NERACOOS	Northeastern Regional Association of Coastal and Ocean Observing Systems
PIC	Particulate inorganic carbon (mmol C m−3)
P m B	Maximum chlorophyll‐normalized primary production (g C [g chl−1] d−1)
POC	Particulate organic carbon (mmol C m−3)
PON	Particulate organic nitrogen concentration (mmol N m−3)
PSDF	Particle size distribution function
SS	Scotian Shelf
SUVA	Specific UV absorption
T TSG	Temperature as derived from shipboard underway thermosalinograph (°C)
T XBT	Temperature as derived from XBT (°C)
VINDTA	Versatile INstrument for Detection of Titration Alkalinity
Ωar	Aragonite saturation state
Ωca	Calcite saturation state
WMCC	Western Maine Coastal Current
XBT	Expendable bathythermograph

Water Mass Influences From Outside the GoM

The proximity of cold fresh Labrador Slope Water (LSW) versus North Atlantic Slope Water (NASW, aka Warm Slope Water (Townsend et al., 2015)) to the Northeast Channel affects the relative concentrations of nitrate and silicate in the water that enters the GoM, thus influencing the blooming of siliceous diatoms versus non‐silicate‐requiring phytoplankton (e.g., dinoflagellates, coccolithophores and other nanoplankton). The southward extension of LSW occurs more during low values of the North Atlantic Oscillation (Petrie & Drinkwater, 1993).

Recently, the reduction of the westward extent of LSW has been shown to be associated with the warm, saline Gulf Stream water migrating north to a region known as the Tail of the Grand Banks, a crossroads of subpolar and subtropical currents (Gonçalves Neto et al., 2021). A climatic‐dipole‐like behavior was observed, associated with an anomalous increase in sea‐surface height (observed with a 25‐year altimetry record), which was associated with an anomalous increase in sea‐surface height. This, in turn, shunted cold, fresher LSW away from the NE continental shelf in 2008, causing warm, saltier NASW to be positioned at the mouth of the Northeast Channel, potentially causing significant warming and salinification within the GoM. This observation provides the potential for a teleconnection between the Atlantic Meridional Overturning Circulation (AMOC) and the GoM hydrographic properties.

Primary Production in the GoM and Its Variability Through Time

Henry Bigelow was the first oceanographer to qualitatively describe the patterns of phytoplankton biomass (“flowering”) in the GoM (Bigelow, 1914, 1924), as well as seasonal patterns of blooming (phenology) for dominant phytoplankton species. By measuring the diurnal variability of oxygen concentrations, Riley (1941) was the first to quantitatively estimate gross and net integrated production over Georges Bank just south of the GoM. He described average rates of net productivity throughout the year, in carbon equivalents, at 400 (±370) mg C m−2 d−1 with a range of −50 mg C m−2 d−1 in January—the negative value signifies net heterotrophy—to 950 mg C m−2 d−1 in April (see his table 6). Subsequent estimates of GoM primary production using 14C are fairly rare, but have shown average integrated rates of 180–250 mg C m−2 day−1 in the GoM (Graziano et al., 2000; O'Reilly & Busch, 1984; O'Reilly et al., 1987). Balch et al. (2008) estimated depth‐integrated primary productivity of the GoM based on a combination of surface 14C productivity and satellite measurements, plus XBT estimates of the mixed layer depth as measured by the Gulf of Maine North Atlantic Time Series (GNATS), a coastal time series that has operated between Portland, Maine, and Yarmouth, Nova Scotia since 1998 (Figure 1). Applying these data to the Howard‐Yoder productivity model (Howard & Yoder, 1997) the GoM integrated productivity was seen to be on par with earlier measurements, 458 (±108) mg C m−2 d−1, which translated to an annual, regional rate of ∼38.12 Tg C y−1. This regional productivity estimate was within 2% of an independent estimate of 37.90 Tg C y−1, based on a regression analysis of other historical 14C data from the MArine Resources Monitoring, Assessment, & Prediction (MARMAP) program (1977–1982). This latter analysis also incorporated vertical diffusion, horizontal advection of nitrate, and potential new production. Historical results have reinforced the importance of nitrogen being limiting to primary production of the Gulf (Townsend & Ellis, 2010; Townsend et al., 2015).

GoM time series studies up to 2010 showed some significant decreases in net primary production in the early 2000s (Balch et al., 2012; Conor McManus et al., 2014). There are two hypotheses for the cause of this decrease in primary production, one related to nutrient limitation and one related to light limitation. The first hypothesis is that increased stratification of surface waters due to climate warming, lower salinities associated with increases in precipitation and river discharge, cause the timing of the spring bloom to change. Indeed, lower salinities over the Nova Scotia Shelf region appears to cause earlier spring phytoplankton blooms in waters upstream of the GoM, which inhibits vertical nutrient fluxes from deeper waters (Ji et al., 2007, 2008). Another hypothesis is that the significant influx of river‐born chromophoric dissolved organic matter (CDOM) following major precipitation events reduces primary production since CDOM preferentially absorbs blue light used in photosynthesis (Balch et al., 2012). Absolute discrimination between these hypotheses is not trivial since there is a significant Gulf‐wide inverse correlation between salinity and CDOM (and both properties are significantly correlated with river discharge (Balch et al., 2012)).

In this paper, we summarize GNATS observations from 1998 through 2018, the first 20 years of the time series. This extended compilation runs through the period of significant temperature increases beginning in 2004, making the GoM one of the more rapidly‐warming water bodies within the global ocean (Pershing et al., 2015, 2018; Record et al., 2019). The GNATS data set now includes some of the warmest years in the GoM on record, including the heatwaves of 2012, 2016 (Pershing et al., 2018), and 2018 (Poppick, 2018). It also includes the major elevated precipitation years of 2006–2010 which marked the beginning of a significant decrease in primary production (Balch et al., 2012). We examine this 20 years data set for hydrographic variability and also chemical (nutrients, ocean acidification), biological (chlorophyll concentration, general algal classes, carbon fixation), biogeochemical, and optical changes, to better understand the broader consequences of change in the GoM.

Methods

GNATS Ship Sampling

The data reported here were collected between September 1998 and December 2018 as part of the GNATS from various ships of opportunity, research vessels, and autonomous gliders. The methods have been described previously (Balch, 2004; Balch et al., 2004, 2008, 2012, 2016b). The reader is referred to these publications for more details on the methods. There are two limitations in the data set that should be noted. The GNATS program exclusively ran from the late spring to early fall between 1998 and 2005; winter cruises were added beginning in 2006 with typically one winter cruise per year. Also, there were no cruises in 2011 due to a funding hiatus.

Continuous Surface Underway Sampling

Briefly, most of the sampling along the GNATS transect involved surface underway measurements. Continuous underway sampling was done for above‐water radiometry from the bows of the various vessels (Satlantic microSAS and SeaBird hyperSAS instruments beginning 2016) consisting of downwelling irradiance, upwelling radiance, and sky radiance measured at 40° from nadir or zenith, respectively, measured at 90°–120° from the solar azimuth. Continuous underway measurements were also made on surface seawater collected from the ships' sea chests between 3 and 5 m depth and run through a shipboard: (a) debubbler (to remove optically‐scattering bubbles from the flow stream), (b) SeaBird temperature/conductivity sensor (to estimate temperature and salinity), (c) WETLabs ac‐9 sensor for spectral measurements of total particulate and dissolved attenuation and absorption (at 412, 440, 480, 510, 550, 610, 630, 686, and 710 nm); water was serially switched between raw streams and filtered through a 1 μm poresize filter followed by a 0.2 μm‐poresize filter in order to estimate the spectral absorption and attenuation of only colored dissolved organic matter (note, the ac‐9 was calibrated with reverse osmosis, polished water to a resistance of 18.2 MΩ), (d) a Wyatt Dawn or EOS volume‐scattering sensor which measured the volume scattering function in the forward and backward directions at 16 angles (both at ambient pH and reduced pH to dissolve calcium carbonate coccoliths), (e) a WETLabs flow‐through fluorometer to measure the fluorescence of chlorophyll a, (f) a WETLabs flow‐through fluorometer to measure fluorescent dissolved organic matter, FDOM, and (g) a HOBILabs Hydroscat II sensor for measuring spectral particulate backscattering (bbp; at 440 and 683 nm).

Discrete Water Sampling

Comprehensive discrete sampling was performed about every 30 km, or during satellite overpasses. Each discrete sampling consisted of first a vertical profile of temperature, made using a Sippican expendable bathythermograph. These profiles typically extended to the sea floor (using T‐4 and T‐10 probes, depending on water depth and ship speed). The top meter was excluded from each cast to account for the probe equilibrating to temperature (see Section 3). From 2001 to 2004 vertical profiles of temperature and salinity were made down to 100 m depth with a Moving Vessel Profiler 200 (a tethered, free‐fall device (AML Oceanographic), equipped with an AML CTD probe and WETLabs chlorophyll fluorometer which sampled these variables in 1 m increments).

At the same time as the above temperature profile, discrete water samples were taken from the seawater flow coming from the ship's sea chest. They were processed for all four parts of the carbon cycle: (a) particulate organic carbon (POC), (b) dissolved organic carbon (DOC), (c) dissolved inorganic carbon (DIC; begun in 2012), and (d) particulate inorganic carbon (PIC). Concentrations of PIC and POC were measured using techniques described in Balch et al. (2019) except volumes used for POC and PIC samples were 1 and 0.2 L respectively. Note, particulate organic nitrogen concentration (PON) was also derived from the same samples. DOC concentrations were measured as described in Balch et al. (2016b). For DIC chemistry, begun in 2008, 200 mL discrete water samples were taken from the flowing seawater line and preserved with 100 μL mercuric chloride until processing ashore for total DIC concentration and total alkalinity (TA). Analytical methods followed standardized protocols (Bates et al., 1996a, 2001; Dickson et al., 2007; Knap et al., 1993). TA was determined using a VINDTA (Versatile INstrument for Detection of Titration Alkalinity) analytical system and AIRICA (Automated Infra‐Red Inorganic Carbon Analyzer). A 25 mL or 10 mL sample volume was used for analysis of DIC samples on the VINDTA and AIRICA, respectively. Samples were also analyzed ashore using a highly precise (0.02%; 0.4 μmol kg−1) VINDTA system (Bates, 2007; Bates et al., 1996a; Bates & Peters, 2007). Both DIC and TA analyses had a precision and accuracy of ∼1 μmol kg−1 (precision estimates were determined from between‐bottle and within‐bottle replicates, and accuracy assessed using calibrated reference materials for shipboard and lab work from A. Dickson, Scripps Inst. of Oceanography, La Jolla, CA). Carbonic acid dissociation constants appropriate for temperate seas, and temperature, and salinity data were used to compute other inorganic carbon parameters (i.e., pCO2, pH, [CO3 2−] and Ω for aragonite/calcite) using the software, CO2calc (Robbins et al., 2010).

Discrete measurements included extracted chlorophyll (JGOFS, 1996), micro‐ and nano‐plankton abundance/classification using a FlowCam (Poulton & Martin, 2010), coccolithophore and detached coccolith concentration using polarized microscopy (Balch & Utgoff, 2009), macronutrient concentrations (nitrate + nitrite (hereafter referred to as DIN), phosphate, silicate) (JGOFS, 1996), primary production/calcification measurements (collected with surface bucket samples except when on high‐speed ferries which required sampling from the underway flow stream; all measurements involved half‐day incubations with 14C using the microdiffusion technique) (Balch et al., 2000, 2008), and scanning electron microscopy (Goldstein et al., 2003). Carbon‐specific intrinsic growth rates of organic matter were calculated by dividing daily photosynthetic carbon estimates by the concentration of POC. Carbon‐specific intrinsic growth rates for PIC were calculated by dividing the calcification rate by the concentration of PIC. Biogenic silica (BSi) was also sampled with all the discrete water samples (Brzezinski & Nelson, 1989).

NERACOOS Temperature/Salinity Sampling

The Northeastern Regional Association of Coastal Ocean Observing Systems (NERACOOS) maintains a series of oceanographic buoys in the Gulf of Maine, with the oldest deployments beginning in 2001 but most measurements were made beginning late 2003, 5 years after GNATS began. Buoys were equipped with sensors measuring temperature and salinity at various depths in the water column, and data are recorded and transmitted hourly via cellular telephone and the GOES satellite system. Data are maintained and served through the NERACOOS website (neracoos.org) available through multiple data interfaces. Here, we used data from Buoy M, located in Jordan Basin (longitude 67.87 W, latitude 43.49 N), close to the deepest part of the GNATS transit line. Data were averaged into monthly time series at each depth, and a monthly anomaly was calculated for both temperature and salinity by subtracting the mean for each month from the value for each month. Monthly anomalies greater than 3°C and monthly salinity anomalies greater than 1.5 psu were filtered out for data visualization.

Autonomous Glider Sampling

Two Slocum Electric Gliders (Webb/Teledyne) have been deployed along the GNATS transect in order to collect depth‐resolved information that was typically unavailable from the ships of opportunity. The first glider was deployed beginning in 2008. The measurements have been outlined previously (Balch et al., 2016b). Briefly, the first glider was a model G‐1 glider (Webb‐Teledyne) equipped with radiometry sensors for downwelling spectral irradiance (E d(λ)) and upwelling spectral radiance ((L u(λ)) at seven wavelengths: 412, 440, 480, 555, 630, 683, and 710 nm). The glider also held a WETLabs Eco‐triplet puck which measured chlorophyll a fluorescence, FDOM fluorescence, and volume scattering at 124°, used to estimate b bp(532 nm). The G‐1 glider was operated within a depth range of 0–180 m and it typically took about 3.5 weeks to transit across the GoM, round trip. The puck optical measurements were used as optical proxies and calibrated to in‐water biogeochemical variables using simultaneous collection of water samples during deployment and recovery. This allowed estimates of chlorophyll a concentration (using chlorophyll a fluorescence), DOC (based on FDOM), and POC (based on b bp532). A second glider was deployed along the GNATS line in 2012 (model G‐2 Slocum Electric Glider). This glider was also equipped similarly to the G‐1 glider above (except only four wavelengths of downwelling irradiance and upwelling radiance were measured (380, 412, 440, and 683 nm)). Additionally, dissolved oxygen was measured with an Aanderaa Oxygen Optode and nitrate was measured with a SeaBird SUNA nitrate sensor. The G‐2 glider was depth rated to reach the bottom in the GoM (maximum depth is 366 m; average depth is 139 ± 75 m). There were no deployments of gliders in 2011 due to a funding hiatus.

Data Merging

All glider and ship data were merged into a single database from which data could be evaluated in space and time. Hövmoller plots were used to explore variability in space and time where time (years) was plotted on the y axis and longitude (for the approximately east‐west GNATS transect line) on the x axis. The statistical krigging techniques used in generating the Hövmoller plots has been described (Balch et al., 2008). The JMP statistical package was used to do all statistical analyses.

Results

GNATS Data Distribution Through Time

For the 204 cruises of GNATS data made between 1998 and 2018, the sampling was biased to the months of June to October (Figure S1 in Supporting Information S1). This is because the time series was started using commercial ferries which operated across the GoM between late spring and early fall, when the weather and seas in the GoM were best. Cruises from late fall to early spring began in 2006 using small research vessels which biased those results to the latter part of the time series. The need for clear‐sky days (for ocean color satellite validation) and safe sea states for the small vessels required that the cruise frequency was reduced between late fall and early spring period.

Hydrographic Variability as Measured In Situ

The surface temperatures from XBT temperature profiles were compared with surface temperatures from the shipboard underway seawater systems for all the vessels used in the GNATS over the 20 years time series. The results demonstrated that there was no statistical difference between the data distributions (see Supporting Information S1). Given the importance of temperature in defining the physics of the GoM water masses, as well as the availability of large amounts of temperature data as a function of depth, our statistical analysis of temperature trends includes partitioning the data by longitude, season, water mass, and depth. Definition of the seasons was done according to the northern hemisphere celestial calendar (where winter begins on 21 December [winter solstice], spring begins on 21 March [spring equinox], summer begins on 21 June [summer solstice], and fall begins on 21 September [fall equinox]). For chemical, biological, bio‐optical, and biological temporal trends discussed later, the larger overall variability of these variables only warranted such partitioning in a few cases. Finally, based on previous work (Balch et al., 2008), water masses henceforth are defined to approximately encompass the following longitudes along the GNATS transect: WMCC west of 69.2°W, ExtEMCC between 69.2°W and 68.8°, Jordan Basin Water (JB) between 68.8°W and 66.7°W, and SS water east of 66.7°W longitude.

Temperature variability, as measured in situ with temperature profiles across the GoM (using XBT, MVP200 and gliders), has shown three distinct epochs (Figure 2): (a) 1998–2005 inclusive, which showed moderate temperatures within 10, 50, 100, 150, and 180 m depth strata, across the GoM and the most negative vertical temperature gradients in the upper 50 m, (b) 2006–2009 showed the coolest temperatures within a given depth strata with relatively moderate temperature gradients in the upper 50 m, and (c) 2010–2018 in which in the eastern Gulf has been consistently warming in all depth strata with near‐zero temperature gradients in the upper 50 m in 2017–2018. Regionally, the 50–150 m temperatures on the west side of the Gulf have not experienced as large a temperature increase as the eastern side. Moreover, vertical temperature gradients in the top 50 m on the eastern side of the GoM in SS waters consistently have been the least negative (approaching zero) over the entire 20 years (Figure 2f). Tables S1–S5 in Supporting Information S1 show the longitude‐ and seasonally‐partitioned 20‐years temperature trends over the time series for 10, 50, 100, 150, and 180 m, respectively, based on XBT, MVP, and Slocum Glider temperature profiles.

Figure 2

Hovmoller plots from the GNATS of in situ temperature measurements (Z axis; °C) over time (Y axis; years) and space (X axis; West longitude) at the following depths: (a) 10 m, (b) 50 m, (c) 100 m, (d) 150 m, and (e) 180 m. These plots include data collected with XBT, MVP200, and Slocum Electric gliders. (f) Temperature gradient over the top 50 m of the water column (negative values mean water is cooling with depth). Areas in white represent where no data are available either due to periods with no measurements or bathymetric constraints. There were no winter cruises funded between 1998 and 2006. Three time epochs described in text are shown along left side of time axes. Horizontal dashed lines mark the summer solstice of each year.

Stratifying Hydrographic Data by Depth, Longitude, and Season

In order to derive quantifiable relationships of temperature change, we illustrate the Gulf‐wide temperature properties measured as a function of time (years) at 10 and 150 m (Figure 3), aggregating all the cross‐gulf data. The least‐square fits showed a decrease in temperature at both 10 and 150 m depths during the second epoch (with the minimum in 2009 at 10 m (Figure 3a) and at 150 m, the minimum was 3 years earlier (2006) (Figure 3b) followed by subsequent warming in the third epoch. The description and statistics of the polynomial fit are given in Figures S2 and S3 in Supporting Information S1. Within the polynomial fit, there was still net increasing temperature across the entire length of the GNATS transect, when averaged over 20 years, with the rate of temperature increase at 10 m of 0.057 ± 0.203°C y−1 and at 150 m of 0.082 ± 0.052°C y−1 (Tables S1 and S4 in Supporting Information S1). Results from 50, 100, and 180 m also showed significant warming over 20 years; see Tables S2, S3, and S5 in Supporting Information S1, respectively.

Figure 3

Temperature (°C) versus year plot for Gulf‐wide GNATS data set, at all longitudes at (a) 10 m depth and (b) 150 m depth. A polynomial line was fit (dark blue line) through the data from each depth given the non‐linear trends in the data. Blue shaded area immediately surrounding each line represents the 95% probability limits for the mean relationship. Lighter blue shading denotes the 95% confidence intervals for individual predicted values. The statistics for the fit in panel (a) and (b) are given in the Supporting Information (Figures S2 and S3 in Supporting Information S1).

It was critical to stratify our analyses of the GNATS data, including the above temperature example, in space (by depth and longitude) and time (by season) because GNATS cruises were not performed regularly as a function of time (Figure S1 in Supporting Information S1) and because of strong differences in environmental conditions associated with the various water masses along the transect. In this and all subsequent examples, a linear fit was used to describe the 20 years temperature trends instead of the polynomial model (as in Figure 3) because the smaller sample size in each data bin made such discrimination of limited value. Moreover, this allowed the calculation of an average temperature increase or decrease (in °C y−1) in each season and longitude bin over the 20 years. At 10 m depth during the winter, summer, and fall, there was either a negligible or statistically significant increase in temperature along much of the GNATS sample line, between 69.75°W and 66.75°W (Figure 4; Table S1 in Supporting Information S1). When viewed over the entire GNATS transect, temperature at 10 m showed statistically‐significant increases during winter, summer, and fall (Table 2). For the spring cruises, however, the 10 m temperature trends were significantly negative; that is, spring temperatures decreased over the 20 years of GNATS (Table 2). The rate of decrease was −0.13 to −0.17°C y−1 in the WMCC, −0.25°C y−1 in the ExtEMCC, −0.17 to −0.19°C y−1 in JB, and −0.08°C y−1 in SS water (Figure 4; Table S1 in Supporting Information S1). Because of potential aliasing of the spring data by the addition of April and May cruises post 2006, we also calculated the average 20‐years temperature trends just for the month of June (which was sampled at a frequency similar to the maximally sampled seasons over the 20‐years data set) in the four water masses of the GoM. There was significant June cooling observed, as well (−0.13°C y−1 in the WMCC, −0.095°C y−1 in the ExtEMCC, −0.10°C y−1 in JB, and −0.08°C y−1 in SS water). Note, there was also statistically significant cooling in May and significant heating from January to April.

Figure 4

Temperature (°C) versus year plots for the GNATS time‐series at 10 m depth, stratified by season and longitude (in 0.5° Longitude bins). Least‐square fit lines are also shown in blue, with 95% confidence limits of the mean shown with darker shading around the line. Lighter shading designates 95% confidence interval for individual predicted values. The statistically‐derived mean trend in temperature (°C y−1) for each regression fit is given in each panel along with a designation of the level of significance (∅ = P > 0.05; * = P < 0.05; † = P < 0.01; ǂ = P < 0.001; § = P < 0.0001). Complete statistics for each regression are summarized in Table S1 in Supporting Information S1.

Table 2

Summary Statistics for GNATS 20‐Years Trends in Various Upper Water Column Variables Sorted by Hydrographical, Chemical/Biogeochemical, Biological, and Bio‐Optical Data Types

Data type	Y variable	Avg. annual change	Slope	Y intercept	RMSE	R 2	F test DoF denominator	F (1, DOF denom)	P value <
Hydrographical	Temperature 10 m (°C) [spring]	0.9721	−0.1815	371.7	2.16	0.133	1775	272.25	0.0001
	Surface salinity (PSU) [summer]	1.0006	0.0202	−8.51	0.68	0.034	762	27.16	0.0001
	Surface salinity (PSU) [spring]	1.0008	0.0255	−19.61	0.88	0.028	317	9.25	0.0025
	Surface salinity (PSU)	1.0011	0.0340	−36.12	0.76	0.076	1563	128.6	0.0001
	Surface density (Sigma Theta)	1.0012	0.0283	−32.78	0.89	0.040	1561	64.5	0.0001
	Surface salinity (PSU) [fall]	1.0015	0.0495	−66.96	0.52	0.298	402	170.88	0.0001
	Temperature 10 m (°C) [fall]	1.0024	0.0284	−44.76	1.73	0.008	1657	13.69	0.0002
	Temperature 10 m (°C) [summer]	1.0036	0.0471	−81.45	2.26	0.009	4453	41.71	0.0001
	Surface salinity (PSU) [winter]	1.0038	0.1229	−214.8	0.89	0.236	74	22.8	0.0001
	Temperature 10 m (°C) [winter]	1.0450	0.2252	−447.7	1.23	0.337	341	173.35	0.0001
Chemical/Biogeochemical	Log PIC/BSi (Molar ratio)	0.9382	−0.0277	55.3	0.56	0.051	729	39.4	0.0001
	Log PIC (μM)	0.9519	−0.0214	42.27	0.38	0.111	1421	176.9	0.0001
	Log DIN/DIP (molar ratio) [Summer]	0.9583	−0.0185	37.12	0.48	0.047	446	21.9	0.0001
	Log DIN/DIP (molar ratio) [Spring]	0.9600	−0.0177	36.1	0.53	0.030	215	6.7	0.0102
	Log BSi (μM) in WMCC	0.9602	−0.0177	35.41	0.47	0.029	158	4.8	0.0303
	Log DOC (mg L−1) in Ext EMCC	0.9683	−0.0140	28.28	0.17	0.138	81	12.9	0.0001
	Log DOC (mg L−1) in WMCC	0.9689	−0.0137	27.77	0.16	0.152	145	26.0	0.0001
	Log POC (μM) in Ext EMCC	0.9763	−0.0104	22.24	0.26	0.066	140	10.0	0.0001
	Log POC (μM) in WMCC	0.9798	−0.0089	19.16	0.22	0.070	252	18.9	0.0001
	Log Silicate (μM)	0.9851	−0.0065	13.25	0.48	0.006	965	5.5	0.0188
	Log POC (μM) in JB	0.9864	−0.0060	13.13	0.2	0.036	764	28.3	0.0001
	Log POC (μM)	0.9865	−0.0059	13.01	0.23	0.028	1322	38.4	0.0001
	Log POC (μM) in SS	0.9947	−0.0023	5.708	0.17	0.007	160	1.2	0.2806
	Log DOC (mg L−1) in JB	0.9989	−0.0005	1.051	0.1	0.001	419	0.2	0.626
	normalized Alkalinity (μequiv kg−1)	1.0017	3.6800	−5252	20.42	0.107	339	40.7	0.0001
	nTA/nDIC	1.0018	0.0020	−2.937	0.02	0.025	339	8.72	0.0034
	nDIC (μequiv kg−1) [summer]	1.0020	3.9110	−5922	27.34	0.058	182	11.2	0.001
	Log Nitrate plus Nitrite	1.0021	0.0009	−2.242	0.81	0.000	918	0.0	0.8483
	Log Phosphate (μM)	1.0136	0.0059	−12.44	0.3	0.012	992	11.8	0.0006
	CO3 = (μequiv kg−1)	1.0179	2.4690	−4839	26.46	0.031	341	10.9	0.001
	Ωar	1.0180	0.0415	−81.44	0.42	0.035	339	12.4	0.0005
	Log DOC (mg L−1) in SS	1.0610	0.0257	−51.53	0.23	0.196	74	18.0	0.0001
	Residual Nitrate (μM)	1.0741	0.0741	−149.7	1.6	0.062	903	60.0	0.0001
Biological	Log Abundance Coccoliths (ml−1)	0.8181	−0.0872	178.6	0.6	0.311	683	308.8	0.0001
	Log abundance 4–12 μm particles (ml−1) [fall]	0.9117	−0.0402	84.45	0.46	0.122	150	20.9	0.0001
	Log abundance diatoms (ml−1) [fall]	0.9194	−0.0365	75.8	0.78	0.039	147	6.0	0.0159
	Log Abundance Plated Coccolithophores (ml−1)	0.9217	−0.0354	73.18	0.27	0.273	639	240.5	0.0001
	Log abundance <4 μm particles (ml−1)	0.9303	−0.0314	67.04	0.35	0.091	516	51.5	0.0001
	Log Primary Prod (μg L−1 d−1)	0.9387	−0.0275	56.89	0.47	0.103	964	110.6	0.0001
	Log abundance 4–12 μm particles (ml−1) [spring]	0.9423	−0.0258	55.44	0.47	0.054	113	6.5	0.0124
	Log 2.3P.m./K (gC m−2 d−1)	0.9436	−0.0252	50.35	0.45	0.098	831	90.7	0.0001
	Log Calcification (μg L−1 d−1)	0.9450	−0.0246	48.76	0.75	0.037	738	28.0	0.0001
	Log μPOC (d−1)	0.9480	−0.0232	46	0.41	0.097	865	92.9	0.0001
	Log P m B (gC gChl−1 d−1)	0.9502	−0.0222	46.21	0.41	0.087	958	91.0	0.0001
	Log Chl (mgC m−3)	0.9847	−0.0067	13.5	0.39	0.012	1561	19.0	0.0001
	Log μPIC (d−1)	0.9957	−0.0019	2.723	0.79	0.000	728	0.1	0.7063
	Log abundance diatoms (ml−1) [winter]	1.1105	0.0455	−89.63	0.53	0.107	53	6.4	0.0147
	Log abundance dinoflagellates (armored and naked; ml−1)	1.1993	0.0789	−155.8	0.7	0.192	585	138.8	0.0001
	Log abundance 4–12 μm particles (ml−1) [winter]	1.2266	0.0887	−175.4	0.59	0.278	55	21.2	0.0001
Bio‐optical	Log a p440 (m−1)	0.9838	−0.0071	13.05	0.35	0.002	1042	16.2	0.0001
	Log b b' (m−1)	0.9867	−0.0058	8.127	0.5	0.006	1088	6.3	0.0126
	PSDF slope	1.0108	0.0420	−88.28	0.77	0.055	605	35.2	0.0001
	Log (a pg440/a Φ440)	1.0111	0.0048	−7.521	0.22	0.014	1377	20.0	0.0001
	Spectral Slope a g412 to ag440 (m−1 nm−1)	1.0185	0.0002	−0.4963	0.01	0.016	633	10.4	0.0013
	SUVA254 (m−2 (g DOC)−1	1.1384	0.1937	−37.55	0.29	0.090	715	71.1	0.0001

Note. In all Cases, the independent variable (X axis) was time in years. Within each category, data have been sorted from decreasing trends to increasing trends. Bold entries are statistically significant at a P‐value < 0.05. Seasonal statistics are given for all longitudes sampled. Statistics for specific water masses are based on the longitudinal range of each watermass (see Results Section 3.2 for specific longitude ranges).

At 50 m depth during the winter (Figure 5; Table S2 in Supporting Information S1) there were statistically‐significant increases in temperature between 70°W and 66.5°W (average 0.256°C y−1). In the spring, less than half of the eight longitude bins showed a statistically‐significant cooling. However, in the summer and fall seasons, all but one of the temperature versus time plots showed statistically significant warming (Figure 5; Table S2 in Supporting Information S1). The average temperature change for all seasons at 50 m was 0.086°± 0.021°C y−1 (±SE).

Figure 5

Temperature (°C) versus year plots for the GNATS time‐series at 50 m depth, stratified by season and longitude (in 0.5° Longitude bins). Least‐square fit lines are also shown in blue, with 95% confidence limits of the mean shown with darker shading around the line. Lighter shading designates 95% confidence interval for individual predicted values. The statistically‐derived mean trend in temperature (°C y−1) for each regression fit are given in each panel along with a designation of the level of significance (∅ = P > 0.05; * = P < 0.05; † = P < 0.01; ǂ = P < 0.001; § = P < 0.0001). Complete statistics for each regression are summarized in Table S2 in Supporting Information S1.

Evidence for statistically significant warming at 100 m depth was sporadic but the average increase in temperature across all seasons was 0.07 ± 0.018°C y−1 (Figure 6; Table S3 in Supporting Information S1). At 150 m depth, significant warming was observed at an average rate of 0.082 ± 0.012°C y−1 across all seasons and water masses (Figure 7; Table S4 in Supporting Information S1). Greatest warming was observed in the winter months in JB water. The 180 m depth data showed an overall average increase in temperature of 0.090 ± 0.013°C y−1 across all seasons and water masses (Figure 8; Table S5 in Supporting Information S1). Due to bathymetry of the GoM, the only water mass represented at this depth spanned three half‐degree longitude bins within JB water. In the winter, the average temperature change was statistically significant at 0.151 ± 0.0050°C y−1, the highest of all the seasons. Twenty‐year trends in warming at 180 m were still significant in the spring, summer, and fall (Figure 8; Table S5 in Supporting Information S1).

Figure 6

Temperature (°C) versus year plots for the GNATS time‐series at 100 m depth, stratified by season and longitude (in 0.5° Longitude bins). Least‐square fit lines are also shown in blue, with 95% confidence limits of the mean shown with darker shading around the line. Lighter shading designates 95% confidence interval for individual predicted values. The statistically‐derived mean trend in temperature (°C y−1) for each regression fit are given in each panel along with a designation of the level of significance (∅ = P > 0.05; * = P < 0.05; † = P < 0.01; ǂ = P < 0.001; § = P < 0.0001). Complete statistics for each regression are summarized in Table S3 in Supporting Information S1.

Figure 7

Temperature (°C) versus year plots for the GNATS time‐series at 150 m depth, stratified by season and longitude (in 0.5° Longitude bins). Least‐square fit lines are also shown in blue, with 95% confidence limits of the mean shown with darker shading around the line. Lighter shading designates 95% confidence interval for individual predicted values. The statistically‐derived mean trend in temperature (°C y−1) for each regression fit are given in each panel along with a designation of the level of significance (∅ = P > 0.05; * = P < 0.05; † = P < 0.01; ǂ = P < 0.001; § = P < 0.0001). Complete statistics for each regression are summarized in Table S4 in Supporting Information S1.

Figure 8

Temperature (°C) versus year plots for the GNATS time‐series at 180 m depth, stratified by season and longitude (in 0.5° Longitude bins). Least‐square fit lines are also shown in blue, with 95% confidence limits of the mean shown with darker shading around the line. Lighter shading designates 95% confidence interval for individual predicted values. The statistically‐derived mean trend in temperature (°C y−1) for each regression fit are given in each panel along with a designation of the level of significance (∅ = P > 0.05; * = P < 0.05; † = P < 0.01; ǂ = P < 0.001; § = P < 0.0001). Complete statistics for each regression are summarized in Table S5 in Supporting Information S1.

NERACOOS Observations‐Independent Temperature and Salinity Observations With Depth

The location of NERACOOS buoy M was just south of the GNATS line in JB (Figure 1) and provided the opportunity to cross‐check GNATS hydrographic trends in JB, at least for the 15 years after late 2003 when the mooring was installed. The results from different depth strata have also shown increasingly positive temperature anomalies with time after about 2010 (Figure 9a). The mean increase in temperature at each depth strata during this period ranged from 0.096 ± 0.070°C y−1 to 0.117 ± 0.099°C y−1 over the entire water column (results summarized in Table S6 in Supporting Information S1). We note positive anomalies in temperature were observed first in the 200–250 m sensors, followed about 2 years later with changes at the surface (Figure 9a). NERACOOS data for salinity at buoy M, in the middle of JB, also showed anomalously saltier water appearing deep after 2010, with those anomalies extending upwards to the surface about 2 years later (Figure 9b). The salinity record at 200 m shows a series of sharp spikes of salinity once per year or once every other year. These are likely some kind of sensor artifact, but that could not be confirmed. The seasonally‐stratified temperature data did not show evidence of spring cooling over the 15 years NERACOOS data set, unlike the 20 years GNATS data (Table S6 in Supporting Information S1). The seasonally‐stratified NERACOOS salinity data showed increasing salinity in 28 of 35 depth/time bins over 15 years. However, annual and seasonally‐averaged changes were not significant at most depths (i.e., P > 0.05; Table S7 in Supporting Information S1) suggesting that these salinity increases were near the limits of detection of the conductivity sensors.

Figure 9

Monthly mean (a) temperature anomalies (°C) and (b) salinity anomalies (PSU) observed at NERACOOS buoy M (latitude 43.49 N longitude 67.87 W) from late 2003 to 2018. This station is just south of the GNATS line in the middle of Jordan Basin. Red regions represent positive anomalies and blue regions represent negative anomalies. Data are shown from sensors at depths of 1–250 m depth. Statistics on seasonal and annual trends to these data are given in Tables S6 and S7 in Supporting Information S1.

Gulf‐Wide Changes in Surface Temperature and Salinity as Measured With the Continuous Flow System

Temperature and salinity were also measured in the continuous underway seawater stream sampled on board the various GNATS vessels from a depth of 3–5 m (Figures 10a and 10b). The Hovmoller plots showed similar warmest periods during the first and third epochs noted above. Surface salinity typically showed reduced values on the eastern and western sides of the Gulf with highest values in Jordan Basin. Temporal changes however showed lowest salinities on the western side of the transect, especially during the second epoch. These salinity changes previously were shown to be associated with several years of extreme precipitation events (Balch et al., 2012). In the central Gulf, however, the salinities have shown increasing values over the time series with a slight drop in 2016–2018 (Figure 10b). The net effect of these temperature and salinity changes is that the surface density showed a broad increase across the GoM over the 20‐years time series (along with lowest values along the western Gulf where most of the fresh water runoff was discharged) (Figure 10c). On the eastern side of the Gulf, the coldest water was usually associated with reduced salinities such that densities in Scotian Shelf water were similar to the mid‐Gulf.

Figure 10

Hovmoller space‐time plots from the GNATS of surface measurements taken at 3–5 m for (a) temperature (°C), (b) salinity (PSS), (c) density anomaly (sigma theta; ×10−3 kg m−3), (d) DIN concentration (μM), (e) silicate concentration (μM), and (f) phosphate concentration (μM). All other features of the plots are the same, as described in Figure 2.

The Gulf‐wide changes in surface salinity have shown a statistically‐significant gradual increase during the period 1998–2018 (Figure 11; Table 2) but with a period of anomalously low salinities observed between 2005 and 2010. The least‐squares linear fit of salinity (S) versus year (Y) is: S = −36.12 + 0.03395 × Y, suggesting an increasing 20‐years salinity trend of 0.034 salinity units per year. This equation has an RMSE of 0.76 salinity units, F(1,1563) = 128.61 (P < 0.0001) but a low overall r 2 of 0.076. Greatest rates of increase of salinity were observed in the winter (0.123 salinity units per year) and the slowest rates were in the summer (Table 2).

Figure 11

Surface salinity along GNATS transect (from continuous shipboard measurements) plotted against time (years). Data are shown for all locations across the GNATS transect. Dashed box highlights results between 2005 and 2010, the period of intense precipitation when anomalously low salinities were observed. Line represents least‐squares linear fit to all the data. Statistics of the fit are also shown in the inset.

Changes in Standing Stocks of Chlorophyll, POC, PIC, DOC, and BSi

Chlorophyll a concentration (proxy for the biomass of phytoplankton) showed an overall decline over the GNATS time series of 0.0067 log units per year (Figure S4 in Supporting Information S1; Table 2), however, in the first epoch, regions of highest phytoplankton biomass were found on either side of the Gulf, in SS water (to the east) and WMCC and ExtEMCC water in the west (Figure 12a). The pattern weakened with less variability across the entire GoM during the wet period (Epoch 2). Then there was a slight gulf‐wide increase in chlorophyll during the third epoch, at least until 2017–2018 when the chlorophyll concentrations fell dramatically between the western side of the Gulf to the middle (WMCC, ExtEMCC, and JB waters).

Figure 12

Hovmoller space‐time plots from the GNATS of concentrations of surface chlorophyll, particulate, and dissolved carbon plus biogenic silica. (a) Chlorophyll a, (b) POC, (c) PIC, (d) POC:PON, (e) DOC, and (f) BSi. All panels have logarithmic color scales. All other features of the plots are the same, as described in Figure 2.

Total POC (which includes living phytoplankton and detritus) showed a somewhat similar pattern to chlorophyll a (Figure 12b) with an overall net decrease over the entire 20‐years time series (Table 2; Figure S5 in Supporting Information S1). When the data were stratified by water mass, significant decreases in POC were observed over time in the WMCC, ExtEMCC, and JB and a slight, but nonsignificant decrease in SS water (Table 2; Figure S6 in Supporting Information S1). Regionally, total PIC was highest on the eastern side of the GoM with Gulf‐wide range of about 1–25 μg L−1 during the first and second epochs. After 2014, the PIC levels in the WMCC, ExtEMCC, and JB generally fell by an order of magnitude to 0.25 μg L−1, the lowest levels observed during the entire GNATS program (Figure 12c, Figure S7 in Supporting Information S1; Table 2). The POC:PON molar ratio across the GoM has, for the most part, stayed near the Redfield Ratio of 6.6 (Redfield et al., 1963). Exceptions to this were during the second epoch, the period of heavy precipitation when the ratio exceeded 20 across the Gulf in 2008. During the third epoch, POC:PON ratios remained slightly elevated over Redfield ratios (Figure 12d).

DOC concentrations on the other hand have shown no net change over the time series but still with large contrasts in space and time. Epoch 1 showed lowest values across the GoM (note, the DOC measurements began in 2004). The wet epoch 2 showed strong inputs of DOC into the Western GoM from Maine rivers. After 2010, the highest DOC concentrations were observed in the SS water of the eastern GoM with concentrations over one order of magnitude greater than in the western GoM (Figure 12e), however, this water did not have strong CDOM absorption (Figure 13a). When the DOC values were stratified by water mass, there was a statistically significant decrease in the WMCC and ExtEMCC, no change in JB and a significant increase in SS water (Table 2).

Figure 13

Hovmoller space‐time plots from the GNATS of surface (a) dissolved absorption at 412 nm (a g412; m−1), (b) total backscattering at 532 nm (b btot532; m−1), (c) particle absorption at 440 nm (a p440; m−1), (d) acid‐labile backscattering at 532 nm (b b’532; m−1), (e) ratio of particulate and dissolved absorption at 440 nm normalized by the phytoplankton absorption at 440 nm, and (f) ratio of acid‐labile backscattering at 532 nm to the total backscattering at 532 nm (dimensionless). All other features of the plots are the same, as described in Figure 2.

The particulate BSi concentration showed no statistically‐significant 20‐years trend over the time series except for a decreasing trend in the WMCC (Table 2). BSi concentrations were almost always highest in the eastern and western water masses leaving the JB with the lowest BSi concentrations (Figure 12f). In the eastern Gulf during Epoch 1, highest BSi values were seen in SS water which then switched to western waters during Epoch 3. However, the lowest BSi values were observed in the central JB during 2017–2018 (Figure 12f).

Changes in Nutrient and Biogeochemical Indices Through Time and Space

Surface nutrient concentrations were typically elevated on the eastern side of the transect (Figure 10d). Moreover, nutrient concentrations were patchy in space and time over the time series but generally elevated during the second epoch relative to the first and third epochs (Figures 10d–10f). Surface DIN was frequently drawn to undetectable levels in the western Gulf in epochs one and three. Silicate concentrations increased in the western part of the transect during the second and third epochs (Figure 10e).

Ratios of nutrients and certain biogeochemical variables provide insights into the limiting nutrients as well as the relative composition of different particles. For example, the ratio of DIN to inorganic phosphate (aka dissolved inorganic phosphate, DIP) (Figure 14a) showed patchy distributions in space and time with mean summer values of 1–2 (range 0.03–20; Figure S8 in Supporting Information S1) as nitrate concentrations were drawn down, suggesting nitrogen limitation (Figure S8 in Supporting Information S1). Typical winter DIN:DIP ratio values averaged around the Redfield ratio of 16 with mean spring values of 4 (range 0.1–20) and mean fall values of 3 (range 0.02–20) (Figure 14a; Figure S8 in Supporting Information S1). There was a decrease in the average DIN:DIP ratio over the 20‐year time series, but only significant in the spring and summer months, suggesting increased nitrogen limitation with time (Figure S8 in Supporting Information S1; Table 2). Ratios of silicate to DIN typically showed values of 1–10 (Figure 14b) and residual nitrate ([DIN]‐silicate; N*) showed values between −2 and 0 (Figure 14c), with N* values increasing over time (Table 2; Figure S9 in Supporting Information S1) suggesting that silicate was in overabundance relative to DIN early on, which changed to more DIN dominance later in the time series. The most negative values of N* were observed during the years 2005–2010 when the low salinity waters were observed on the western side of the GNATS transect (Figures 10b and 14c). The PIC:BSi ratio provides the relative balance of the two major forms of ballast found in particles in these surface waters (Figure 14d). The average PIC:BSi ratio decreased through the time series from 0.7 to 0.3 (range 0.04–100; Figure S10 in Supporting Information S1) suggesting that in terms of ballast‐driven export, BSi would have played an increasing role on a molar basis especially in the western GoM (Table 2). The PIC:POC ratio consistently reached 0.1 in the eastern side of the Gulf, off Yarmouth, NS over the 20 years of the measurements (Figure 14e). In the western and central portions of the Gulf, however, more typical PIC:POC ratios were 0.01–0.03, with some of the lowest ratios observed since 2015 in the mid and western parts of the GoM (0.001–0.003). There is a significant decrease in the mean PIC:POC ratio over the 20‐year time series (Table 2; Figure S11 in Supporting Information S1). Significant decreases in PIC:POC were only seen in Spring through Fall seasons, and in WMCC and JB water masses (data not shown). The POC:Chl a mass ratio indicates the relative health of the phytoplankton community (where POC here includes both living and detrital organic fractions). The lower the value, the healthier the phytoplankton population (values of <10–30 can be found in algal cultures in logarithmic growth, while values of 300–1000 are more characteristic of nutrient‐stressed, slow‐growing populations found in oligotrophic areas (Geider, 1987)). The overarching pattern in our POC:Chl a data was of low values at the eastern end of the GNATS transect (SS waters), moderate values in the WMCC and ExtEMCC waters, with highest values in Jordan Basin waters (Figure 14f). There was no significant change in the mean POC:Chl a of ∼150 g POC (g Chl a)−1 over the 20‐year GNATS, however.

Figure 14

Hovmoller space‐time plots from the GNATS of derived biogeochemical‐relevant ratios/differences, taken at 3–5 m: (a) (DIN):DIP molar ratio, (b) silicate:DIN molar ratio, (c) residual nitrate DIN‐silicate, (d) PIC:BSi molar ratio, (e) PIC:POC ratio, and (f) POC:Chlorophyll mass ratio. All other features of the plots are the same, as described in Figure 2.

Decreases in Carbon Fixation and Growth Rates Over the Time Series

Highest primary productivities (1000–1700 mg C m−3 d−1) were observed in the first epoch of the GNATS program (2001–2006), on either side of the GoM, in WMCC, EMCC, and SS water. These remain the highest productivities of the study (Figure 15a). Low values of production of 16 mg C m−3 d−1 typically were found in the central JB waters during this time. Further, during 2012–2015, the productivity of even the SS water showed low values (Figure 15a). A plot of primary productivity versus time showed that the mean production value has decreased from 80 to 30 mg C m−3 d−1 between 2001 and 2018, when primary production was sampled (Figure S12 in Supporting Information S1; Table 2). Calcification rates were overall about two orders of magnitude lower than primary production rates (Figure 15b). As with primary productivity, highest calcification values were seen in 2001–2002 and such high values have not been seen since. There was a sharp drop‐off in calcification rates seen from 2008 to 2010, virtually across the entire GoM. Calcification values in the third epoch rebounded a bit from the second epoch, but there is still an overall net decrease in calcification rates from 0.4 to 0.2 mg C m−3 d−1 during the time series (Figure S13 in Supporting Information S1, Table 2). For calcification:photosynthesis ratios, with the exception of values approaching 1 for parts of 2001, ratios of 0.001–0.1 were typical, with little significant geographic or temporal variability (Figure 15c). The maximum chlorophyll‐normalized primary production rates (P m B) also showed little geographic variation across the GoM, with typical range of ∼10–400 mg C mg Chl a −1 d−1 during the first epoch of the program, dropping somewhat after 2007. In 2018, a few samples showed elevated values similar to the early GNATS program but otherwise GNATS P m B values have never recovered following the decline in 2008 (Figure S14 in Supporting Information S1; Figure 15D). There was a significant decrease in the mean P m B from ∼60 to ∼20 over the GNATS study (Table 2; Figure S14 in Supporting Information S1).

Figure 15

Hovmoller space‐time plots from the GNATS of surface (a) primary productivity (mg C m−3 d−1), (b) calcification, (c) calcification:photosynthesis ratios, (d) chlorophyll‐normalized primary production (mg C (mg Chl a)−1 d−1), (e) carbon‐specific growth rates (d−1), and (f) carbon‐specific intrinsic growth rate for PIC (d−1). All other features of the plots are the same, as described in Figure 2.

The POC‐specific growth rate (μPOC) showed Gulf‐wide rates of 0.04–2 d−1 from 2001 to 2007 after which rates lowered to 0.1–0.2 d−1 except in the SS water from 2008 to 2011 (Figure 15e). Thereafter, μPOC rates have shown subtle increases but not as high as the rates during the early part of the time series, except in the eastern side of the GoM in SS waters. Overall, μPOC values have decreased, on average, by 60.9% over the 20 years GNATS (Figure S15 in Supporting Information S1; Table 2). On the other hand, PIC‐specific rates of calcite precipitation (μPIC) ranged between 0.1 and 1.0 d−1 for much of the study with patchy (in space and time) lower values of 0.03–0.1 d−1 from 2008 to 2014 but the mean values of 0.1 d−1 showed no significant change over the 20 years time series (Table 2).

Optical Properties Across the GNATS Transect

Inherent optical properties generally demonstrated strong contrast of the second epoch on the overall GoM optical properties. For example, dissolved absorption at 412 nm (a g412; note, 412 nm is the wavelength for many space‐based reflectance measurements dominated by CDOM) was highest during this time near the Maine coast (out to 69°W) while after this epoch, values dropped to levels seen in the early parts of the study (Figure 13a). Overall, there was no significant 20 years trend in the mean a g412 values over the time series. Total backscattering at 532 nm also showed this same relative pattern as well as consistently elevated values in the SS water in the eastern GoM (Figure 13b) but with no significant 20 years trend, overall. Acid‐labile backscattering from PIC (b b’) showed consistently high values in the SS waters over this two decade study (Figure 13d) with a similar pattern to PIC (Figure 12c). Early‐on in the program, there were elevated values of bb’532, gulf‐wide, associated with coccolithophore blooms. The bb’532 values otherwise remained low to moderate (Figure 13d) with PIC accounting for 10%–30% of the total particle backscattering across the GoM (Figure 13f). There was a slightly negative, statistically significant, 20‐years trend in b b’ (Table 2). Particulate absorption at 440 nm (a p440) also showed enhanced values on the western side of the Gulf during the second epoch. However it also showed uniformly low values in the JB waters over the 20‐years time series, consistently enhanced values at 68.5°W–69°W (ExtEMCC), slightly enhanced values at 66.8°W (boundary between JB and SS waters), and increasing values eastward of 66.5°W (SS waters; Figure 13c). Averaged over the entire 20 years of GNATS, a p440 showed a significant negative trend (Figure S16 in Supporting Information S1; Table 2). The expected phytoplankton absorption at 440 nm (a Φ440) was calculated based on equation 2 from Bricaud et al. (1998), a Φ440 = 0.0378 Chl0.627. The ratio of the total measured particulate and dissolved absorption at 440 nm (a pg440) to a Φ440 is shown in Figure 13e. The dynamic range of the a pg440:a Φ440 ratio was 0.4–16, with highest values of 20 observed during the second epoch near the coast. Otherwise, values of the a pg440:a Φ440 ratio on the western and eastern sides of the GoM hovered near 4–6 for most of the 20‐years time series and values in JB waters ranged from roughly 1–2.5. There was a significant positive trend in a pg440:a Φ440 over 20 years (Table 2). Such values show the dominance of CDOM absorption over phytoplankton absorption in the GoM.

Significant Variability in Carbonate Chemistry Over Events, Seasons, and Years

Carbonate chemistry and ocean acidification (OA) parameters have only been measured in GNATS since 2012 (i.e., only seven years, inclusive, of data are available here), thus the long‐term conclusions of the results are more limited than for the other GNATS variables sampled for 20 years. Alkalinity values were highest in the JB and ExtEMCC waters and lowest in WMCC waters with moderate values in the SS waters (Figure 16a). There was strong correlation of TA with salinity (Figure S19 in Supporting Information S1) and strong seasonality observed in salinity‐normalized total alkalinity (nTA; normalized to 32 salinity units), salinity‐normalized DIC (nDIC), and their ratio (nTA:nDIC; Figure 17). The nTA and nDIC varied seasonally and out of synch, plus there was a west/east difference in their seasonal pattern. That is, nTA showed peak concentrations in the spring and minimum values in the late fall whereas nDIC showed a winter peak and a well‐defined minimum in late summer. Moreover, western and central waters (WMCC, ExtEMCC, and JB) showed lower nTA values and similar seasonal oscillations while eastern SS waters showed higher nTA and nDIC values, which decreased over the calendar year. The annual range of nTA was greater in western and central waters (45–65 μmol kg−1) than in eastern SS waters (25 μmol kg−1). The range of nDIC was greatest in JB waters (105 μmol kg−1) and less in WMCC, ExtEMCC, and SS waters (70–75 μmol kg−1). The resultant ratios of nTA:nDIC showed highest values in the summer and lowest values in fall and winter, whereas in SS waters, the ratios increased over the year. The pH showed overall seasonal variability with highest values in the summer and lowest values typically in the winter and late fall (Figure 16c). Overall pH variability of about 0.2 pH units was seen in the western water masses over all seasons whereas eastern SS water showed more muted pH variability (Figure 17b). Lowest CO3 = values of <60μeq kg−1 were observed in the winter and highest values (>140 μeq kg−1) were observed in the summer and early fall (Figure 16d). Aragonite saturation showed a seasonal summer high (associated with warm temperatures) and winter low, and, for the most part never showed undersaturation (Figures 16f and 17d). Seasonal pCO2 values showed that the western Gulf is a slight CO2 sink or close to neutral for atmospheric CO2 and the eastern Gulf is a slight CO2 source (Figure 17f) relative to the mean global atmospheric CO2, which averaged 401.5 ± 5.38 ppm for the years 2012–2018 (Dr. Pieter Tans, NOAA/GML (gml.noaa.gov/ccgg/trends/)).

Figure 16

Hovmoller space‐time plots from the GNATS of surface OA‐related variables: (a) total alkalinity (μequiv kg−1), (b) DIC (μequiv kg−1), (c) pH, (d) CO3 = concentration (μequiv kg−1), (e) Ωcalcite, and (f) Ωaragonite. All other features of the plots are the same, as described in Figure 2.

Figure 17

Surface OA variables shown as a function of season for all years sampled (where W, S, S, and F on the X axes indicate winter, spring, summer, and fall, in their chronological order, respectively). Data is stratified by water mass (WMCC = Western Maine Coastal Current; ExtEMCC = Extension of the eastern Maine Coastal Current; JB = Jordan Basin water; SS = Scotian Shelf Water). (a) Normalized total alkalinity (nTA), (b) pH, (c) normalized dissolved inorganic carbon (nDIC), (d) saturation state for aragonite (Ωar), (e) ratio of nTA:nDIC, and (f) pCO2. All data are fit with a cubic spline curve. Panels (a), (c), and (e) show horizontal lines which designate the upper and lower range of the spline fits and the change is designated with a Δ symbol. In panel (d), the value of Ω ar = 1.6 is shown with the horizontal coarse dashed line and the value of one is designated with the horizontal fine dashed line (see text for explanation). In panel (f), the horizontal dashed line designates a pCO2 of 401.5 μatm, the global average pCO2 for the time period when the OA variables were being measured, in order to define periods where the GOM was a net source or sink for CO2 (see text for details).

There were statistically significant positive trends in nTA, the ratio of nTA:nDIC, and saturation state (Table 2; Figure 18). There was no significant trend, however, for nDIC, pH, or pCO2. Concentrations of nDIC generally showed high values Gulf‐wide during winter months but this was especially so during two cruises in December 2014 and January 2015 winter (Figures 17 and 18c). Gulf‐wide, there was no significant annual trend in DIC although just in winter months there indeed was a significant increase (Table 2; Figure 18c). There was no overall trend in pH, either, whether data were stratified by water mass or season (Figure 18b). There was a statistically significant Gulf‐wide increase in CO3 = from 2012 to 2018 (Figure S17 in Supporting Information S1; Table 2). These OA parameters resulted in the calcite saturation state (Ωcalcite) of the GoM surface waters varying from 1.5 to 4.5 and with lowest values observed gulf‐wide during the winter of 2015 (Figure 16e). The aragonite saturation state (Ωaragonite) showed the same relative pattern as Ωcalcite with values of the latter ranging from 1 to 2.5, with lowest values observed in the winter and highest values in the summer months. As with Ωcalcite, the absolute lowest values of Ωaragonite (<1) were observed gulf‐wide during the winter 2014–2015 event (Figures 16f and 17d). There was a statistically significant increasing seven‐years trend in Ωaragonite over the GNATS time series, albeit with a low correlation coefficient (Table 2; Figure 18d) Ωaragonite was much more significantly correlated with temperature (measured in °C; Figure 19). The least‐squares linear fit to the data was: Ωaragonite = 0.0757 × temperature + 1.142 (RMSE = ±0.30; r 2 = 0.49; P < 0.0001).

Figure 18

OA variables for surface waters plotted as a function of year. Least square linear fits for average trends are shown in dark blue colors. Confidence limits for individual data points are shown in lighter blue shading. (a) Normalized total alkalinity (nTA), (b) pH, (c) normalized dissolved inorganic carbon (nDIC), (d) saturation state for aragonite (Ωar), (e) ratio of nTA:nDIC, and (f) pCO2. Arrows below X axes designate two consecutive cruises in December 2014 and January 2015 where waters had anomalously low nTA, nTA:nDIC ratios, pH, Ωar, and high nDIC and pCO2 were observed. Statistical trends are given in the upper left insets.

Figure 19

Surface saturation state for aragonite (Ωar) shown as a function of temperature over the GNATS program. The value of Ωar of 1.6 is shown with the horizontal blue‐dashed line. Such values are considered corrosive to aragonite shell formation (Salisbury et al., 2008). Least squares fit to the data given in the upper left of the panel.

Changes in Particle Counts Through Space and Time

The FlowCam results showed fundamental changes in the broad particle types over the 20‐year time series (Figure 20). Concentrations of the smallest particles (<4 μm diameter) were not reliable pre‐2008 (due to the use of a 5X objective on the instrument instead of the 10X objective; therefore, the data are not presented). Post‐2008, the 10X objective was used with the FlowCam and <4 μm diameter particles could be reliably enumerated. However, the image resolution still did not allow any sort of classification into general algal classes of this size class. The FlowCam results for <4 μm diameter particle abundance showed a range of concentrations from 1.6 × 103 to 4 × 104 particles mL−1 (Figure 20a). Their geographic variability across the GoM generally showed somewhat elevated concentrations in the WMCC, ExtEMCC, and SS waters, with slightly lower concentrations in the JB waters. During 2009, peak concentrations were observed, ∼15% higher than other times (Figure 20a). The overall mean abundance of this particle size class significantly decreased through the 20‐years GNATS, from 10,800 to 4,600 particles mL−1 (Figure S18 in Supporting Information S1; Table 2).

Figure 20

Hovmoller space‐time plots from the GNATS for surface waters showing the broad characterization of the particle assemblages as measured with the FlowCAM: (a) concentrations of <4 μm particles mL−1 (note, through the end of 2007 the enumeration of <4 μm particles was not reliable due to the use of a 5X objective, therefore no data are shown for this time period), (b) concentrations of 4–12 μm particles mL−1, (c) concentrations of coccolithophores mL−1, (d) concentrations of diatoms mL−1 (pennate + centric), (e) concentrations of dinoflagellates mL−1 (naked + armored), and (f) slope of the particle size distribution function. All other features of the plots are the same, as described in Figure 2.

For 4–12 μm particles (“nanoeukaryotes”), the overall variability was 30–30,000 particles mL−1 (Figure 20b). Concentrations of 4–12 μm particles increased an order of magnitude by 2009 in the summer months, and stayed high through 2017 then decreased about an order of magnitude during 2018. There was low variability in the concentrations of these particles across the entire transect (with relatively constant concentrations across the GoM) and no significant overall mean trend over the GNATS study. Seasonally, however, there was a significant increasing trend in the winter months, and decreasing trend in spring and fall (Table 2).

In the first epoch, the concentration of coccolithophores (as determined by polarized microscopy) were greatest in WMCC, ExtEMCC, and SS waters along the GNATS transect, with lowest concentrations in JB waters (Figure 20c). After 2009, coccolithophore concentrations dropped about a factor of 4–10 and remained at that level until the end of 2017 followed by a slight increase in concentration to 150 mL−1 in 2018 (Figure 20c). Overall, the concentrations of plated coccolithophores and detached coccoliths have significantly decreased in the GoM (Figures S22 and S23 in Supporting Information S1; Table 2).

Diatoms showed a range of concentrations from 10 to 104 mL−1 between 2004 and 2018. Diatom concentrations showed a pronounced increase to 104 mL−1 from 2008 to 2009, with another peak year in 2014 and 2016. Regionally, diatoms showed highest concentrations in WMCC, ExtEMCC, and SS waters and lowest values in JB waters (Figure 20d). Seasonally, diatom concentrations have increased significantly in the winter and decreased in the fall (Table 2).

The mean concentration of dinoflagellates (both naked and armored) varied over two orders of magnitude (from 200 to 4,200 mL−1) with low variability in concentrations across the different water masses of the Gulf (Figure 20e). The biggest change in dinoflagellate concentration happened in 2008 across the Gulf when average dinoflagellate concentrations rose from several hundred mL−1 to several thousand mL−1 (Figure 20e). The change in mean dinoflagellate concentration was highly significant overall (Table 2) and was observed over all seasons and water masses. The dinoflagellate concentrations remained high through 2018.

The slope of the particle‐size distribution function was typically −3 to −4 across the entire GoM, punctuated with two periods where the slope decreased to −2 in 2007–2008 and 2010 (Figure 20f). The mean particle size distribution function (PSDF) slope significantly increased overall from −4.1 to −3.6 (Table 2) and those changes were significant in ExtEMCC and JB waters but not WMCC and SS waters.

Discussion

Hydrography‐Incursions of Warm Water at Depth

Taken as a whole, the surface GoM showed a curvilinear change in temperature over the 20‐year of the GNATS (Figure 3a). The 10 m surface water in the GoM, when averaged over all seasons cooled from 12.2 to 8.9°C between 1998 and 2009 (cooling rate of −0.29°C per year), however, this estimate was potentially biased because we did not begin winter sampling until 2006. Only including the pre‐2006 data in the analysis (with no winter sampling) gave a slower cooling rate of −0.19°C per year (p < 0.0001). Only including post‐2006 GNATS data in the analysis of the data in Figure 3a, the average surface temperatures increased from 10.0°C in 2006 to 12.6°C in 2018 over the 12 years for a net warming rate of 0.22° y−1, almost identical to the SST increases observed by satellite between 2004 and 2013 of 0.23°C y−1 (Pershing et al., 2015). The 0.01°C y−1 difference between GNATS and satellite‐SST data is small given (a) the GNATS XBT temperature data were averaged over the top 10 m of the water column while the Pershing et al. (2015) data were from satellite‐derived infrared measurements, representative of the surface few millimeters of the water column (Minnett, 2003), (b) the average time period of these two calculations was not exactly the same, and (c) GNATS data are representative of the transect between Portland, ME and Yarmouth, NS whereas the Pershing et al. data were representative of the entire GoM. Regardless, averaged temperature data from GNATS are consistent with the observation that the GoM is warming at a rate similar to the three other rapidly warming ocean regions on Earth documented between 2004 and 2013 (waters west of Australia, western Pacific, and Barents Sea regions) (Pershing et al., 2015). The overall GoM temperature trend observed in the Pershing et al. data set also showed cooling between 1999 (first full year of GNATS) to 2004 (Pershing et al., 2015; see their Figure 1a), just as we observed (Figure 3). Friedland et al. (2020) also observed significant warming of surface waters in the eastern Gulf of Maine of about 0.1°C y−1 between 2004 and 2018 (based on their Figure 2). They describe a change point in the Gulf of Maine beginning in 2008 for bottom waters. For the GNATS data, the mean temperature at 150 m depth (Figure 3b), was coldest in mid‐2006, at the beginning of epoch 2. The mean fit to 10 m GNATS data reached coldest temperatures in mid‐2009. Our data, as well as the NERACOOS buoy M data suggest that the deeper temperatures started increasing 2–3 years earlier than the surface temperatures.

Seasonally‐binning the temperature data presented a more complicated picture of only surface warming of the GoM, however (Figure 4). Seasonal data showed most pronounced surface warming beginning in 2009, but still a net increase over the 20 years duration of GNATS in surface temperature of 0.22°C y−1 in winter, 0.05°C y−1 in summer, 0.03°C y−1; in the fall (statistically significant in many but not all longitude bins). This was contrasted by net cooling in the spring seasons over the 20‐years time series (−0.18°C y−1; significant in all longitudinal bins). Removing any potential aliasing in the spring by focusing only on the month of June (which was well sampled over the entire 20‐years time series) gave the same result. A recent study of the entire NE continental shelf (Friedland et al., 2020) revealed slower warming in spring months, compared to fall months. This is in contrast with the GNATS record, which showed significant cooling in spring months and warming in fall months in the surface and 50 m layers (Figures 6 and 7; Tables S1 and S2 in Supporting Information S1), albeit spring cooling at 50 m was more muted in the GNATS record than in the top 10 m. The fact that the NERACOOS temperature data showed no evidence of spring cooling while the GNATS data do, suggests that the inclusion of the pre‐2003 GNATS data were critical in driving this seasonal trend. Indeed, when we eliminated all pre‐2003 data from the GNATS data set, the 20‐year spring cooling trend for surface waters disappeared and was replaced by a significant, but small, 15‐year trend of increasing spring temperatures (0.061°C y− 1; P < 0.0001; results not shown). Bigelow (1927) (p. 546) commented on the importance of vernal warming, in which the Gulf warms from below, during the early spring as warmer, saltier, deep slope water entering through the Northeast Channel can be mixed upwards. He also discussed (p. 547) how vernal warming could be slowed or reversed in the eastern GoM associated with the shallow, colder, fresher SS current flooding into the Gulf of Maine past Cape Sable. Feng et al. (2016) also demonstrated (using altimetry and NERACOOS buoy salinity data), the importance of upstream, near‐surface, cold, fresher SS water periodically entering the GoM and influencing the EMCC. Thus, the cooling vernal temperatures in GNATS that we observed could indicate an increasing vernal influence of the SS current flooding into the GoM over this 20 year period.

While surface waters showed average warming in the GoM during winter, summer and fall, and cooling in the spring, deeper waters showed more consistent warming with decreasing variance at depth regardless of season (Figures 4, 5, 6, 7, 8; Tables S1–S5 in Supporting Information S1). The fact that spring cooling was fastest at 10 m and only moderately discernible at 50 m suggests surface currents/meteorological forcing might be driving this phenomenon. Moreover, consistent warming from 100 to 180 m, regardless of season, and salinification at depth (Figure 9; Tables S6 and S7 in Supporting Information S1) is probably best explained by an increase in the advection of warmer North Atlantic Slope Water (NASW) into the GoM through the Northeast Channel over the 20‐years time series. Notable in these observations is that the warming rates were often (but not exclusively) faster on the eastern side of the GoM within a season (closer to the source of NASW) (Friedland et al., 2020) (Tables S1–S5 in Supporting Information S1; Figures 4, 5, 6, 7, 8). NASW is typically saltier than GoM water (Gonçalves Neto et al., 2021). Combining the observations of warming surface water with increasing salinity (+0.034 salinity units y−1 for surface waters, except during extremely wet second epoch; Figure 11), would also suggest that over annual winter mixing cycles, the deep intrusions of warm, salty NASW ultimately had been mixed into the surface layer. The NERACOOS data at Buoy M also support this conclusion (Figure 9; Tables S6 and S7 in Supporting Information S1). Alternatively, increasing temperature and salinity could also have been due to mixing of NASW into SS water outside the GoM, with subsequent advection into the GoM (Townsend et al., 2015).

The balance of deep NASW at the entrance to the NE Channel over Labrador Sea Water has been thought to affect the balance of dinoflagellates over diatoms, respectively (Townsend et al., 2014, 2015). Recently this balance has been attributed to the relative position of the Gulf Stream (Seidov et al., 2021) which has shown a northward excursion after 2005, especially in the summer months. Additionally, when the Gulf Stream moves north to a region known as the Tail of the Grand Banks, it blocks the cold, fresher Labrador Sea water from its southwesterly flow along the Scotian Shelf, instead forcing it eastward, away from the North American continental shelf and thus allowing warm, saltier NASW nearer to the entrance of the NE Channel where it can flow into the Gulf of Maine. This creates a dipole‐like pattern in sea surface height, increasing ∼10.2 cm in 2008 which brought anomalously warm NASW into the Gulf of Maine 14 months later (Gonçalves Neto et al., 2021). In GNATS, this was observed between 10 and 150 m (Figures 4, 5, 6, 7).

Chemistry

Changes in Nutrients and Nutrient Indices

The time‐space variability of DIN (Figure 10d; note logarithmic scale) showed that surface DIN was consistently drawn‐down to values < 0.1 μM in the summer months of the first and third epochs, in all water masses except the SS waters. During the second, high precipitation epoch, the DIN concentrations generally were only drawn‐down to 1–3 μM in the summer months. Nitrogen limitation of phytoplankton populations of the GoM is consistent with previous reports (Townsend et al., 2014). Interestingly, SS waters rarely were depleted of DIN, DIP or Si(OH)4 suggesting a continuous upwelled source, previously postulated at the border between SS and JB waters (Balch et al., 2008, 2012).

The molar ratios of DIN:DIP (Figure 14a, Figure S8 in Supporting Information S1; Table 2) were almost always <16 (Redfield ratio), except during winter months during periods of intense mixing, again consistent with overall nitrogen limitation during non‐winter months. The small but significant decreasing ratios through time suggest increasing nitrogen limitation, at least in spring and summer (Figure S8 in Supporting Information S1). Further, the ratio of Si(OH)4:DIN was usually well‐above 1 (Figure 14b) and the N* generally stayed remarkably close to zero for the entire 20‐years time series (Figure 14c) but showed a small but significant increase over time (Table 2; Figure S9 in Supporting Information S1). This would likely result from low absolute concentrations of the two nutrients in surface waters such that their relative variability could be substantial, yet differences between them stayed close to zero. These results suggest that silicate limitation of diatoms was unlikely in GoM surface waters early‐on but potential silicate limitation has increased through time. It has previously been shown that diatoms cannot completely deplete silicate from their growth medium below concentrations of about 1 μM (Paasche, 1973a, b). Our GoM data are consistent with this (Figure 10e).

The only large perturbation of the N* in the GoM was during the wet years (Epoch 2), on the western side, when massive river discharge caused negative N* values of −4 to −6, presumably associated with massive silicate efflux from rivers associated with major river discharge events (following terrestrial weathering; Figure 14c). The overall significant decline in silicate over the 20‐years time series is consistent with the declining influence of Labrador Slope Water and increased influence of NASW being advected into the Gulf (Gonçalves Neto et al., 2021; Townsend et al., 2014).

Biogeochemistry‐Variability in Ballast Minerals

Following the wet years (second epoch), the Gulf of Maine has seen a gradual diminishing of the ratio of coccolithophore versus diatom ballast minerals (PIC:BSi ratio) from 0.7 in the early time series, to 0.3 by the third epoch, especially as PIC standing stock decreased (Figure 12c, Figures S7 and S10 in Supporting Information S1; Table 2). Notably, there have been few, if any, high‐reflectance coccolithophore blooms in the GoM since 2000 (Balch, 2004). It is also noteworthy how low the PIC:BSi ratios are in the GoM relative to other ocean regions. Virtually all of the GNATS PIC:BSi values are lower than values observed across the entire Atlantic Meridional Transect (from 50°N to 50°S at ∼20°W), illustrating how dominant BSi is as a ballast mineral in the GoM (Balch et al., 2010). This has direct ramifications to lower sinking velocities of ballasted particles in the GoM due to BSi being less dense than PIC (Balch et al., 2010; Klaas & Archer, 2002).

Ocean Acidification and Productivity‐Changes Over Annual and Event Time Scales

There were two striking aspects of the variability in ocean acidification observed in the seven years after 2012 when such data were taken during the GNATS. First, the effect of the dilution of the coastal total alkalinity signal by river outflow was particularly obvious in the Hovmoller plot (centered on 2014 and 2017; Figure 16a). This has been observed before off of Maine rivers (Salisbury et al., 2009). This also corresponded to salinity and density minima on the western side of the GoM (Figures 10b and 10c). Earlier in GNATS (2007–2010), strong river discharge was associated with salinity and density minima in the western gulf, too (Balch et al., 2012). On the eastern side of the GoM, in SS water, there was a similar drop in total alkalinity suggesting that incoming SS water was also being diluted by river discharge (except this would have presumably originated from the Gulf of St. Lawrence (GoSL) discharge (Smith et al., 2001).

Normalization of TA and DIC to a salinity of 32 allowed other patterns to be discerned. Given the circulation patterns of waters in the GOM (Figure 1) then one can assume that the nTA and nDIC properties of SS waters reflect the surface waters entering the GoM from the NW Atlantic while the western water masses (WMCC, ExtEMCC) and central water mass (JB) reflect the seasonal and annual changes happening within the GoM. Homogenization of the waters along the GNATS transect for temperature and salinity is associated with increases in fall mixing (Balch et al., 2004) (see their figure 7) which would presumably increase further during strong winter mixing. For example, compare the impact of the deep mixing event in December 2014 and January 2015 (Figure 2f), the effects of which were experienced Gulf‐wide for all the OA variables regardless of water mass (Figure 17). Seasonal changes show that by late summer, WMCC, EMCC, and JB waters were, on average, ∼83 μmol kg−1 lower for nDIC over about a 6‐month time period. Moreover, nTA was, on average, 57 μmol kg−1 lower in WMCC, EMCC, and JB waters over the same 6‐month time period. If we use a nTA:nDIC ratio of 2:1 for calcification (Bates et al., 1996b), and correction for changes in alkalinity due to nitrate uptake (Broecker & Peng, 1982), calcification probably accounted for about ∼37% of DIC drawdown in the GoM. Regional calcification can be calculated by extrapolating the average seasonal drawdown of 57 μmol kg−1 nTA over the average 50 m euphotic zone of the GoM and the 90,700 km2 area, to get a regional GoM calcification of 6.32 Tg C y−1. The GOM net community production (NCP) can also be estimated by taking the annual nDIC drawdown, reducing it by the fraction associated with calcification and extrapolating it over the 50 m euphotic depth and aerially over the GoM. This gave a GOM‐wide NCP of 5.85 Tg y−1.

A previous GoM estimate of net primary productivity (NPP) for the Gulf of Maine—based on a harmonic regression technique of field‐based 14C data from 1977 to 1982, horizontal advection and vertical diffusion of nitrate, and potential new production (Bisagni, 2003)—was 37.9 Tg y−1. Another completely independent estimate of GOM NPP—based on GNATS 14C data, applied to a primary productivity model—was ∼38.1 Tg y−1 (Balch et al., 2008), within 2% of the Bisagni estimate. GNATS data show that the mean GoM NPP has dropped ∼50% since 2008 (Figure S12 in Supporting Information S1) which would adjust the 2018 gulf‐wide NPP levels to ∼19 Tg y−1. The seasonal nitrate drawdown in GNATS averages ∼10 μM which happens over about the same 6‐month period as above (not shown). Assuming a Redfield molar ratio of C:N (6.625) (Redfield et al., 1963) then the regional new carbon production based on this nitrate drawdown would be 7.21 Tg C y−1. To be consistent with the same overall NPP estimated by Bisagni (2003) and Balch et al. (2008) (as corrected above by 50%), then the observed nitrate drawdown cited here, would require an average GoM F ratio of 0.38, in good agreement with the estimated F ratio of 0.4 cited by Townsend (1998). Finally, if the eastern SS waters reflect a source‐water into the GoM surface layer (used as an end‐member), the productivity of the GoM drives the waters from being a source of ∼10 μM CO2 to the atmosphere to a neutral/modest sink.

Over the 7‐year time scale that OA variables were measured in GNATS (since 2012), nTA, nTA:nDIC, Ωaragonite, and Ωcalcite have all increased (Figure 18). Salisbury and Jönsson (2018) previously demonstrated that in the Gulf of Maine, during years of marine heat waves, the Ωaragonite responded to not only ocean acidification but also thermodynamic aspects (net heating) and chemical aspects (TA) such that there was a net increase in Ωaragonite even though the pH was decreasing. This is consistent with the positive relationship between Ωaragonite and temperature for GNATS (Figure 19) as well as the increasing trend of Ωaragonite over 7 years as temperature also increased (Figure 18d).

One striking feature of the GNATS OA results was the increase in nDIC and decrease in nTA, pH, Ωaragonite, Ωcalcite (Figure 18), and CO3 = (Figure S17 in Supporting Information S1) during a single winter event which spanned two months, December 2014 to January 2015, and which stretched completely across the GoM. The highly‐significant correlation between temperature and carbonate solubility (Figure 16b), as well as the cross‐gulf response suggest that this probably resulted from deep mixing forced by an atmospheric event, which brought cooler, saltier, more corrosive water to the surface. It also should be noted that alkalinity varied quasi‐conservatively with salinity with the two variables being statistically correlated (Figure S19 in Supporting Information S1; Alkalinity = 51.56 (S) + 523.2; RMSE = 19.22 μequiv kg−1; r 2 = 0.768; F(1, 341) = 1126; P < 0.0001) thus one end member was the deep corrosive water and the other end member was from coastal, river‐impacted waters. Wang et al. (2013) have previously shown that subsurface waters in the Gulf of Maine can show high DIC plus corrosive, low pH and low Ωaragonite values. Our results here show that while such events are not common, surface waters across the GoM can occasionally have Ωaragonite values in the winter <1.6, low enough to be corrosive to aragonitic organisms (Salisbury et al., 2008).

Significant changes in OA variables over the 7‐year period were seen for nTA, nTA/nDIC, CO3 =, and Ωar (Table 2; Figures 18a, 18d, and 18e). Even excluding the data from the deep‐mixing event noted above, nDIC, pH and pCO2 did not show statistically‐significant changes (results not shown). This is probably due to the shorter seven‐years time period that such data have been taken as part of the GNATS, and the significant natural variability inherent to those variables. As rising temperatures and salinities in the GOM suggest increasing influx of NASW, given the co‐variation of salinity and TA (Figure S19 in Supporting Information S1), then we would attribute the increasing nTA and nTA:nDIC ratio, to increasing influx of NASW, as well. There also could have been reduced alkalinity‐drawdown associated with reduced calcification (Table 2; Figure S13 in Supporting Information S1). This increase in nTA has likely influenced the GOM to be a net sink for atmospheric CO2 while the nDIC results show no evidence for the GOM being a net source/sink of CO2. Further, with the increase in nTA by ∼20 μmol kg−1 between 2012 and 2018 (Figure 18a), this geochemical mass balance (ignoring other factors potentially at play) is consistent with the observed decrease in calcification by about half (Figure S13 in Supporting Information S1).

Variability in Phytoplankton Communities, Carbon Fixation, and Factors Controlling Them

In terms of the standing stocks of different generic phytoplankton groups, one of the most striking aspects of the 20‐year GNATS was the transition from the second to the third epoch after 2010. That transition included a strong decrease in coccolithophore concentration (Figure 20c) and the smallest <4 μm particles (Figure 20a). In contrast, the FlowCam data showed strong increases in the 4–12 μm particles and naked/armored dinoflagellates with several years of strong increases in diatoms post 2010 (Figures 20b, 20d, and 20e). There was also a logarithmic steepening of the PSDF slopes from values of −1 to −2 in the central and western Gulf, pre‐2010 (characteristic of more larger phytoplankton species) to gulf‐wide slopes of −3 to −4 (representative of more common nanoplankton‐dominated communities; Figure 20f).

This assemblage transition was roughly concurrent with the arrival of the warm water to depths as deep as 180–250 m (Figures 2 and 9a) as well as higher temperature, salinity and density surface water Gulf‐wide (Figure 10) as observed by others (Friedland et al., 2020; Gonçalves Neto et al., 2021). While this correlation in species trends and hydrography does not indicate causation, the timing is compelling, and could ultimately suggest a connection of the Atlantic Meridional Overturn Circulation to floral changes in the GoM (Gonçalves Neto et al., 2021). Despite these dramatic changes in phytoplankton groups, the change in chlorophyll biomass was not so much observed in the quantity, but instead a change in the spatial pattern, transitioning from elevated concentrations on eastern and western sides of the Gulf and lower concentrations in JB waters (pre‐2010), to relatively constant chlorophyll concentrations across all the water masses of the Gulf after 2010, whether high chlorophyll (e.g., 2016) or low chlorophyll concentrations (e.g., 2017–2018; Figure 12a).

Carbon fixation showed the most profound changes after the beginning of the second epoch (2006), when the primary productivity and calcification dropped to levels <30% of previous values (Figures 15a and 15b). These changes appeared to be associated with an intensification of the hydrological cycle, with several years of anomalously high rainfall during the second epoch (Balch et al., 2012). As of 2018, however, mean carbon fixation had not recovered to values observed prior to 2005 in the ExtEMCC and SS waters; indeed, Gulf‐wide, mean carbon fixation was still about 30% of the mean values in the early 2000s (Figures S12 and S13 in Supporting Information S1). Specific water masses, like WMCC waters, showed less recovery to values observed pre‐2006 (Figures 15a and 15b). Through all of these transitions, there has been little coherent change in the ratio of calcification to photosynthesis. The biomass‐normalized productivity showed little variability among water masses, and but was highest in epoch 1 (100 mg C (mg chl a)−1 d−1), after which it fell to about 30 and has not varied since then, save for a few stochastic anomalously high or low values.

The factors causing these changes in productivity, its efficiency and carbon specific growth rates are, at least qualitatively, related to both nutrients and light. For example, DIN values were drawn down to near zero values during Epoch 1, when productivity and carbon specific growth rates were highest (Figures 10 and 15) and the ratio of DIN:DIP was well below the Redfield ratio of 16:1 (Figure 14a (Redfield et al., 1963);) ultimately suggesting nitrogen‐limited growth. During the period of dramatically‐reduced productivity during epoch 2, nutrient drawdown was less apparent, even in summer months (Figure 10d). Moreover, the DIN:DIP ratios were more frequently elevated, closer to 16, suggesting less evidence during epoch 2 for limitation by DIN relative to phosphate (Figure 14a). Instead, during epoch 2, the concentration of CDOM increased dramatically in the western Gulf, associated with the enhanced runoff (Figure 13a) such that the absorption of light by all the particulate and dissolved material in the blue absorption peak of chlorophyll (Soret band) was as much as 16X greater than the actual absorption by chlorophyll, extending well offshore (Figure 13e) (Balch et al., 2012). This trend has continued since our previous report. Absorption of blue light by CDOM exceeded that by chlorophyll most of the time (except in JB) however, the amount of CDOM during Epoch 2 was particularly pronounced. In epoch 3, when carbon fixation increased somewhat (but nowhere to the level seen in epoch 1), patchy DIN draw‐down was the norm and DIN:DIP ratios were more frequently below Redfield ratio values, again suggesting a resumption of nitrogen limitation (Figure 14a).

Without much knowledge about the factors controlling primary production in the GoM, we gained a more quantitative appreciation of these factors by using a backward stepwise linear regression approach to examine primary production (log10‐transformed to accommodate the log normal distribution for phytoplankton production and biomass variables (Campbell, 1995)). The independent variables that we started with in this analysis were: chlorophyll (log10 transformed; μg L−1), POC (log10 transformed; μM), DIN, silicate, phosphate, N*(μM), DOC concentration, a pg412, a g412, temperature, salinity, b b531, pH, total DIC, ΩAR, bicarbonate, carbonate, and alkalinity. Surface primary production was most correlated to log10 POC, log10 chlorophyll a, plus temperature and N* (Table S8 in Supporting Information S1). These four independent variables accounted for only 36% of the variance in surface primary production, however. The best‐fit prediction expression was: Log10 Prim prod = 0.89 + 0.448 Log10Chl + 0.344 Log10POC – 0.065N* + 0.022 SST (with r 2 = 0.36; all factors significant below P < 0.0001; RMS error of 0.406 log units; ANOVA significant at P < 0.0001) (Table S8 in Supporting Information S1). This relation suggests that while phytoplankton biomass (as both chlorophyll a and POC) is most correlated to primary production, the C:chlorophyll ratio may also play a role given the different coefficients for POC and chlorophyll a, respectively. Moreover, highest primary production values would be associated with negative N* (i.e., negative N* means silicate exceeded nitrate, conditions stimulatory for diatoms). Finally, log10 primary production is positively correlated with temperature, consistent with the observations of phytoplankton growth rates increasing with temperature (Eppley, 1972).

Hydrographic and Optical Connections Through Space and Time

As noted above, the optical properties of the GoM changed significantly over the 20‐years of GNATS time series (Balch et al., 2016a). Overall, the quantity of DOC in the GoM has slightly decreased (Table 2) but the quality of the optically‐active CDOM has shown significantly increased specific UV absorption (SUVA) overall (Figure S20 in Supporting Information S1) suggesting increased aromaticity of the DOC (Kellerman et al., 2018). There is interest in whether DOC concentrations are increasing globally with climate change. In North America, the results are mixed (Monteith et al., 2007). Predictions for the GoM rivers are for DOC concentrations to increase (Huntington et al., 2016). Our results are also mixed in that we saw significant decreases in DOC concentrations in WMCC and ExtEMCC over the 14 years that this variable was measured, no statistical change in JB waters and significant increase in DOC in SS waters (Figure S21 in Supporting Information S1; Table 2). For studies in the NE North America, both Clair et al. (2008) and Fahey et al. (2005) showed decreasing concentrations of DOC associated with watersheds. We suspect that the increasing DOC trends in the SS water represent changes happening outside the GoM, possibly reflecting changes in the GoSL.

The GoM is primarily considered Case 2 dissolved (where inherent optical properties (IOPs) are dominated by materials other than phytoplankton, such as colored dissolved organic matter from rivers) and sometimes Case 2 particulate (Balch et al., 2004) (in which the IOPs are dominated by non‐algal particulate matter such as sediments). It is no surprise then that qualitatively, one of the stronger hydrographic connections to ocean optical properties is seen when comparing surface salinity (Figure 10b) or density (Figure 10c) with a g412 (Figure 13a). In the west, the positions of the WMCC and the ExtEMCC are clearly demarcated by all these variables, especially from 2006 to 2010. Interestingly, the distribution of total particulate scattering (Figure 13b) was not as clearly defined, suggesting that freshwater sources were not as large a source of particulate backscattering in the WMCC and ExtEMCC. Viewing the time/space distribution of b bptot532 reveals that in epoch 3, the period of warm NASW intrusion (Friedland et al., 2020; Gonçalves Neto et al., 2021), the low b bptot values previously observed in just JB waters, extended all the way to 70°W by 2018, concurrent with the lowest chlorophyll values seen in the western Gulf over 20 years, and synchronous with low concentrations of several of the algal classes (<4 μm, 4–12 μm, coccolithophores and diatoms; Figure 20). Note, however, this is different for dissolved scattering. It has been shown previously that total a g412 and total b g412 are well correlated in the GoM (Balch et al., 2016a), so at least for fractions <0.2 μm, river‐supplied scattering materials can be important (Zhang & Gray, 2015).

The elevated optical CDOM signal in the western GOM during 2006–2010 (Figure 13a) also showed a strong DOC maximum (Figure 12e), however the converse was not always true, that is, the strong DOC maximum in the SS water during 2016–2018 only showed a minor enhancement in CDOM. Thus, there is a fundamental disconnect between the CDOM and DOC composition of the water on the east side (SS water) and west side (WMCC and ExtEMCC) of the GoM, simply belying their different primary sources, GoSL and Maine watershed, respectively. The differences also are due to their history. CDOM‐rich water coming from the GoSL takes about a year to flow to the Gulf of Maine, and photooxidation likely will decrease the CDOM during this period (Smith et al., 2001). In contrast, water from the Maine watershed takes about a month to be transported to the WMCC and ExtEMCC (Balch et al., 2016b) so the CDOM likely still remains highly blue‐absorbing.

The slight (but significant) decrease in b b’ over the 20‐years GNATS is consistent with the decreases in the concentrations of PIC (Figure S7 in Supporting Information S1), coccolithophores (Figure S22 in Supporting Information S1) detached coccoliths (Figure S23 in Supporting Information S1; Table 2), as well as the decrease in calcification (Figure S13 in Supporting Information S1; Table 2). Nonetheless there was no significant change in the fraction of backscattering contributed by PIC in the GoM over the two decades. The mean value of this fraction was about 9% (95% CI of 1%–70%). The SS waters still remain the region with the largest amount of acid‐labile backscattering on average (Figure 13d).

Concluding Remarks

This 20‐year retrospective of the GNATS has shown a number of highly significant trends over time in the GoM, trends that can be grouped into several categories: (a) hydrographic trends in temperature, salinity, and density show overall net increases (seen both in GNATS and NERACOOS data at buoy M) which we suspect are consistent with deep intrusions of NASW into the GoM (as opposed to LSW or SS water). We note that the 20‐years surface temperature trends were curvilinear—initially decreasing from 1998 to 2010 then steadily rising thereafter—but nonetheless an overall net temperature increase over 20 years. (b) A multiple‐linear stepwise correlation approach of many variables within the GNATS data set show that the overall decrease in algal primary production was most significantly correlated to decreases in standing stocks of carbon and chlorophyll, increases in N*, and temperature‐dependent growth, consistent with the intrusions of warm NASW. The significant decrease in GoM primary production will likely have negative impacts on the overall productivity of the higher trophic levels as well as fishery yields (Moore et al., 2018). However, the exact nature of the negative impacts on fishery yields are likely to vary over different marine ecosystems as a function of factors such as latitude, phytoplankton biomass, particle export and mesozooplankton productivity (Friedland et al., 2012; Sherman et al., 2013). (c) While surface DIN showed no significant change over the 20 years, there were other mixed (but nonetheless statistically significant) trends in nutrients and their relative variability; concentrations of phosphate increased, silicate decreased, residual nitrate increased and DIN:DIP decreased. (d) while DOC has slightly decreased gulf‐wide, the optically‐active fraction of DOC, CDOM, has increased in the GoM over the 20 years: (a pg440:a ϕ440, a g254, and SUVA254), suggesting more land‐to‐sea transfer of aromatic humic materials which had previously been associated with changing precipitation patterns (Balch et al., 2016b). (e) Finally, despite some dramatic seasonal drops in Gulf‐wide OA‐related variables such as pH, CO3 =, Ωar, Ωca, associated with deep winter mixing of corrosive water to the surface (e.g., winter 2015), the overall trend in these variables was nonetheless positive, again associated with the deep intrusions of warm, salty NASW. If the intrusions of NASW into the GoM (measured over a 25 years altimetric record) are indeed associated with changes in the AMOC (Gonçalves Neto et al., 2021; Record et al., 2019), then our results suggest that the changes in the hydrography, chemistry, biology and bio‐optics of the GoM are fundamentally tied to changes well outside the GoM, with changes of circulation of the warming North Atlantic Ocean.

Dedication

We dedicate this paper to Dr. George Aiken, USGS, Boulder, CO (1951–2016) who provided the critical expertise and analyses of DOM samples over the GNATS time series, fundamentally allowing us to connect the terrestrial carbon cycle to the marine carbon cycle in the Gulf of Maine.

Conflict of Interest

The authors declare no conflict of interest relevant to this study.

Supporting Information S1

Click here for additional data file.