Thirty-Three Years of Ocean Benthic Warming Along the U.S. Northeast Continental Shelf and Slope: Patterns, Drivers, and Ecological Consequences.

The U.S. Northeast Continental Shelf is experiencing rapid warming, with potentially profound consequences to marine ecosystems. While satellites document multiple scales of spatial and temporal variability on the surface, our understanding of the status, trends, and drivers of the benthic environmental change remains limited. We interpolated sparse benthic temperature data along the New England Shelf and upper Slope using a seasonally dynamic, regionally specific multiple linear regression model that merged in situ and remote sensing data. The statistical model predicted nearly 90% of the variability of the data, resulting in a synoptic time series spanning over three decades from 1982 to 2014. Benthic temperatures increased throughout the domain, including in the Gulf of Maine. Rates of benthic warming ranged from 0.1 to 0.4°C per decade, with fastest rates occurring in shallow, nearshore regions and on Georges Bank, the latter exceeding rates observed in the surface. Rates of benthic warming were up to 1.6 times faster in winter than the rest of the year in many regions, with important implications for disease occurrence and energetics of overwintering species. Drivers of warming varied over the domain. In southern New England and the mid-Atlantic shallow Shelf regions, benthic warming was tightly coupled to changes in SST, whereas both regional and basin-scale changes in ocean circulation affect temperatures in the Gulf of Maine, the Continental Shelf, and Georges Banks. These results highlight data gaps, the current feasibility of prediction from remotely sensed variables, and the need for improved understanding on how climate may affect seasonally specific ecological processes.

Introduction

Ocean habitats are warming at an accelerating rate, with the past two decades accounting for over half of the increase in ocean heat content that has occurred since preindustrial times (Gleckler et al., 2016). Along the U.S. Northeast Continental Shelf, several studies have documented multidecade or interannual increases in sea surface temperature based on ship‐based records (Forsyth et al., 2015), shore‐based measurements (Fulweiler et al., 2015; Oczkowski et al., 2015; Rheuban et al., 2016), remote sensing (e.g., Chen et al., 2014; Maynard et al., 2016; Pershing et al., 2015), and multiplatform reconstructions (Friedland & Hare, 2007; Shearman & Lentz, 2010). These studies document both long‐term changes (Nixon et al., 2004; Friedland & Hare, 2007) as well as accelerated warming in the past decade (Forsyth et al., 2015; Pershing et al., 2015; Rheuban et al., 2016), the latter associated with both atmospheric (Chen et al., 2014) and ocean circulation (Chen et al., 2015) changes. Regional forcing may lead to variation in rates of warming; for example, sea surface temperatures in the Gulf of Maine warmed faster than 99% of the ocean over the past decade (Pershing et al., 2015).

Ocean warming over the last several decades has affected organism abundances, distributions, and phenology as well as trophic transfer, disease incidence, and biogeochemistry (Doney et al., 2012; Nye et al., 2009; Pershing et al., 2015). In the North Atlantic, warming has been associated with declines in cold water species and increases in warm water species (Gawarkiewicz et al., 2013; Montero‐Serra et al., 2015). Phytoplankton blooms are occurring earlier (Ji et al., 2010; Hunter‐Cevera et al., 2016) potentially leading to trophic mismatch (Edwards & Richardson, 2004) and reduced trophic transfer. Distribution shifts have been noted from zooplankton to fish (e.g., Friedland et al., 2013; Genner et al., 2010; Pinsky et al., 2013) with species range extents moving poleward.

Benthic and demersal species such as American lobster, sea scallop, and Atlantic cod play important ecological and economic roles in the New England Shelf system (e.g., Hare et al., 2016; Pershing et al., 2015). Benthic species comprise over 92% of the landing value of northeast regional commercial fisheries (1.23 billion USD in 2015, NEFSC; Gledhill et al., 2015); thus, understanding species vulnerabilities to changing ocean conditions is a high priority for resource monitoring and management (e.g., Hare et al., 2016; Stortini et al., 2015). Temperature is a dominant driver for predicting pelagic fishery health and likely a strong driver in the benthic environment. Gawarkiewicz et al. (2012) reported that during the winter of 2011, benthic temperature anomalies from ship‐based measurements exceeded +4°C over Georges Bank. These warm temperature anomalies were associated with northward shifts in the distribution of Atlantic cod (Friedland et al., 2013) and changes in the phenology of American lobster that dramatically impacted fisheries revenue (Mills et al., 2013). Physiological thresholds and tolerance of temperature variability are likely to vary by species and population. Furthermore, mechanistic drivers of temperature shifts are likely not uniform over the domain, and thus the capacity to develop short‐term forecast models or early warning systems for management is limited.

While regional variability of surface temperature signatures is discernible from remote sensing, the spatial patterns of benthic temperature variability are less well characterized from in situ observations. Benthic temperatures may be decoupled from surface values because of physical transport or stratification. The NE Shelf is a western boundary confluence zone, with Labrador Current/Scotian Shelf waters transporting cold, fresher water south and the Gulf Stream transporting warm, saline water north. Equatorward flow is strong along the coast of Gulf of Maine and weakens toward the mid‐Atlantic Bight. Cross Shelf heat transport can occur through Gulf Stream rings (Joyce et al., 1992) as well as Shelf break processes and frontal meandering (Gawarkiewicz et al., 2004). In the spring and summer, surface warming due to solar insolation, reduced heat exchange to atmosphere, and freshwater inputs can result in strong vertical stratification until winter mixing and cooling destratifies the water column. In the Gulf of Maine and Slope regions, winter temperature gradients are inverted with cold, freshwater overlying warmer, saline water, with the inverted or weak vertical stratification supported by salinity (Deese‐Riordan et al., 2009). Thus, spatial and seasonal variability challenge the capacity to predict benthic temperatures from readily available satellite products or surface underway data.

An understanding of the spatial and temporal variability of bottom water temperatures will be critical to appropriately managing ecosystems in changing environment. Such an understanding can also illuminate vulnerabilities and inform improved sampling to reduce uncertainties. Here we utilize a statistical approach to integrate remote sensing and in situ observations, estimate uncertainties to determine where and when statistical models can predict or interpolate between sparse data points, determine spatial variability of temporal trends, and illuminate dominant interannual forcing factors. After demonstrating the skill of the statistical model in space and time, we ask the following questions relevant to climate forcing and benthic ecology:

How have benthic temperatures changed along the NE Shelf over the past 34 years?

What are the patterns of spatial and seasonal variability in rates of benthic warming?

Finally, how do the major drivers of interannual benthic warming (e.g., surface warming, increased winds, and circulation) change regionally?

Methods

CTD, Bathymetry, and Wind Data

The domain of our study was limited to the Continental Shelf and upper Slope from Cape Hatteras to the Gulf of Maine and depth less than 500 m. CTD casts from 1982 to 2015 were obtained from the Northeast Fisheries Science Center's Oceanography Branch. High‐resolution CTD casts extracted from the World Ocean Database (WOD; Boyer et al., 2013) were used to supplement the NEFSC data. The WOD data were filtered such that there were no duplicate CTD casts. Data were binned in both time and space by averaging all profiles by month onto a 1/3° block‐averaged grid, with greater spatial and seasonal density occurring after 1990. Spatial variation of sampling density varied, with the greatest coverage occurring on Georges Banks (∼50% of months sampled) and data density limited south of Cape Hatteras (∼10% of months sampled; Figure 1). Within grid cells, the in situ data record varied in length, with coverage ending in 2014 for most of the domain, but some nearshore records ending before 2013 and the northern Gulf of Maine records persisting until 2015 (supporting information Figure S1). We omitted records that ended prior to 2013 and truncated records to 2014.

Figure 1

Data density for in situ benthic temperatures from 1982 to 2012. Perfect coverage would include data from 372 months in each pixel (1/3° × 1/3°). White circles indicate regions where statistical model was validated through time.

Bathymetric data were obtained from the National Geophysical Data Center, U.S. 3 arc sec Coastal Relief Model (CRM) volumes 1 and 2 and binned to the same grid as the CTD casts. Surface temperatures were obtained from the binned casts as the average of the first 2 m of the profiles. Bottom temperatures were assumed to be the last measurement in the profile if the pressure reported in the CTD cast was within 25 m of the bathymetry from the CRM. Because we were interested in nearshore processes, we limited our analysis to only CTD profiles that were in waters shallower than 500 m. Wind stress was calculated from monthly mean vector wind data from the North American Regional Reanalysis for the New England Shelf (Mesinger et al., 2006). For this analysis, we only consider boundary layer mixing and net heat surface exchange, not wind stress curl and subsequent Ekman pumping. Wind stress data are not normally distributed and were log10‐transformed prior to analysis.

Grid cell‐specific seasonal climatologies for surface and benthic temperature were calculated using in situ data from 1991 to 2010. This span covered a period of relatively high and spatially uniform seasonal data density but excluded recent warming (i.e., 2012; Chen et al., 2014). Missing spatial data were interpolated using a nearest neighbor technique restricted to the adjacent pixel; thus some, but not all spatial gaps were filled. The seasonal cycle was then calculated in a two‐step process by subtracting the long‐term variability (Loess filter, 24 month span), and then adding the long‐term mean of the smoothed data to the anomalies. In some regions (e.g., Gulf of Maine), observed interannual variability exceeds that of the seasonal cycle; the Loess filter processing was necessary to account for the within year variability appropriately. A smoothed (3 month running mean) seasonal cycle was fitted per pixel, excluding cells with <6 calendar months represented to avoid seasonal biasing.

We used the NOAA 1/4° daily Optimum Interpolation Sea Surface Temperature (OISST) (Reynolds et al., 2007) in our predictive model. The OISST record is constructed by combining observations from different platforms (satellites, ships, and buoys) on a regular global grid that minimizes sensor and spatial biases by tuning to in situ data. Preliminary comparisons (supporting information Figure S3) suggested that the AVHRR 5.2 SST data (Casey et al., 2010) introduced a cold bias in our region, a common problem associated with imperfect cloud masking. This cold bias was minimized using the assimilated OISST product. OISST daily data were downloaded, and monthly means were computed per pixel.

Climate Drivers

We investigated the effect of three drivers of interannual variability, the North Atlantic Oscillation (NAO; Hurrell, 1995), the Atlantic Meridional Oscillation (AMO; Enfield et al., 2001; Trenberth & Shea, 2006), and the Pacific Decadal Oscillation (PDO; Mantua et al., 1997). During positive NAO phases, there is a northward shift and increase in westerly winds and cooler temperatures off of Labrador and northern Newfoundland but warmer temperatures off of the mid‐Atlantic U.S. The position of the Gulf Stream shifts north during positive NAO phases, although the relationship appears to have weakened over the recent decade (Meinen et al., 2010) and the southward transport of cold, Labrador Subarctic Slope Water (LSSW) is low. Conversely, during a negative NAO phase westerly winds decrease and shift southward. Water temperatures are warmer off of Labrador and Newfoundland but cooler off of the eastern U.S. LSSW penetrates to the mid‐Atlantic Bight, displacing Atlantic Temperate Slope Water.

The AMO signal is based on spatial patterns in SST variability after removing the effects of anthropogenic forcing on temperature, revealing natural long‐term variability patterns in SST. The AMO is thought to be related to variations in the North Atlantic branch of the deep thermohaline circulation (Wang & Zhang, 2013). The Atlantic Meridional Overturning Circulation, or AMOC, is a large‐scale overturning circulation that transports heat from the South Atlantic to the North Atlantic high latitudes. Anomalies in AMOC can lead to a surface temperature anomaly pattern that mimics that of the AMO (Wang & Zhang, 2013).

The PDO is defined by leading empirical orthogonal function (EOF) of monthly sea surface temperature anomalies over the North Pacific. The PDO consists of a warm and cool phase that alter upper level atmospheric winds; teleconnections to the Atlantic occur as a result of the propagation of atmospheric Rossby waves. Shifts in the PDO phase thus can have significant implications for global climate, affecting Pacific and Atlantic hurricane activity, SST, and precipitation. PDO phases are associated with opposite winter conditions in the Pacific Northwest than in the eastern United States where negative phases are associated with warmer winter temperatures and decreased precipitation along the southeast and mid‐Atlantic coasts. Recent efforts have suggested that these patterns extend into the Gulf of Maine also (Pershing et al., 2015).

Regionalization

Grid cells were aggregated into seasonally dynamic classes and statistical analyses were conducted in a nested fashion. As a first step, the normalized vertical gradient (GRAD) was calculated across monthly climatological temperatures (absolute temperature, T, in K): with x and y corresponding to the longitude and latitude coordinates of the grid cell and m corresponding to the climatological month.

Dynamic regions were classified from two variables, normalized vertical gradient monthly climatology and bathymetry, using a probabilistic self‐organizing mapping algorithm (Kavanaugh et al., 2014) that groups data in space and time simultaneously. Six regions explained ∼75% of the variance of the input data (supporting information Figure S2), and the move from five to six regions resulted in a step increase in variance explained. The next apparent step increase was at nine regions; however, there were insufficient observations in each class at nine regions to populate some of the regression models. Regions represent widely varying conditions from well mixed (gradient ∼0), to stratified (gradient >0), and thermally inverted (gradient <0) as well as water columns in deeper (>200 m) and shallower water (Figures 2b and 2c). Region 1 described shallow and well‐mixed conditions that were present on Georges Bank year round and in the mid‐Atlantic Bight during fall and winter. Region 2 described thermally inverted conditions in the deeper waters of the Gulf of Maine and the Continental Slope that occur through the winter. Region 3 described the stratified conditions present in the same geographic areas as Region 2 during the summer months and early fall. Region 4 also described a stratified water column during summer, but in moderately shallow areas, primarily in the mid‐Atlantic Bight. Region 5 describes thermally inverted conditions that occur in shallow regions in the nearshore and mostly southern part of the domain during the winter. Finally, Region 6 describes a variable water column, but is isolated to very shallow nearshore regions and embayments.

Figure 2

Seascape classification used to partition regression analysis by dynamic region. (a) Climatological seasonal progression of classes; (b) median of thermal stratification values (normalized vertical temperature gradient, equation (1); and (c) bathymetry across seascape regions. Blue lines denote 25th and 75th percent quantile, with whiskers representing 1 and 99%. Outliers are shown in red.

Multiple Linear Regression

A multiple linear regression model was built using in situ data and then subsequently evaluated using satellite SST in lieu of in situ SST. First, each data point was ascribed a region category based on grid cell and month, so that each month's classification in a grid cell matched that of the climatological month's classification (Figure 2). Then within each of the six stratification classes, benthic temperatures were predicted for each grid cell and time, BTx,y,t, as a function of monthly surface temperature (SST), wind stress (tau), bathymetry (Z), the climatological seasonal benthic temperature (BT x,y, ), the climatological seasonal normalized gradient (GRAD), AMO, NAO, and PDO using stepwise multiple linear regression with interactions:

Model parameters were chosen based on minimizing root‐mean‐square error (RMSE), with sequential inclusion (p < 0.05) and removal of variables (p ≥ 0.2). After the best model was fit using in situ sea surface temperatures as part of the predictor suite, benthic temperatures were estimated using monthly means of satellite SST in place of in situ SST. A merged data set was created by substituting the predicted benthic temperatures on a per grid cell basis where in situ data were missing, interpolating between in situ observations.

To assess how the coefficient estimates were being inflated by multicollinearity in the benthic temperature models, we calculated the variance inflation factors (VIFs; O'Brien, 2007; Saba et al., 2015) for each predictor variable used: where is the variance explained in a linear regression model with variable (i) used as a response and all other predictors used as explanatory variables.

VIFs were calculated within and across all regions (supporting information Table S1). A typical rule of thumb is that multicollinearity can be ignored for a particular predictor variable if its VIF is less than 4 (O'Brien, 2007) which all variables met.

Trends and Drivers: Interannual variability was calculated by removing the mean seasonal cycle from the merged data set (in situ data and statistical model) for each grid cell, and then regression analyses were conducted on the anomalies to quantify trends and characterize possible drivers. Linear fits were calculated as a function of time using all data or only December through March for winter trends. Spatial patterns of slopes (p < 0.1) were mapped and compared. Drivers and response of benthic temperature anomalies were assessed using multiple linear regression of the standardized anomalies on a per‐grid‐cell basis. Response and driver variables were standardized by subtracting the mean and dividing by the standard deviation (z‐scores), with driver variables including surface temperature anomalies, wind stress anomalies, and climate indices (NAO, AMO, and PDO). After standardization, regression slope coefficients represent relative influence or effect sizes of a given driver after considering the effect of the other drivers in the regression model. Effect sizes were then compared among each other and spatially for strength of influence.

Results

Regional Validation, Time Series

Our statistical model explained nearly 90% of the variation of observed in situ benthic temperature data providing a viable interpolant across regions with varying thermal profiles (Figure 3 and Table 1).

Figure 3

Calibration and validation of the statistical model for (a) in situ data (training set) from 1982 to 2012; (b) validation using in situ data not included in original statistical model; and (c) statistical model applied to OISST data. Colors denote different seascape regions as shown in Figure 2.

Table 1

Multiple Linear Regression Models to Predict Benthic Temperatures Within Dynamic Regions

	Region 1		Region 2		Region 3		Region 4		Region 5		Region 6
	EST	SE	EST	SE	EST	SE	EST	SE	EST	SE	EST	SE
Intercept (°C)	0.05	0.14	3.08	1.36	−5.10	2.19	–4.17	4.49	10.29	4.42	8.86	6.22
SST (°C)	0.20	0.07	1.04	0.20	0.08	0.11	–0.23	0.18	0.54	0.44	1.71	0.53
bathy (m)	1.03	0.01	0.15	0.47	0.93	0.94	2.63	1.98	–3.36	2.17	–2.41	2.51
seas_bt (°C)	229.09	23.14	–0.36	0.24	1.32	0.14	1.45	0.25	0.39	0.55	–1.05	0.58
seas_grad	0.05	0.05	93.85	76.77	57.97	42.98	74.17	78.22	–80.04	54.94	−804.77	187.42
tau (N m−2)	0.58	0.35	1.40	0.47	−1.93	0.86	–1.32	1.39	6.25	2.01	2.35	2.30
NAO	3.19	0.96	−1.55	0.53	–1.24	0.69	0.68	0.78	2.18	0.84	–0.05	1.82
AMO	0.26	0.10	–1.48	2.03	−3.51	1.61	−4.32	0.94	2.18	3.06	–0.90	2.12
PDO	0.05	0.14	0.06	0.29	–0.26	0.22	−0.36	0.16	0.58	0.42	–0.12	0.80
SST:bathy	–0.12	0.07	−0.34	0.09	0.09	0.05	0.15	0.07	−0.73	0.12	–0.34	0.21
SST:seas_bt	0.03	0.01	−0.01	0.01	–0.02	0.01	−0.02	0.01	−0.02	0.01
SST:seas_grad	−10.66	2.78	−14.72	2.09	−1.79	0.80					−12.88	4.87
SST:tau							−0.07	0.04	−0.56	0.16	0.33	0.19
SST:NAO	−0.07	0.03	0.19	0.05	−0.08	0.02	0.08	0.03			−0.27	0.15
SST:AMO			0.75	0.16					0.38	0.20
SST:PDO			−0.05	0.02					−0.15	0.04	0.22	0.06
bathy:seas_bt			0.34	0.09					0.69	0.15	0.43	0.23
bathy:seas_grad	−82.98	11.43	−72.48	30.87	−42.42	19.77	−96.33	36.21	63.24	36.11	129.28	74.95
bathy:tau					0.47	0.38	0.96	0.60	–1.77	0.95
bathy:NAO	−0.60	0.14			0.38	0.29			−1.45	0.57	1.00	0.77
bathy:AMO	−0.92	0.37	1.05	0.69	1.22	0.62			–2.13	1.42	2.51	1.82
bathy:PDO			–0.21	0.13	0.16	0.10
seas_bt:seas_grad	−5.34	0.91			3.77	2.50	8.00	4.04			14.71	4.66
seas_bt:tau			−0.18	0.06	0.15	0.05	0.12	0.07	0.26	0.20	−0.41	0.20
seas_bt:NAO	−0.06	0.01	−0.19	0.06	0.09	0.05					0.26	0.16
seas_bt:AMO	−0.28	0.03	−1.24	0.21	−0.21	0.09			−0.53	0.23
seas_bt:PDO	0.01	0.01	0.08	0.02					0.07	0.05	−0.16	0.06
seas_grad:tau			−58.88	16.80							−247.39	59.38
seas_grad:NAO	−15.39	5.78	−75.69	16.86	26.83	10.06	–27.47	15.42	65.68	26.00	60.12	40.89
seas_grad:AMO	−82.16	13.74	−324.38	54.70			126.98	24.40
seas_grad:PDO	7.33	2.46			−8.31	2.31	6.27	4.17			−38.76	18.59
tau:NAO	−0.54	0.10	−0.58	0.21			0.70	0.24
tau:AMO	−0.82	0.28	−1.30	0.64	−0.93	0.36
tau:PDO	0.30	0.05									0.44	0.22
NAO:AMO	–0.24	0.14			−0.64	0.23
NAO:PDO	−0.24	0.03	–0.13	0.09			0.22	0.07			0.59	0.21
AMO:PDO	−0.72	0.10	−1.43	0.26							1.63	0.56
N	8,286		1,961		2,096		1,683		214		1,054
RMSE	1.16		1.38		1.05		1.72		2.21		2.26
R2	0.92		0.66		0.61		0.51		0.85		0.75
RMSE (full domain)	1.5		1.52		1.05		1.76		3.24		2.58
R2 (full domain)	0.84		0.54		0.5		0.45		0.64		0.65

Note. Listed are estimates of regression coefficients (EST) and their standard error (SE). Final model for each region is a result of a stepwise forward and backward elimination process to optimize model root‐mean‐square error (RMSE). Coefficients are significant (p < 0.05) unless italicized. Units are omitted from interaction terms (:) for clarity. Variance explained (R2) and RMSE are reported for individual models within dynamic regions and for the best fit model across the domain.

The statistical model varied in complexity (i.e., number of coefficients) depending on bathymetry and whether the stratification class was thermally inverted, well mixed or stratified (Figure 2 and Table 1). Model complexity and skill (Table 1, R2 = 0.93) was highest in the well‐mixed case (Region 1), where benthic temperatures were related to SST, bathymetry, wind stress, the NAO, and the PDO. In the thermally inverted cases (Regions 2 and 5, R2 = 0.74 and 0.85, respectively), bottom temperatures were correlated with surface temperatures and wind stress; climate indices only affected benthic temperatures when interacting with other variables. The model had moderate skill in predicting bottom temperatures over thermally stratified summer conditions in both deep and very shallow water (Regions 3 and 6, respectively; R2 = 0.71 and R2 = 0.74). Region 3 bottom temperatures were affected by SST, the seasonal gradient, wind stress, both the NAO and the AMO, but not bathymetry except in the interaction of bathymetry with wind stress and the climate indices. Region 6 bottom temperatures were affected by SST, the seasonal mean and gradient, bathymetry, interactions between the NAO and SST, wind stress and PDO, and interactions between the NAO and AMO with the PDO. Despite relatively high data density (N > 1,500), the model had the least skill predicting in thermally stratified conditions over the mid‐Atlantic Bight (Region 4, R2 = 0.51), with minimal effects of SST, but significant roles of other variables including wind, NAO, AMO, PDO and interactions.

In general, classification improved skill and reduced RMSE by 0.13°C across all regions. The improvement was most dramatic in Region 5 (inverted temperature gradient, nearshore mid‐Atlantic in winter), where RMSE was reduced by nearly 1°C. Here, data density (∼220 observations) was less than 1/5 that of any other region, and less than 1/40 of Region 1 (>8,000 observations), where model skill was highest. All regions benefited by having a climatological benthic temperature or gradient seasonal cycle.

Despite being trained on data that retained a seasonal cycle, the model was able to recover interannual variability across all regions and months (actual anomalies versus statistical benthic temperature model anomalies: R2 = 0.68). However, the model failed to capture the amplitude of the interannual signal in the Gulf of Maine (Figure 4a). The model also failed to capture some episodic highs at Cape Hatteras, possibly caused by warm core rings impinging upon the Shelf (Figure 4c). With these regional caveats in mind, we used the statistical model as an interpolant between measured data in time in the analyses of the spatial patterns of interannual to decadal trends.

Figure 4

Validation of satellite‐based predictive model against monthly time series data at individual sites. (a) Gulf of Maine; (b) Georges Bank; and (c) Cape Hatteras (see white circles in Figure 1).

Spatial Trends

Benthic warming, calculated from 1982 to 2014, was evident throughout the region (Figure 5a), with highest rates associated with inshore and nearshore regions (e.g., Chesapeake Bay and Delaware Bay: nearly 0.6°C/decade) and George's Bank (0.3°C/decade). Despite its depth, warming was also evident in the Gulf of Maine, albeit to a lesser degree (0.2°C/decade). Benthic temperatures increased throughout the domain with the exception of south and east of Cape Hatteras following the general pattern of surface warming (Figure 5b). Rates of benthic warming generally were between 0.1 and 0.2°C/decade whereas surface warming rates generally were between 0.2 and 0.4°C/decade, although exceeded 0.5°C/decade in the Gulf Stream and along the continental margin. The only exception to the pattern of warming occurred to the south and east of Cape Hatteras, where significant cooling was apparent in surface and, to a lesser degree, benthic waters.

Figure 5

Interannual trends of temperature from 1982 to 2014 for (a) benthic environments within the 500 m depth contour and (b) surface water. Rates are °C yr−1. Trends (p < 0.1) are shown.

Benthic warming rate was greater and on average 1.3 times faster in winter than over the entire year (Figure 6a). The opposite was true in the surface (Figure 6b), where rates of warming were much slower in winter than the rest of the year. On Georges Bank and through southern New England, rates of benthic temperature rise, particularly in winter, exceeded that seen in the surface (Figure 7), outpacing surface warming by a factor of 1.5 to 2.

Figure 6

Ratio of wintertime (December–March) to rest of year warming rates from 1982–2014 for (a) benthic and (b) surface habitats (dimensionless). Trends (p < 0.1) are shown.

Figure 7

Ratio of benthic to surface warming rates from 1982 to 2014 for (a) winter and (b) all seasons (dimensionless). Trends (p < 0.1) are shown.

Interannual Drivers of Benthic Temperature Patterns (SST, Winds, AMO, PDO, and NAO)

SST increased throughout the domain with the exception of south of Cape Hatteras, where surface cooling was evident (Figure 5b). Wind stress generally increased in the southern and central domain, but decreased in the northern domain, including the northern Gulf of Maine (Figure 8a). The AMO has been in a positive, warm phase since the early 1990s (Figure 8b). The NAO, however, exhibited large interannual variability over the course of our study (Figure 8c). Negative excursions, which favor a warmer Gulf of Maine relative to the mid‐Atlantic, occurred in 2000, 2010, and 2012 and were generally stronger in winter. The PDO also exhibited large interannual variability with brief, negative phases (associated with warmer conditions on the U.S. east coast occurring) in 1991 and 1995, and longer negative phases in 2000, 2008–2009, and 2011–2013.

Figure 8

Temporal trends of (a) wind stress (N m−2 yr−1), (b) AMO, (c) NAO, and (d) PDO from 1982 to 2014.

Sea surface temperature was the dominant correlate of interannual variability of benthic temperatures, with effect sizes contributing to greater than 0.5 standard effect size (unitless) in the northern Gulf of Maine (Figure 9e). This indicates that for every 1 standard deviation rise in surface temperatures, benthic temperatures rose by 0.5 standard deviations in that grid cell. Wind stress had weak, patchy, and mixed effects on bottom temperatures, with positive effects of wind stress in central Gulf of Maine and negative effects in some shallow and slope waters (Figure 9a). We note that SST and wind stress were used in the interpolative model. However, the sign or significance of the SST and wind in the interpolative model varied as a function of the seasonal stratification classification (Table 1) whereas on interannual scales, the SST effect is consistently positive (Figure 9e) and the wind effect is spatially variable (Figure 9a). The AMO was also positively correlated with benthic temperatures, although its influence was limited to shallower regions on Continental Shelf and Georges Bank (Figure 9b). This is consistent with modeling efforts that suggest that the warm phase of the AMO is linked to a stronger and northward Gulf Stream, although there may be enhanced mesoscale activity or other physical processes responsible for bringing Gulf Stream water from the Shelf break to the shallow coastal waters. The NAO had moderate effects and followed the canonical NAO pattern, with negative NAO associated with increased temperatures in the Gulf of Maine and neutral to positive phases associated with warmer temperatures in the mid‐Atlantic south to Cape Hatteras (Figure 9c). Benthic temperatures were negatively correlated with the PDO throughout the domain (Figure 9e).

Figure 9

Standardized effect sizes of different drivers of benthic temperature as determined by multiple linear regression. Drivers included (a) wind stress, (b) AMO, (c) NAO, (d) PDO, and (e) SST. Both drivers and response (benthic temperature) were standardized using z‐scores prior to analysis; effect sizes (slopes, p < 0.05) are unitless.

Discussion

Continental shelf ecosystems in the temperate North Atlantic have undergone rapid warming, with increases in average surface temperatures of up to 1.3°C in the last three decades (Pachauri, 2007; Shearman & Lentz, 2010; Sherman et al., 2009). Benthic temperatures in this domain have also risen over the past 34 years, and likely will continue to rise for at least the near‐ to mid‐term future (this study, Maynard et al., 2016; Rheuban et al., 2017; Saba et al., 2016). Although the rate of benthic temperature increase on the New England Shelf is generally less than that of the surface (0.2°C/decade or 0.68°C over the last 34 years), some benthic environments either meet (George Bank, 1.3°C in 30 years) or exceed rates of surface warming (off of Delaware and Chesapeake Bays). Warming in these regions is most evident during the winter, where benthic temperatures are increasing on average 1.6 times faster than the rest of the year. Thus, management strategies of marine resources must consider the effect of increased benthic temperatures overall, but also where and when rates may be inferred from surface warming, and the effects of seasonal asymmetry on community dynamics.

Warming can trigger changes in abundance, range, phenology and body sizes of marine organisms (e.g., Doney et al., 2012; others below), including commercially important species. While studies have demonstrated range shifts on pelagic fisheries (Montero‐Serra et al., 2015; Pinsky et al., 2013), evidence for demersal fish community shifts are somewhat more limited (Genner et al., 2010; Simpson et al., 2011) and have generally focused on benthic environments as refugia from surface warming. However, benthic warming will likely lead to lower recruitment, suboptimal growth conditions and reduced fishery productivity of several species of invertebrates and fish (Hare et al., 2016), including Atlantic cod (Drinkwater, 2005; Fogarty et al., 2008; Planque & Frédou, 1999; Pershing et al., 2015). Increased benthic temperatures may exacerbate the risk of predation as cod move to even deeper water (Linehan et al., 2001) and increase energy demands as the fish transitions from benthic to pelagic prey (Sherwood et al., 2007). For the high value American Lobster fishery, Wahle et al. (2015) found significant declines in the abundance distribution and size composition of juveniles associated with increased temperatures in Narragansett Bay. Habitats previously considered suitable are now consistently surpassing a physiological threshold of 20°C (McLeese & Wilder, 1958; Pearse & Balcom, 2005). While some variability between northern and southern populations may minimize effects on some parts of early life history fitness (e.g., larval duration, Quinn & Rochette, 2015), the effects on benthic phases are likely to be profound and continued.

Rates of benthic warming in wintertime exceeded summertime warming in many regions, particularly on Georges Banks, but also through the Gulf of Maine. A recent modeling effort found similar patterns, with fall warming rates exceeding that of spring and faster rates occurring on Georges Bank and the mid‐Atlantic (Kleisner et al., 2017). While the absolute rise in temperature may still allow for thermal refugia, the effect of warmer winter temperature on disease incidence and energetics should not be ignored. Warmer winters are associated with increased disease incidence and parasite loads in organisms due to increased development and transmission, reducing overwintering mortality, and reduced fitness of host (Harvell et al., 1999). Rising benthic temperatures have been linked to the spread and prevalence of epizootic shell disease, ESD (Maynard et al., 2016). ESD is caused by chitinolytic bacteria and affects American lobsters from Long Island Sound to Buzzards Bay, MA. Prevalence of ESD is rising, and is particularly high in gravid females (Maynard et al., 2016). The spread and survivorship of fungal and protozoan pathogens of Eastern Oyster have also been linked with warmer winter temperatures (Harvell et al., 1999). Cold wintertime temperatures are also critical for overwintering copepods. Increased temperatures are associated with reduced dormancy (Pierson et al., 2013), increased consumption of lipid stores (Runge et al., 2014) and reduced fecundity, negatively affecting trophic transfer to both pelagic and demersal fish. Organisms tend to integrate multiple stressors over a longer period of time (Pörtner, 2012) and species thermal adaptations affect their responses. For example, using biannual data from the NEFSC Trawl Survey, Kleisner et al. (2017) found that the range of suitable thermal envelopes increased for species with southerly distribution and decreased for species with more northerly distributions. However, while the seasonal surveys gather community structure and environmental data concurrently, the frequency of measurements may be insufficient to determine the duration or variability of a given thermal envelope (Kleisner et al., 2017).

Drivers of ecosystem variability do not operate in isolation (Nye et al., 2013; Gröger & Fogarty, 2011), and may have spatially distinct fingerprints in how they influence ecological communities (Bopp et al., 2013; Boyd et al., 2014). We evaluated the effects of interannual shifts in SST, winds, the AMO, NAO and PDO concurrently to understand the relative effect of local factors (e.g., winter mixing of warming surface waters to depth), regional circulation driven by changes in the AMO and the NAO, or atmospheric teleconnections (PDO).

SST had a strong positive influence on benthic temperature across the domain. Wind effects on benthic temperature were patchy but primarily positive in shallow regions of the mid‐Atlantic and southern New England shelves, and negative on the slope and deeper regions of the Gulf of Maine. Together these results suggest potential for mixing to be an important factor for increased benthic temperatures, particularly in the southern and central regions of the domain. Wind stress from 1982 to 2014 increased in the mid‐Atlantic and southern New England, but decreased in the northern Gulf of Maine, which is associated with a weak inverted thermal gradient in winter, with warmer temperatures on the bottom. Together with increased influx of warmer water associated with changes in circulation (discussed below), decreased winds may exacerbate winter heating as less heat is ventilated to the atmosphere.

Benthic temperatures were positively correlated with the AMO in the nearshore and Shelf waters, with the exception of the Chesapeake Bay. These patterns have important implications for ecosystems (Nye et al., 2013) including commercially important fisheries such as menhaden or cod (Buchheister et al., 2016). The AMO has been associated with coast‐wide recruitment patterns of menhaden, and has opposing effects in Chesapeake Bay (negative) compared to Southern New England (positive) (Wood & Austin, 2009), possibly due to the warmer winter temperatures required for spawning (>17°C; Light & Able, 2003). The AMO was associated with warming in the Gulf of Maine and subsequent declines in Atlantic cod (Pershing et al., 2015). Our study also found AMO‐associated benthic warming in the Gulf of Maine, but in addition noted strong effects of the NAO.

NAO influence is strongest in winter (Drinkwater et al., 2003). We found winter benthic temperatures have warmed over the past three decades, with the fastest rates occurring over Georges Bank, southern New England, but also in the Gulf of Maine. The NAO effect on benthic temperatures was strong (and negative) in the Gulf of Maine, strong (and positive) in the south Atlantic bight and neutral in the mid‐Atlantic, where modeling studies report only minor influence of circulation changes upstream on the Labrador Shelf (Chen et al., 2016). Ship‐based measurements of bottom temperatures over Georges Bank were four °C warmer than average during winter of 2011 (Gawarkiewicz et al., 2012). The temperature anomalies were associated with northward shifts in the distribution of Atlantic cod (Friedland et al., 2013) and increased abundances of warm water species reported by commercial fisherman off southern New England (Gawarkiewicz et al., 2014).

Benthic temperatures were negatively correlated with the PDO throughout the domain. This pattern is in agreement with earlier studies that found strong ties between the negative phase of the PDO and increased SST in the Gulf of Maine (Pershing et al., 2015). Negative PDO phases are associated with increased SST and air temperatures off of the U.S. southeast and mid‐Atlantic (Mantua et al., 1997). While this PDO phase is associated with increased upper atmospheric winds, surface winds have also increased in this region (Figure 8a). Thus, despite cooler temperatures predicted for the mid‐Atlantic from the negative NAO excursions (and increased inflow of LSSW), increased water and air temperatures associated with the PDO, decreased summer ventilation (Chen et al., 2014) and increased mixing of heat downward possibly superseded any cooling of benthic water in the mid‐Atlantic. In the Gulf of Maine, a negative PDO and negative NAO act synergistically to increase warming.

How the benthic thermal environment will change over the next century in response to future climate forcing is the focus of several studies (Kleisner et al., 2017; Rheuban et al., 2017; Saba et al., 2016). However, depending on the phase of AMO, an enhanced or reduced large‐scale warming could be observed in the future, and may interact with the projected anthropogenically induced warming. Closely associated with the AMO is the AMOC (Wang & Zhang, 2013), which is predicted to weaken by 25% over the next century (Meehl et al., 2007; Schmittner et al., 2005). For the Northwest Atlantic Shelf, this could result in a poleward shift of the north wall of the Gulf Stream and a reduction of transport of cooler water from the Scotian Shelf (Nye et al., 2013).

As anthropogenic warming continues, the interannual oscillations of the NAO will also play a large role in bottom temperatures and the dynamics and persistence of marine benthic species in the Gulf of Maine and throughout the mid‐Atlantic. How the frequency and magnitude of NAO will shift under climate projections is uncertain. In the early 1980s, the NAO shifted from several decades of predominantly negative to largely positive (Visbeck et al., 2001). This shift was attributed to warmer tropical SSTs and strengthened stratospheric vortex (Visbeck et al., 2001). The NAO has been mostly neutral from 2000 to 2010, which may have contributed to a weakening of the NAO and ecosystem state in the Gulf of Maine (Hare & Kane, 2012). Strong negative excursions in 2010 and 2012 that have coincided with negative PDO may have led to ocean heatwaves and the rapid warming in the Gulf of Maine occurring over the past ∼15 years (Persching et al., 2015; this study). Thus it is likely important to consider the interaction of the short‐term and long‐term phasing of the NAO with the PDO and AMO. Recent models have dramatically improved the skill of predicting NAO phasing and magnitude (Dunstone et al., 2016) by considering the short‐term and long‐term interactions of the NAO with sea‐ice and the strength of the Arctic vortex. Integrating these newer climate models into ecological forecasting tools will be an important step for managing benthic species in a changing climate.

We constructed a synoptic time series that merged predictions based on satellite remote sensing data and in situ observations. Synoptic, near real‐time satellite data are ideal to inform spatial management of dynamic species and regions (e.g., Kavanaugh et al., 2016; Lewison et al., 2015); however, satellite‐driven analyses are not without caveats. Satellite data only see the surface, thus we interpolated benthic temperatures using regionally specific statistical models of varying complexity. While most regions were well characterized, the model performed only moderately well in the mid‐Atlantic summer (Region 4), where strong vertical stratification may further decouple surface signals from bottom temperatures. Opaque clouds (if undetected) can introduce a negative bias in estimated surface temperatures, potentially resulting in underestimated rates of surface warming and predicted benthic warming, particularly if cloud cover increases in response to warming oceans. Previous investigations using the more cloud biased AVHRR Pathfinder version 5.2 data found no warming in the Gulf of Maine (Maynard et al., 2016), however their time series also ended prior to the 2011–2012 period. Nevertheless, satellite‐driven analyses (this study; Maynard et al., 2016) present comparable trends to modeling efforts (Saba et al., 2016), which predict increases in benthic temperatures of 0.4–0.6°C/decade under a doubling of preindustrial CO2.

Coupling ecological information with higher‐resolution dynamical and predictive ocean‐climate models will provide a powerful means of managing living marine resources in a changing climate (Kleisner et al., 2017; Rheuban et al., 2017; Tomassi et al., 2017). However, while based on known, mechanistic, physical relationships, many models fail to characterize seasonal and subseasonal variability due to issues and uncertainties representing subgrid cell process (e.g., mesoscale and submesoscale features), surface and lateral boundary forcing, and emergent properties of a local or regional ocean climate system. Forecasting skill also largely depends on initializing a model using information specific to the current conditions. Skillful statistical models, while empirical, may provide information on the spectrum of variability in addition to accurate assessment of the current conditions. Our satellite‐based model predicted seasonal and interannual variability of benthic temperatures with a relatively high degree of skill (R2 > 0.7). Future efforts may benefit by integrating multiplatform data from additional observational efforts including, for example, innovative collaborations of scientists and the commercial fishing community (e.g., Commercial Fisheries Research Foundation, http://science.whoi.edu/users/seasoar/cfrfwhoi/) and the Ocean Observing Initiative (OOI; http://oceanobservatories.org/) Coastal Pioneer array, where moorings and profiling gliders may contribute to both mechanistic understanding as well as fill in the gaps where model skill or resolution are lacking. Together, with regional and statistically downscaled models, these tools can provide an integrated approach for managers to understand and mitigate the effects of warming oceans on commercially and ecologically important benthic species.

Supporting Information S1

Click here for additional data file.