Propagation of the 2014–2016 Northeast Pacific Marine Heatwave Through the Salish Sea

Key Points

- A multi-year simulation of the Salish Sea successfully reproduced the warmer conditions, stronger circulation, and higher biological activity observed during Northeast Pacific marine heatwave.

– Sustained high temperatures may actually have increased stratification, lowered circulation, and reduced nutrients in the photic zone, and reduced diatom growth resulting in a slight lowering of algal growth relative to reference conditions.

– Sensitivity tests showed that increased biological activity observed during the heatwave was due to the dominant effect of stronger estuarine exchange, higher freshwater inflows, and higher nutrient loads, and counterintuitively not the warmer conditions as previously assumed.

Introduction

Anomalous high sea surface temperatures occupied a large region off the coast of North America during the winter of 2013. This mass of warm water, referred to as the “the blob,” was ≈3°C warmer than normal in February 2014 and expanded its spatial extent to reach the U.S. Pacific Northwest continental shelf toward the end of 2014 (Bond, 2014). The occurrence of this Northeast Pacific marine heatwave is attributed to decreased surface cooling and lower equatorward Ekman transport in the Gulf of Alaska due to a persistent atmospheric high-pressure ridge (Wang et al., 2014; Bond et al., 2015; Hartmann, 2015). The heatwave persisted over multiple years, from the end of 2014 to the end of 2016, before dissipating in 2017, but remains a concern for potential impacts to estuarine ecosystems and possible recurrences. Nearshore sea surface temperature anomalies along the Washington, Oregon, and California coasts reached a maximum of 6.2°C off Southern California and only abated seasonally during spring upwelling-favorable wind stress (Gentemann et al., 2017). This event has been linked to numerous incidences of complex ecological impacts over the continental shelf waters directly exposed to the heatwave (Kintisch, 2015). Whitney (2015) pointed out that just prior to the onset of the heatwave, the winter of 2013/2014 was accompanied by anomalous winds from the south that weakened nutrient transport in the eastern North Pacific, resulting in substantial decreases in phytoplankton biomass. This was followed by a coastwide bloom of the toxigenic diatom Pseudo-nitzschia in spring 2015 that resulted in the largest recorded outbreak of the neurotoxin, domoic acid, along the North American west coast. McCabe et al. (2016) demonstrated that this was caused by the anomalously warm ocean conditions. Similarly, Peterson et al. (2017) concluded that warm ocean conditions observed from late September 2014 through 2016 resulted in an overall drop in populations of copepods and potential collapse of the usual food chain up to forage fish (smalt, herring, and lance). Du and Peterson (2018) noted a shift in phytoplankton abundance and diversity from diatoms that favor lower temperatures to dinoflagellate that favor higher temperatures for optimal growth. More frequent and intense ocean warming events may have complex impacts on the food webs as concluded by Jones et al. (2018) based on an assessment of massive mortality of planktivorous seabirds following the heatwave.

As may be expected, the inner waters of coastal estuaries were also affected. This includes the Salish Sea, an inland fjord-like estuary comprised of the Puget Sound, the Strait of Juan de Fuca, and the Strait of Georgia basins in the U.S. and Canadian waters (see Figure 1). The Salish Sea is one of the world’s largest and most biologically rich inland seas that supports numerous species of mammals, birds, fish, and invertebrates, which in turn are vital to the regional economy, culture, and quality of life. There is continuing urgency to improve our understanding of the vulnerability or resilience of these coastal estuarine waterbodies to stressors such as anthropogenic nutrient and carbon pollution, sea-level rise, climate change, and now marine heatwaves.

FIGURE 1

Figure 1. Oceanographic regions of the Salish Sea (Puget Sound, the Strait of Juan de Fuca, and the Strait of Georgia). Inset shows the entire model domain. Monitoring data stations include Ecology monthly monitoring sites (green) and Puget Sound zooplankton monitoring sites (red). Transects A, B, and C are locations where volume flux (exchange flow) calculations were conducted.

Impacts on Salish Sea temperatures were seen throughout 2015 and 2016 before mostly dissipating in 2017. During this period, depth-averaged temperature increases up to 1.5°C were seen at the Washington State Department of Ecology monthly monitoring stations relative to 2013. Unlike typical conditions in such fjord-like estuaries, where variations are largest in surface waters, these high temperatures influenced by the higher-than-normal temperature of incoming exchange flow affected the entire water column. The largest annual average increase in surface water temperature of 1.89°C was recorded at the Admiralty Inlet station, while the largest annual average increase in bottom water temperature of 1.98°C was seen in Hood Canal. A daily average sea surface temperature increase as high as +2.32°C, was seen in 2015 at the National Oceanic and Atmospheric Administration (NOAA) National Data Buoy Center Seattle station. The presence of warmer and lower density water in the Strait of Juan de Fuca likely lowered the density gradient in this exchange pathway. Consequently, density-driven exchange flow would also have been reduced to a certain extent and may have affected overall circulation. Occurrences of early blooms of Alexandrium spp. in Hood Canal in April, of 2015 followed by an increase in 2016, and numerous examples of increased bacterial growth effects such as Vibrio-contaminated oysters, and impacts on higher trophic levels including herring, seabirds (e.g., rhinoceros auklets), and some species of marine mammals, may have been due to these warmer conditions (Bond, 2021). However, understanding the biogeochemical and ecological response of the Salish Sea to this marine heatwave thus far has remained qualitative and somewhat inconclusive in terms of quantifying the causes, extent, severity, and magnitude of projected impacts to the nearshore estuarine ecosystems. This could be partly because ecological monitoring and modeling assessments that were conducted during the marine heatwave years had a different focus, and the results simply did not show a particularly anomalous hydrodynamic or biogeochemical response inside the Salish Sea domain. These include modeling of years 2015–2017 by Olson et al. (2020) who successfully reproduced key aspects of the nitrogen cycle in the Strait of Georgia, region of the Salish Sea including the seasonal cycle and regional differences and highlighted mechanisms through which tidal currents lead to nutrient supply to the surface layer. They acknowledged that 2015 and 2016 were somewhat anomalous years in the Strait of Georgia due to the influence of the incursion of warm water associated with the heatwave. Similarly, Moore-Maley and Allen (2021) analyzed a 5-year hindcast record (2015–2019) along with model-data validation to examine wind-driven upwelling effects on surface nutrient delivery in the Strait of Georgia. MacCready et al. (2021) conducted a 3-year hindcast (2017–2019) covering mostly a period post marine heatwave, focused on circulation mixing and residence times in the Salish Sea. A marine heatwave impact or signature was hard to detect from the data and analysis presented in the above assessments as the warmer water had already occupied the Salish Sea or that marine heatwave had already passed by the time simulations were initiated. But a study by Jarníková et al. (2021) to characterize physical drivers of productivity dynamics did span years 2013–2016, and noticed potential inhibition of water column mixing by higher thermal stratification of the system, in 2015 relative to prior years which they attributed to possible marine heatwave effects.

This potential for marine heatwave-related impacts is particularly important to address given region-wide concerns associated with climate change and the general understanding that an increase in extreme event frequency, including marine heatwaves, may be expected (Pearce and Feng, 2013; Bond et al., 2015; Chen et al., 2015; Di Lorenzo and Mantua, 2016; Hobday et al., 2016; Scannell et al., 2016). Although climate change is not considered to be a major cause of the heatwave, findings by Laufkötter et al. (2020) validate the concern that the frequency of marine heatwaves has increased nearly twentyfold due to anthropogenic climate change effects. Marine heatwaves, which typically occurred once in hundreds to thousands of years in preindustrial times, are projected to occur on a decadal to annual basis if the global average air temperature rises by 1.5–3.0°C. The elevated temperatures during this marine heatwave serve as a preview of the type of conditions and estuarine biogeochemical stress we may expect under future climate conditions. The sustained high temperatures experienced here during the marine heatwave are comparable to the projected average high temperatures (+2.53°C on the shelf and + 1.5°C in the Salish Sea; Khangaonkar et al., 2019) corresponding to the high-emissions future scenario (the Intergovernmental Panel on Climate Change’s representative concentration pathway scenario, RCP8.5). The conditions experienced provide an opportunity for indirect validation of projected future ecological impacts.

In this paper, supported by synoptic field data, we present an assessment of the effects of the marine heatwave on the Salish Sea as a whole, including impacts on density-driven fjord-like estuarine circulation, followed by conventional water-quality indicator variables and parameters of concern such as temperature (T) and salinity (S) [nearshore habitat], dissolved oxygen (DO) [hypoxia], pH [ocean acidification], and nutrients and phytoplankton biomass [eutrophication], and zooplankton biomass [food web].

Materials and Methods

We used a modeling-based approach to isolate and distinguish the effects of the high-temperature pulse that propagated through the Salish Sea. A biophysical model of the Salish Sea was first validated over the 5-year marine heatwave window from 2013 to 2017. “Existing” conditions were compared with “reference” conditions—defined as existing conditions without the anomalous heatwave forcing. This companion “reference” simulation was created using climatological ocean boundary and average meteorological heat load from pre and post marine heatwave years.

The Salish Sea Model—Marine Heatwave Years 2013–2015

The Salish Sea Model (Khangaonkar et al., 2018) is a comprehensive biophysical model that uses the Finite Volume Community Ocean Model (FVCOM, Chen et al., 2003) framework to solve Reynolds-averaged Navier-Stokes equations for turbulent flows in coastal ocean environments. The model domain encompasses Vancouver Island completely, allowing tides to propagate into the Salish Sea around Vancouver Island through Johnstone Strait and the Strait of Juan de Fuca (Figure 1). The finite volume model grid consists of 16,012 nodes and 25,019 triangular elements. The vertical configuration of the model uses 10 sigma-stretched layers distributed using a power law function with an exponent P-Sigma of 1.5, which provides more layer density near the surface. The lateral resolution varies from cell size of 250 m near river mouths to 800 m in Puget Sound increasing to 3 km resolution in the straits and 12 km over the continental shelf. The hydrodynamic solutions are forced using ocean boundary tides and stratification (T, S) provided by the Hybrid Coordinate Ocean Model (HYCOM) and meteorological conditions from the University of Washington (UW) Weather Research Forecast model data. Tidal forcing at the open boundary is based on tidal constituents from the Eastern North Pacific (ENPAC) tidal database (Szpilka et al., 2018). The model includes a total of 99 wastewater discharges and stream flows from 161 watersheds developed by Ecology through a combination of monitoring data and multi-variate regression analysis (Mohamedali et al., 2011; Ahmed et al., 2019).

The biogeochemical component of the Salish Sea Model uses FVCOM-ICM developed originally by Kim and Khangaonkar (2012). It is based on the kinetics of CE-QUAL-ICM (Cerco and Cole, 1994, 1995) water-quality model that was further developed for use with FVCOM framework through external coupling. The baseline biogeochemical model includes 23 state variables: temperature, salinity, two species of phytoplankton (diatoms P1 and dinoflagellates P2), labile and refractory dissolved organic carbon (DOC) and particulate organic carbon (POC), NH₄, nitrite + nitrate (NO₂ + NO₃), labile and refractory dissolved organic nitrogen and particulate organic nitrogen, PO₄, labile and refractory dissolved organic phosphorous, DO, dissolved inorganic carbon (DIC), and total alkalinity (TA). The model was updated to include a sediment diagenesis module that allows directly coupled interaction between the water column and sediments through the processes of organic sediment settling, burial, remineralization, and carbonate chemistry with DIC, TA, pCO₂, and pH (Bianucci et al., 2018). New constituents such as inorganic suspended solids, turbidity, zooplankton, and submerged aquatic vegetation were recently added (Khangaonkar et al., 2021). Ocean boundary values for most water-quality variables were set based on available climatological data from a combination of world ocean atlas and the Canadian Department of Fisheries and Oceans monitoring database except for DO, NO₂ + NO₃, TA, and DIC variables that were specified using regressions to salinity developed by the UW (personal communication from Ryan McCabe and Parker MacCready of UW).

A 5-year long simulation was set up spanning the marine heatwave window from 2013 to 2017. Model performance was validated through comparison to monthly monitoring data from 23 water-quality stations collected by Ecology and bi-weekly zooplankton data collected from 16 sites beginning in 2014 provided by the University of Washington (Figure 1). A notable model upgrade made as part of this assessment was the change in bed friction from uniform roughness factor (z_o) to spatially varying with localized higher values near channel entrances and islands. This resulted in improved prediction of water surface elevation with relative root mean square error (RMSE) of ≈< 6% (0.27 m). This change to the hydrodynamic circulation and the need to ensure overall acceptable model performance required minor fine-tuning of some of the biogeochemical parameters. Specifically, maximum photosynthetic rates (PM1 and PM2) for diatom and dinoflagellates were increased from 250 to 300 g C g^–1 Chl d^–1 to provide a better match to observed primary production and ensure sufficient prey was available to match observed zooplankton biomass. The maximum prey consumption ration of microzooplankton was reduced from 0.83 to 0.55 g prey C g^–1 zooplankton C d^–1. Settling rates of particulate organic matter were increased from 5.0 to 7.5 m/d to ensure a desired match with observed DO levels in all years. Tables 1 and 2 list the major algal and zooplankton kinetic parameters including those updated through this effort.

TABLE 1

Table 1. A list of major algal kinetics (diatoms and dinoflagellates) model parameters as part of the FVCOM-ICM model implementation for the Salish Sea Model.

TABLE 2

Table 2. A list of major zooplankton parameters (microzooplankton and mesozooplankton as part of the FVCOM-ICM model implementation for the Salish Sea Model.

It is important to note the prior calibration effort attempts at explicit zooplankton simulation was based on a single year of data from 2014. The update here reflects a compromise that allows satisfactory model performance compliant with targets established by the local user community (RMSE_T < 1°C, RMSE_S < 1 psu, and RMSE_DO < 1 mg/L) over the 5 years of continuous simulation. Table 3 provides skill assessment results (error statistics and skill scores) for key model constituents, tabulated by each year, prepared as part of the model validation demonstration. Figure 2 provides an example of multi-year time series comparison between model predictions and measured data at a representative station (Elliot Bay) from the Puget Sound Basin region of the Salish Sea.

TABLE 3

Table 3. The Salish Sea Model overall error statistics and skill score for major constituents with calibration data spanning the Northeast Pacific marine heatwave period from 2013 to 2017.

FIGURE 2

Figure 2. Salish Sea Model validation: a time series comparison of model predictions and observed data for T, nitrate (NO₃ + NO₂), phytoplankton biomass as chlorophyll a, zooplankton biomass, DO, and pH from the surface (outflow) and bottom (inflow) layers from a representative station (Elliot Bay) in the Central Puget Sound region of the Salish Sea from 2013 through 2017.

Typically, the algal bloom occurs around spring (late March) and reaches peak biomass in summer, followed by a decline in the fall, noticeable particularly in sub-basins such as Hood Canal, South Puget Sound, Bellingham Bay, and Skagit Bay. The recently added zooplankton module improves the ecosystem prediction capability by facilitating linkage of primary production externally to the food web models which often utilize zooplankton as the lowest trophic level at the start of their simulations. The zooplankton module in FVCOM-ICM simulates biomass variability of two aggregated zooplankton groups (i.e., micro and mesozooplankton) that consume produced phytoplankton and particulate organic matter. The effect of temperature on the growth rates of phytoplankton as well as zooplankton is incorporated through a factor f(T) that is defined by the equations that have the following form:

$f\left(T\right)=e^{-KTg1{(T-T_{opt})}^2}\mathit{\hspace{1em}}\mathrm{when}T\le T_{opt}$

$f\left(T\right)=e^{-KTg2{(T_{opt}-T)}^2}\mathit{\hspace{1em}}\mathrm{when}T>T_{opt}$

where T = temperature (°C), T_opt = optimal temperature for plankton production (°C), KTg₁ = effect of temperature below T_opt on production (°C ^–2), and KTg₂ = effect of temperature above T_opt on production (°C ^–2). These parameters (KTg₁, KTg₂, and T_opt) are specified independently for each phytoplankton species and each zooplankton species (microzooplankton, and mesozooplankton). Through the model calibration process, optimum temperatures for phytoplanktom production in the model were set at 14 and 20°C for diatoms and dinoflagellates respectively. Optimum temperatures for zooplankton production in the model were set at 20°C for micro as well as meso zooplanktom.

Figure 3 shows a consolidated 4-year time series of mesozooplankton data from University of Washington from Puget Sound and San Juan Islands vicinity over the period 2014 to 2017 (Keister et al., 2017, 2019). The corresponding prediction of mesozooplankton biomass concentrations is plotted over the same period spanning the marine heatwave. The figure shows that the intra and interannual variability of mesozooplankton biomass in Puget Sound is reproduced reasonably well. The model predictions show a significant increase in zooplankton biomass for the year 2014 following 2013. Zooplankton data as well as model results show continued increase from 2015 to 2017 including the peak heatwave years. The shaded region in plot represents the 95th and 5th percentile range of predicted values that encompass the variation in observed zooplankton among various stations. Measurements of mesozooplankton show ≈21% average increase in biomass concentration over the period 2015–2017 compared to 2014 at the 16 measurement stations in Puget Sound that is reproduced reasonably well in the model results at an average increase of ≈20% at the same locations.

FIGURE 3

Figure 3. Predicted time series of mesozooplankton biomass from 2014 to 2017 along with measured data from 16 monitoring stations in the Salish Sea. The dark blue line is the mean of predicted zooplankton values from model nodes corresponding to the station locations. The shaded region represents the 95th and 5th percentile range of predicted zooplankton at the station locations at each time step. A specific color of a data point is associated with a unique monitoring station.

Simulation of Reference Condition—Years 2013–2015 Minus Heatwave Forcing

The challenge of assessing impacts from extreme events such as the marine heatwave is in distinguishing anomalous heating from natural interannual and seasonal variability. Warming of estuarine waters is controlled by three major sources of heat: (a) heat flux from the exchange with open ocean waters; (b) atmospheric heat flux through the sea surface; and (c) heat flux from river and wastewater inflows. In coastal ocean modeling, in the absence of data, it is a common practice to use climatological average ocean boundary conditions based on historical monitoring data in place of real-time observations or global model-derived data. This practice is considered reasonable as interannual variations in ocean chemistry at deep ocean boundary waters are relatively small and seasonal variations in temperature and salinity profiles follow a consistent pattern year-to-year except during extreme events that are recorded as anomalies. This has also been the case with this Northeast Pacific marine heatwave of 2014–2016, where warming of ≈3°C over the continental shelf is often described as an anomalous increase relative to 1981–2010 climatological normal conditions. In the case of the Salish Sea Model, the open ocean boundary condition is derived from global HYCOM prediction interpolated to the continental shelf ocean boundary nodes. For this assessment, a climatological reference ocean state was derived from 10-year HYCOM model predictions averaged monthly over a period from 1999 to 2009. Figure 4 shows a vertical transect along the Salish Sea Model boundary, starting from the southernmost node near Waldport, Oregon, U.S. to the northern end at Calvert Island in British Columbia, Canada. Heatwave related perturbation at the ocean boundary is illustrated through temperature differences relative to reference conditions in March 2014 prior to the heatwave and during the peak impingement in March 2015 as an example.

FIGURE 4

Figure 4. Temperature difference contour plot along the Salish Sea Model boundary at the continental shelf edge to highlight the incoming marine heatwave through the ocean boundary. Distance is measured along the ocean boundary from the southernmost boundary node near Waldport, Oregon, to the northern most at Calvert Island in British Columbia, Canada at the north of Vancouver Island. (A) Temperature difference along the ocean boundary in March 2014 prior to the marine heatwave impingement and (B) during peak marine heatwave impingement in March 2015. The difference is relative to climatological monthly average from 1999 to 2009 selected as the reference condition at the boundary.

The oceanic heat flux was accompanied and somewhat preceded by a sharp increase in net heat flux into the Salish Sea from atmospheric loading at the sea surface. Net surface heat flux is a summation of (a) shortwave solar radiation (downward-upward), (b) longwave solar radiation (downward-upward), (c) latent heat flux, and (d) sensible heat flux. In addition to meteorological conditions including air temperature and relative humidity, the overall magnitude is also dependent on sea surface temperature. The Salish Sea Model uses net heat flux input from the UW Weather Research and Forecasting (WRF) model (12 and 4 km versions) with albedo factors adjusted as part of temperature calibration. Figure 5 shows net heat flux values for the respective years. To isolate and remove the marine heatwave-induced increase in atmospheric heat flux, we focus on the change relative to ambient conditions immediately prior to and after the heatwave. An average annual net flux value of 45.2 W/m² was computed using WRF net heat flux data from 2012–2013 and 2017–2018 to generate the reference value. Average atmospheric net heat flux during the marine heatwave years (2014–2016) was higher (+19%) relative to the reference value of 45.2 W/m². The time series of hourly net heat flux values in each year were then scaled to the reference value. This provided an annual average net heat flux value of 45.2 W/m² for each year of the 5-year period spanning the heatwave from 2013 to 2017 for use in the reference condition simulation while retaining scaled hourly temporal characteristics.

FIGURE 5

Figure 5. Net heat flux input to the Salish Sea Model based on the UW WRF—12 km model with albedo factor calibration to Salish Sea temperature profile monthly monitoring data from 2013 to 2017. The reference level net heat flux of 45 W/m² is the average of pre- and post-marine heatwave years 2013 and 2017.

Hydrological heat flux from river and wastewater flow is controlled by total inflow and inflow temperature. A summary of total annual average freshwater flow rates into the inner Salish Sea domain is shown in Figure 6. These do not include flows to the model domain along the Washington and Vancouver Island Pacific Ocean coastlines or the flows into the northern portion of Johnstone Strait that is considered outside of the Salish Sea basin. The freshwater flows during the period of interest 2013–2017 are unique in that Year 2013 was a particularly low-flow year followed by a continuous stretch of average-to-high flow years in 2014, 2015, 2016, and 2017. The high flow magnitudes are comparable to prior high flow years. Based on determination by Bond et al. (2015) that the effect of warmer temperatures on precipitation was negligible, we adopted the assumption that river flows during the marine heatwave were not altered significantly due its effects. Rather than eliminating the interannual variability associated with local weather, the flows and river temperatures for reference conditions were left unchanged.

FIGURE 6

Figure 6. Summary of annual average inflows into the inner Salish Sea covering Georgia Basin, Puget Sound, and the northwest straits surrounding San Juan Islands and the Strait of Juan de Fuca. The shaded region represents the study period from 2013 to 2017, and the peak marine heatwave years 2015 and 2016 are indicated.

Results

Propagation of the marine heatwave into the Salish Sea/Puget Sound region was examined using depth-time profile contour plots of temperature difference (ΔT) at selected stations starting from (1) the ocean boundary location directly across from Neah Bay at the mouth of the Strait of Juan de Fuca, to (2) Neah Bay, through (3) Port Townsend at Admiralty Inlet, to (4) Elliot Bay in the main basin of Puget Sound, and on to (5) Gordon Point at the southern end of Puget Sound as shown in Figure 7. The maximum ΔT (≈3–4°C) occurred during the first impingement of the marine heatwave relative to reference conditions was noted at the sea surface and over the continental shelf boundary. The ΔT reduced mostly to < 2°C as the heatwave propagated into the Salish Sea but increased in South Puget Sound in the range of 2–3°C. Results show that the marine heatwave peaked over the continental shelf and simultaneously reached the Salish Sea entrance at the mouth of the Strait of Juan de Fuca in mid-May 2014. Within a span of about 3 months by mid-September 2014, the heatwave had propagated and dispersed over the entire Salish Sea, reaching the innermost regions such as South Puget Sound.

FIGURE 7

Figure 7. Depth-time profile contour plots of temperature difference (ΔT) at selected stations relative to reference condition normalized relative to the depth of the water colum (z/D): (1) ocean boundary directly across from Neah Bay at the mouth of Strait of Juan de Fuca, (2) Neah Bay, (3) Port Townsend, (4) Elliot Bay, and (5) Gordon Point.

Figure 7 panel (1) shows the presence of the warmer sea surface temperature over the continental shelf from 2014 through 2016. The effects of downwelling-induced mixing are noticeable, particularly during the winter months of 2015 and 2016, where the temperature anomaly is seen to extend throughout the water column. Figure 7 panels (2) and (3) show the warmer waters in the Strait of Juan de Fuca and Admiralty Inlet (Port Townsend), mixed over the entire water column, being transported into the Puget Sound region of Salish Sea [Figure 7 panel (4) Elliot Bay] and South Puget Sound [Figure 7 panel (5) Gordon Point]. This early heat transfer into the Salish Sea is attributed primarily to estuarine exchange flow, given minimal atmospheric heat flux during the winter months. The effect of the marine heatwave in the Salish Sea was compounded by atmospheric heating during the summer. Figure 7 panel (1) shows that by mid-2016, the marine heatwave had already dissipated from the continental shelf. However, the resident heat in the Salish Sea remained trapped until the end of 2016, before being flushed out and replaced by cooler waters in 2017. It is important to note that in the deeper waters and the innermost regions of the Salish Sea, the effects of the heatwave lingered almost up to 2018 (Psemp Marine Waters Workgroup, 2020). These results are consistent with findings by Jackson et al. (2018) that warm conditions (0.3–0.6°C warmer than the normal monthly average) persisted below the surface mixed layer through at least March 2018 in adjacent Canadian waters to the north.

Sustained warming of Salish Sea waters was due to the combined effect of higher atmospheric heat flux (+ 19% relative to reference 45.2 W/m²) and conveyance of warmer waters from the shelf via strong estuarine exchange flow into the Salish Sea. The magnitude of the estuarine flow is estimated in the range of 100–150 × 10⁵ m³/s (Sutherland et al., 2011; Khangaonkar et al., 2017, 2018; Olson et al., 2020), MacCready et al., 2021). It is a function of mean water depth (∝H³) and depth-averaged density (salinity) gradient ($\propto \frac{\partial \overline{S}}{\partial X}$) as presented analytically by Hansen and Rattray (1965), Dyer (1973), and MacCready (2004), MacCready (2007) for partially mixed estuaries, by Rattray (1967) for fjords, and by Khangaonkar et al. (2011) for Salish Sea sub-basins. The possibility that warmer waters occupying the Strait of Juan de Fuca may have adversely impacted density gradients and therefore the strength of exchange flow was examined.

The definition of exchange flow used here is “tidally averaged volume flux” across the selected transect that enters the estuary below the depth of net zero motion. Selected transects include (A) Strait of Juan de Fuca (inflow to the Salish Sea), (B) Haro Strait (inflow to Georgia Basin), and (C) Admiralty Inlet (inflow to Puget Sound). Table 4 presents exchange flow magnitudes computed at these transects for existing conditions from 2013 to 2017 including the marine heatwave years. Also presented are exchange flow magnitudes for reference for the same period along with percent difference. Results as expected show that during the marine heatwave, likely due to the reduced density gradient from higher temperatures in the Strait of Juan de Fuca channel, the strength of the exchange flow was reduced relative to reference conditions. This reduction was small (<−1% change) and seen in the Strait of Juan de Fuca and Haro Strait transects. Exchange flow into Puget Sound showed a small increase (< 1% change). This leads to the conclusion that the circulation and exchange flows within the Salish Sea were not significantly affected by the marine heatwave and remained strong, resulting in efficient warming of the Salish Sea by transport of warmer continental shelf waters. What is striking though is the ≈ 8% increase in magnitude of exchange flow that occurred during the period 2015–2017 relative to 2013. This increase in exchange flow is attributed entirely to increase in salinity/density gradient from the ≈10% higher average freshwater inflow in years 2014–2017 relative to 2013 and 2014. Note that high freshwater loads in 2014 occurred late in the year and as shown in Table 4, did not influence 2014 exchange flows significantly but contributed to the 2015 higher values. By the same argument, Year 2013 exchange flows should have benefitted from the high flows in 2012. But a cold start of the simulation was performed using 2013 flows and as a result 2012 flow effects are not included and 2013 exchange flows are likely under predicted.

TABLE 4

Table 4. Exchange flow magnitudes computed at transects (1) Strait of Juan de Fuca, (2) Haro Strait, and (3) Puget Sound for existing conditions from 2013 to 2017 including marine heatwave years and a comparison to reference conditions for the same period.

Figure 8 shows plan view contour plots of depth-averaged temperature difference ΔT, for the months of January, April, July, and October, representing winter, spring, summer, and fall seasons. Temperature response within the Salish Sea is complex as temperatures are affected by river flows and strong circulation effects modulated by reflux and mixing between numerous sills. Those effects are included in existing and reference simulations and likely cancel out as the effects of the heatwave on circulation were shown to be small (≈<1% relative change in magnitude). Therefore, the difference plots isolate direct effects of warming during the heatwave and the effect is seen clearly from fall of 2014 through the fall of 2016. The year 2013 temperatures prior to the marine heatwave were a bit cooler relative to reference, and year 2014 ΔT plots show warmer temperatures from April of 2014. In the summer and fall of 2014, as the marine heatwave propagated into the Salish Sea, inner sub-basins such as South Puget Sound, Lynch Cove, Sinclair Inlet, and Bellingham Bay and selected nearshore regions of the Strait of Georgia showed more pronounced warming effects. This is likely due to a combination higher net heat flux during the heatwave distributed over shallower depths in selected sub-basins. This anomalous local heating is also noticeable in Figure 7 Gordon Point station located in South Puget Sound. The noticeable warm waters over the continental shelf in January of 2015 and 2016 highlight the persistence of anomalous warmer waters at the Strait of Juan de Fuca entrance to the Salish Sea during winter months.

FIGURE 8

Figure 8. Temperature difference (ΔT) by season relative to reference conditions (without marine heatwave effects), depth-averaged contours with a focus on Salish Sea region for years 2013 to 2017.

The impact of the marine heatwave on water quality/biogeochemistry due to sustained higher temperatures in the entire Salish Sea during the marine heatwave are of interest. Given the range of seasonal variation in temperature and concentrations of water-quality constituents, the incremental effect of the marine heatwave is not immediately noticeable in the time series plot of temperature (Figure 2 top panel) from the Elliot Bay region and likely varied in magnitude site to site. Of most interest and concern is the potential exacerbation of algal blooms, hypoxia, and ocean acidification due to the marine heatwave, although major toxic blooms were largely confined to the open coasts. Figure 9 presents the difference between modeled existing and reference conditions in annual average levels of major biogeochemical constituents such as algal biomass, zooplankton biomass, DO, and pH presented in the form of depth-averaged, annual average contour plots. The top row provides an annual summary of the difference in elevated temperatures relative to the reference, showing higher average temperatures in 2015 and 2016. The second and third rows show the difference in algal and zooplankton biomass relative to the reference condition without the temperature increase. Contrary to expectations based on the observation data, despite sustained higher temperatures during the heatwave, the modeled effect on algal biomass and zooplankton biomass relative to reference conditions was small. Some of the shallow embayments and nearshore regions in Puget Sound and around Discovery Islands in the north of Georgia Basin, show an increase in algae and zooplankton biomass. But in most of the larger basins, the predicted change does not show a dramatic increase in biomass relative to reference conditions, but contrary to expectations, shows a small reduction (gradients toward blue shade). The effect of warmer waters on DO as seen in the model was therefore primarily due to lower saturation at higher temperatures from 2014 onward. Most Salish Sea regions show a drop in DO relative to reference conditions except selected locations of noticeably improved DO levels such as the Strait of Juan de Fuca, and surprisingly, Hood Canal, a region known for its annual formation of the hypoxic volume of water. Observations and model results show noticeably lower DO levels in 2015 and 2016 relative to 2013 and 2014 (interannual variation see Figure 10). But small improvement is noted during the marine heatwave relative to reference for the same years. The simulation results for pH also show a similar change relative to the reference conditions in most of the Salish Sea. Many of these regions showed a drop in pH from 2014 to 2016 but small improvement in pH levels is noted in the Strait of Juan de Fuca and Hood Canal relative to reference for the same years.

FIGURE 9

Figure 9. Model biogeochemical response difference by year (algae and zooplankton biomass, DO, and pH) relative to reference conditions (without marine heatwave effects); depth-averaged contours with a focus on the Salish Sea. Years 2015 and 2016 represent the peak marine heatwave period.

FIGURE 10

Figure 10. Time series of (A) algal (phytoplankton) biomass and (B) zooplankton biomass from the south Salish Sea region including Puget Sound, Strait of Juan de Fuca, and Georgia Strait south of the U.S.-Canada border. The blue line represents model results for the existing condition including the marine heatwave warming. The red line represents reference conditions without marine heatwave-related warming.

Discussion

Although the difference plots shown indicate that changes in biological activity and biomass due to the heatwave warming effects alone were relatively small, data and model results clearly show that biomass increased in 2015, 2016, and 2017. As summarized by Ummenhofer and Meehl (2017), the reported impacts of this Pacific Northwest heatwave have included the largest algal bloom on record that negatively impacted shellfish along the western coast of North America (NOAA, 2016), indicating increased biological response in the ecosystem. Available years of monthly monitoring of algal biomass by Ecology from 21 of 23 stations in the Salish Sea that was used in model validation also showed increased biomass: average measured algal biomass concentrations were higher in 2014, 2015, and 2016 by 4, 14, and 19% respectively relative to pre-marine heatwave year 2013. The data show that in 2017, measured concentration of phytoplankton biomass dropped back down to 2013 levels. Similarly, zooplankton monitoring data shown in Figure 3 from Puget Sound available from 2014 onward showed that average biomass concentrations in 2015, 2016, and 2017 were higher by nearly 30, 10, and 21% relative to 2014, matching an increased abundance of food sources associated with phytoplankton and particulate organic carbon. A review of monthly monitoring DO data from 20 Ecology stations in the Salish Sea showed that depth averaged concentrations increased in 2014 but were on average lower during the marine heatwave years 2015, 2016 by 2–3% relative to 2013 before recovering in 2017. The DO response also qualitatively reflects deterioration due to a combination of increased algal biomass in 2015 and 2016 and lower saturation due to higher temperatures before recovering in 2017.

The concurrence of (1) warmer waters from the heatwave and higher atmospheric net heat flux and (2) evidence of increased phytoplankton and zooplankton biomass from monitoring programs in the Salish Sea led to the expectation/speculation that the two may be directly related. However, this expectation precludes the possibility that interannual variability in freshwater inflows, associated nutrient loads, and magnitude of exchange with the Pacific Ocean can also result in a similar response and may be a dominant contributor resulting in this observed increase in biological response. Past work in the Salish Sea in connection with nutrient pollution management has shown that the system is extremely sensitive to nutrient loads from rivers, wastewater discharges, and oceanic loads through exchange flow (Khangaonkar et al., 2018; Ahmed et al., 2019). Figure 10 shows simulated algal biomass in the photic zone (25 m) from the simulations conducted from 2013 through 2017. The annual algal growth and die-off pattern show similarities but also include subtle variations. For 2013, the total biomass (area under the curve) is larger than in 2014, but 2014 shows a higher peak. The peak drops off a bit in 2015 and 2016 before increasing in 2017. This pattern in Figure 10A appears to match the observed net freshwater inflow (and associated nutrient) loading to the system but is difficult to interpret in a stand-alone manner because a part of algal biomass is consumed by zooplankton. The transmitted effect is more visible in increasing zooplankton biomass in years 2014, 2015, and 2016 (Figure 10B).

The year 2014 was a higher-flow year immediately following 2013, which was a low-flow year. It is important to note that higher nutrient loads from 2014 also benefited 2015 algal growth as most of the 2014 high loads occurred during late fall and winter likely after the algal activity had ceased. The effect of sustained higher flow/nutrient loads in 2015, 2016, and 2017 is reflected in higher algal biomass in Figure 10A. This, in turn, results in higher zooplankton biomass, as shown in Figure 10B. The blue line in Figures 10A,B represents existing conditions with full marine heatwave effects included. The red line represents simulated reference conditions, cooler waters using climatological ocean boundary and reduced heating, and all other inputs remaining the same. Two features stand out in this result: (1) interannual variability and increased biological productivity have been reproduced in both simulations and appear to be dominated by nutrient loads and estuarine circulation/exchange and (2) heatwave effects in the Salish Sea appear to have caused a minor reduction in biological response from what may have otherwise occurred, in that the reference simulation biomass is a little higher (≈0–4% higher algal biomass in reference condition years 2014 to 2017) than the existing condition simulation with heatwave effects included. Relative to existing conditions, the zooplankton biomass in the reference simulation is also correspondingly lower (≈0–3% higher algal biomass in reference condition years 2014 to 2017).

Figure 11 presents time series plots of the volume of hypoxic water with DO < 2 mg/L, and volume of corrosive water with aragonite saturation state.

FIGURE 11

Figure 11. Time series of (A) volume of hypoxic water (DO < 2 mg/L) and (B) acidic water (Ω_A < 1) in the Salish Sea region including Puget Sound, Strait of Juan de Fuca, and Georgia Strait. The blue line represents model results for the existing conditions including the marine heatwave warming. The red line represents reference conditions without marine heatwave-related warming. The gray line represents the difference between existing and reference conditions showing the perturbation due to the heatwave.

${\text{Ω}}_A=\frac{[\left(C{a}^{2+}\right)]\cdot [\left(C{O}_3^{2-}\right)]}{[CaC{O}_3]}<1.$

Figure 11A shows the time series of hypoxic volume in the Salish Sea for the existing and reference conditions without heatwave effects. Results show that the volume of hypoxic water increases steadily over the years 2013, 2014, and 2015 and maintains its high level through 2016 and 2017, consistent with the effects anticipated from increased primary production over these years. The peak of hypoxia in the Salish Sea occurs during the fall months of September and October at the culmination of the algal growth period and does not appear to be significantly affected by the heatwave perturbation. The relative impact of the heatwave on hypoxic volume is plotted as the difference between existing and reference conditions. The difference relative to the reference, plotted as a gray curve, indicates that the exposure to hypoxic water (volume-days) is generally reduced during 2015, 2016, and 2017 compared to conditions that might have occurred without the heatwave. This is also consistent with improved DO in Hood Canal presented in Figure 9 relative to conditions without the heatwave warming.

Figure 11B shows the time series of the volume of acidic water in the Salish Sea for the existing and reference conditions without heatwave effects. As shown in Figure 11B, peak acidic volume occurs during the winter months due to lower pH associated with freshwater inflows. Unlike DO, the time series of the volume of water with Ω_A < 1 does not show a noticeable trend over the years 2013 to 2017. The predicted difference relative to reference, plotted as a gray curve, indicates that the exposure to corrosive water (volume-days) is likely reduced during the 2015, 2016, and 2017 heatwave compared to conditions without heatwave effects.

While the biogeochemical response with respect to heatwave perturbation was counter to our expectations, the behavior when examined in its entirety correlates quite well with the interannual variation in conventional key variables that the system is sensitive to. Figure 12 summarizes the Salish Sea conditions for key nutrient loads, exchange flow magnitudes, algal biomass, and zooplankton biomass for the years 2013 through 2017. Annual average nutrient loads, shown in Figure 12A, increased during 2014 to 2017 relative to pre-marine heatwave conditions of 2013, mirroring the increase in freshwater inflow (see Figure 6). Most of the 2014 freshwater inflow increase occurred late in the year, the effects of which carried over to 2015. The increased freshwater-induced salinity gradients resulted in higher exchange flows to the Salish Sea in 2015, 2016, and 2017, relative to 2013 and 2014 as shown in Figure 12B and Table 4. During this period, the Salish Sea experienced an increase in nutrient flux from land-based freshwater sources as well as from the Pacific Ocean. The higher loads and stronger tidal exchange resulted in the availability of more nutrients in the photic zone that led to higher primary productivity, which was transmitted to zooplankton through predation on algae, as shown in Figures 12C,D. Simulation results show a Salish Sea wide average increase in biological activity (+14%, and 6% Δ phytoplankton biomass, respectively) during the peak heatwave year 2015 and 2016 propagating toward higher zooplankton biomass (+14%, +18% Δ zooplankton biomass) relative to 2014. Biomass calculations are volume weighted average values computed over the entire Salish Sea domain for each year. Figures 12B–D show a comparison with simulated results without the effect of heatwave-related warming which once again demonstrate that the heatwave effect on algal growth was small relative to interannual variability in nutrient loads and estuarine processes.

FIGURE 12

Figure 12. Comparison of modeled annual average (A) nutrient loads (estimated using observed flow based regressions), (B) magnitudes of estuarine exchange flow, (C) phytoplankton biomass, and (D) zooplankton biomass. For (B–D), existing conditions simulation with the marine heatwave are compared to the reference condition simulations without the anomalous warming.

Model Limitations

All models have errors and limitations that arise from a combination of simplifying complex hydrodynamic and biogeochemical processes in the mathematical formulation, errors in the solution scheme discretization, lack of adequate site-specific data, and temporal and spatial resolution in model inputs and forcing parameters. Understanding model limitations is essential to ensure that application results are not misused or applied beyond their intended performance design and the deliverables presented are correctly interpreted. Some of the limitations are associated with assumptions and simplifications made as part of application to this study and are listed below.

• Definition of Reference Conditions. Distinguishing marine heatwave effects from normal/reference conditions has been done using an arbitrary definition of reference scenario. Use of a combination of climatological average for the ocean boundary along with a pre/post heatwave average of atmospheric heating is an approximation for normal conditions. This smooths out interannual variability of oceanic and atmospheric heat fluxes in the reference conditions and is a limitation.

• Assumption that hydrology and freshwater inflows were not affected by the marine heatwave. This assumption allows us to retain interannual variability in loading and keep inflows the same between existing and reference conditions. However potential error in this simplification could affect key conclusion in this study with respect to marine heatwave impacts on exchange flow as well as freshwater and nutrient loading.

• Accuracy of global operational models during extremes and marine heatwave conditions. The inherent assumption that operational global (HYCOM) and regional (WRF) models used to specify ocean and atmospheric forcing have accurately captured the anomalous effects of the heatwave conditions is a limitation.

• Lag in the initiation of zooplankton growth. A lag between observed initiation of zooplankton growth and model prediction is noticeable in Figures 2 and 3. There could be many reasons for why the model response lags observed initiation of growth including (a) specification of optimum temperature for growth, (b) prey availability timing (c) nutrient availability in the photic zone, and (d) variations in the response of sub-domains affected by inter basin exchanges etc.

• Fixed biogeochemical model coefficients during the heatwave. The biogeochemical model reaction rates and constants were unchanged for existing as well as reference simulation is a gross simplification. Although temperature dependence of various parameters and reaction rates is included in the model formulation, the potential change in biogeochemical behavior with sustained elevated temperatures for 2 years is unaccounted for. This includes complex response such as a shift toward ecological species that favor higher temperature with potential to drive significant ecological changes especially during the warmest times of the year. Similarly, potential replacement by diatom species adapted to warmer waters or by other phytoplankton types was not considered.

Conclusion

The results of the 5-year simulation spanning the heatwave validated against monitoring data from the same period showed clearly that there was a notable increase in biological activity during the years 2015–2017 relative to pre-heatwave conditions of 2013. The data and simulation results relative to 2014 showed an increase in biological activity (+ 14%, and + 6% Δ phytoplankton biomass, respectively) during the peak heatwave year 2015 and 2016 propagating toward higher zooplankton biomass (+14%, +18% Δ mesozooplankton biomass). This was accompanied by lower DO levels and an increase in the volume of hypoxic water. This behavior was consistent with expectation that sustained warmer waters in the Salish Sea would result in higher biological activity. However, sensitivity test simulations of reference conditions without the heatwave perturbation showed that the warming with all other conditions unchanged would potentially have caused an opposite effect on algal growth in the Salish Sea that is dominated by diatoms which favor cooler temperatures. The increase in biological activity experienced by the Salish Sea was most likely due to year-to-year variation in key parameters that Salish Sea’s biogeochemistry is most sensitive to: (a) nutrient loads from land-based sources and (b) magnitude of the estuarine exchange with the Pacific Ocean and associated oceanic supply of nutrients. Both quantities are directly influenced by freshwater inflow to the Salish Sea. Freshwater during the marine heatwave went through a sharp change from a low-flow condition in 2013 to a sustained period of higher than average flows from 2014 to 2017. Based on these results, it appears this increase in hydrologic loads was the primary driver of increased biological activity during the passage of the 2014–2016 Northeast Pacific marine heatwave through the Salish Sea. However, we also acknowledge that the question whether this sustained multi-year higher-than-normal freshwater inflow to the Salish Sea, was influenced by the marine heatwave due to factors other than precipitation such as higher snow melt, remains unresolved.

Relative to reference conditions each year, the sensitivity tests show a small reduction in primary production during the heatwave. This result is contrary to expectation of higher biological activity at higher temperatures and must be treated with caution. But, several factors have contributed to this result. Stronger stratification during the marine heatwave likely reduced the diffusive flux of nutrients to the photic zone. While some sub-basins experienced higher temperature increases (≈2.0°C), the average temperature increase in the Salish Sea was only ≈0.6°C thanks to the mitigating benefit provided by strong circulation between the sills present in this fjord like estuary. As expected, there was a small reduction in the growth of diatoms at higher marine heatwave temperatures, but a corresponding increase in higher temperature favoring dinoflagellates was not realized due to stronger control exerted by light availability later in the summer.

We note that most literature and reported impacts from the heatwave were based on studies from coastal regions directly exposed to the heatwave over the continental shelf. The results are consistent with our prior finding (Khangaonkar et al., 2019, 2021) that strong physical circulation and reflux mixing processes within the Salish Sea attenuates temperature impacts from the continental shelf.

Data Availability Statement

Publicly available datasets were analyzed in this study. This data can be found here: https://ecology.wa.gov/Research-Data/Data-resources/Environmental-Information-Management-database.

Author Contributions

All authors listed have made a substantial, direct, and intellectual contribution to the work, and approved it for publication.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We would like to express our appreciation for technical input and contributions from our colleagues and collaborators— Ben Cope and Angela Adams of U.S. Environmental Protection Agency who have supported the development of the Salish Sea Model through an Interagency Agreement with Pacific Northwest National Laboratory; Anise Ahmed and Cristiana Figueroa-Kaminsky of Washington State Department of Ecology who provided us with land-based inflows and nutrients data; and collaborators from King County and the University of Washington for providing us with zooplankton data without which this work would not have been possible. We would also like to acknowledge support from Dr. Joel Baker and the Puget Sound Institute for funding a portion of this research. The monitoring data used in the model setup and calibration and the data generated through model application and sensitivity tests are available upon request and will be posted on the Salish Sea Model website: http://salish-sea.pnnl.gov.

References

Ahmed, A., Figueroa-Kaminsky, C., Gala, J., Mohamedali, T., Pelletier, G., and McCarthy, S. (2019). Puget Sound Nutrient Source Reduction Project: Model Updates and Bounding Scenarios. Washington, DC: Department of Ecology.

Google Scholar

Bianucci, L., Long, W., Khangaonkar, T., Pelletier, G., Ahmed, A., Mohamedali, T., et al. (2018). Sensitivity of the regional ocean acidification and carbonate system in Puget Sound to ocean and freshwater inputs. Elem. Sci. Anth. 6:22. doi: 10.1525/elementa.151

CrossRef Full Text | Google Scholar

Bond, N. (2014). The Blob: Warm Water off the Coast of the PNW and What it May Mean for Our Summer Weather - A message from the State Climatologist. Available online at: https://climate.washington.edu/wp-content/uploads/2014Jun.pdf (accessed 22 September, 2021).

Google Scholar

Bond, N. (2021). “The blob: what were its impacts on the salish sea?,” in The State of the Salish Sea, eds K. L. Sobocinsky, G. Broadhurst, and N. J. K. Baloy (Bellingham, WA: Salish Sea Institute, Western Washington University). doi: 10.25710/vfhb-3a69

CrossRef Full Text | Google Scholar

Bond, N. A., Cronin, M. F., Freeland, H., and Mantua, N. (2015). Causes and impacts of the 2014 warm anomaly in the NE Pacific, Geophys. Res. Lett. 42, 3414–3420. doi: 10.1002/2015GL063306

CrossRef Full Text | Google Scholar

Cerco, C., and Cole, T. (1994). Three-dimensional Eutrophication Model of Chesapeake Bay. Washington, DC: U.S. Army Corps of Engineers, Vicksburg, Mississippi. Technical Report EL-94-4.

Google Scholar

Cerco, C. F., and Cole, T. M. (1995). User’s Guide to the CE-QUAL-ICM Three-Dimensional Eutrophication Model, Release Version 1.0. Washington, D.C: U.S. Army Corps of Engineers.

Google Scholar

Chen, C., Liu, H., and Beardsley, R. C. (2003). An unstructured, finite-volume, three-dimensional, primitive equation ocean model: application to coastal ocean and estuaries. J. Atmos. Ocean. Technol. 20, 159–186. doi: 10.1175/1520-0426(2003)020<0159:AUGFVT>2.0.CO;2

CrossRef Full Text | Google Scholar

Chen, K., Gawarkiewicz, G., Kwon, Y.-O., and Zhang, W. G. (2015). The role of atmospheric forcing versus ocean advection during the extreme warming of the Northeast U.S. continental shelf in 2012. J. Geophys. Res. Oceans 120, 4324–4339. doi: 10.1002/2014JC010547

CrossRef Full Text | Google Scholar

Di Lorenzo, E., and Mantua, N. (2016). Multi-year persistence of the 2014/15 North Pacific marine heatwave. Nat. Clim. Change 6, 1042–1047. doi: 10.1038/nclimate3082

CrossRef Full Text | Google Scholar

Du, X., and Peterson, W. T. (2018). Phytoplankton community structure in 2011–2013 compared to the extratropical warming event of 2014–2015. Geophys. Res. Lett. 45, 1534–1540. doi: 10.1002/2017GL076199

CrossRef Full Text | Google Scholar

Dyer, K. R. (1973). Estuaries: A Physical Introduction. New York, N.Y: John Wiley.

Google Scholar

Gentemann, C. L., Fewings, M. R., and García-entem, M. (2017). Satellite sea surface temperatures along the west coast of the United States during the 2014–2016 Northeast Pacific marine heatwave. Geophys. Res. Lett. 44, 312–319. doi: 10.1002/2016GL071039

CrossRef Full Text | Google Scholar

Hansen, D. V., and Rattray, M. (1965). Gravitational circulation in straits and estuaries. J. Mar. Res. 2, 104–122.

Google Scholar

Hartmann, D. L. (2015). Pacific sea surface temperature and the winter of 2014. Geophys. Res. Lett. 42, 1894–1902. doi: 10.1002/2015GL063083

CrossRef Full Text | Google Scholar

Hobday, A. J., Alexander, L. V., Perkins, S. E., Smale, D. A., Straub, S. C., Oliver, E. C. J., et al. (2016). A hierarchical approach to defining marine heatwaves. Prog. Oceanography 141, 227–238. doi: 10.1016/j.pocean.2015.12.014

CrossRef Full Text | Google Scholar

Jackson, J. M., Johnson, G. C., Dosser, H. V., and Ross, T. (2018). Warming from recent marine heatwave lingers in deep British Columbia fjord. Geophys. Res. Lett. 45, 9757–9764. doi: 10.1029/2018GL078971

CrossRef Full Text | Google Scholar

JarníkováT., Olson, E. M., Allen, S. E., Ianson, D., and Suchy, K. D. (2021). A clustering approach to determine biophysical provinces and physical drivers of productivity dynamics in a complex coastal sea. J. Ocean Sci. doi: 10.5194/os-2021-66

CrossRef Full Text | Google Scholar

Jones, T., Parrish, J. K., Peterson, W. T., Bjorkstedt, E. P., Bond, N., Ballance, L. T., et al. (2018). Massive mortality of a planktivorous seabird in response to a marine heatwave. Geophys. Res. Lett. 45, 3193–3202. doi: 10.1002/2017GL076164

CrossRef Full Text | Google Scholar

Keister, J., Winans, A., and Hermann, B. (2017). Salish Sea Marine Survival Project: Zooplankton Monitoring Program 20142015 Final Report. Washington, DC: University of Washington, School of Oceanography, Seattle.

Google Scholar

Keister, J. E., Herrmann, B., and Winans, A. K. (2019). Puget Sound 2014-2018 Zooplankton Monitoring Report. Washington, DC: Prepared for Long Live the Kings, Seattle.

Google Scholar

Khangaonkar, T., Long, W., and Xu, W. (2017). Assessment of circulation and inter-basin transport in the Salish Sea including Johnstone Strait and Discovery Islands pathways. Ocean Model. 109, 11–32. doi: 10.1016/j.ocemod.2016.11.004

CrossRef Full Text | Google Scholar

Khangaonkar, T., Nugraha, A., Xu, W., and Balaguru, K. (2019). Salish sea response to global climate change, sea level rise, and future nutrient loads. J. Geophys. Res. Oceans 124, 3876–3904. doi: 10.1029/2018JC014670

CrossRef Full Text | Google Scholar

Khangaonkar, T., Nugraha, A., Xu, W., Long, W., Bianucci, L., Ahmed, A., et al. (2018). Analysis of hypoxia and sensitivity to nutrient pollution in Salish Sea. J. Geophys. Res. Oceans 123, 4735–4761. doi: 10.1029/2017JC013650

CrossRef Full Text | Google Scholar

Khangaonkar, T., Yang, Z., Kim, T., and Roberts, M. (2011). Tidally averaged circulation in puget sound sub-basins: comparison of historical data, analytical model, and numerical model. Estuar. Coast. Shelf. S 93, 305–319. doi: 10.1016/j.ecss.2011.04.016

CrossRef Full Text | Google Scholar

Khangaonkar, T. K., Premathilake, L., Nugraha, A., Keister, J., and Borde, A. (2021). Projections of algae, eelgrass, and zooplankton ecological interactions in the inner Salish Sea – for future climate, and altered oceanic states. Ecol. Modell. 441:109420. doi: 10.1016/j.ecolmodel.2020.109420

CrossRef Full Text | Google Scholar

Kim, T., and Khangaonkar, T. (2012). An offline Unstructured Biogeochemical Model (UBM) for complex estuarine and coastal environments. Environ. Modell. Softw. 31, 47–63. doi: 10.1016/j.envsoft.2011.11.010

CrossRef Full Text | Google Scholar

Kintisch, E. (2015). ‘The Blob’ invades Pacific, flummoxing climate experts. Science 348, 17–18. doi: 10.1126/science.348.6230.17

PubMed Abstract | CrossRef Full Text | Google Scholar

Laufkötter, C., Zscheischler, J., and Frölicher, T. L. (2020). High-impact marine heatwaves attributable to human-induced global warming. Science 369, 1621–1625. doi: 10.1126/science.aba0690

CrossRef Full Text | PubMed Abstract | Google Scholar

MacCready, P. (2004). Toward a unified theory of tidally-averaged estuarine salinity structure. Estuaries 27, 561–570. doi: 10.1007/bf02907644

CrossRef Full Text | Google Scholar

MacCready, P. (2007). Estuarine adjustment. J. Phys. Oceanography 37, 2133–2145. doi: 10.1175/jpo3082.1

CrossRef Full Text | Google Scholar

MacCready, P., McCabe, R. M., Siedlecki, S., Lorenz, M., Giddings, S. N., Bos, J., et al. (2021). Estuarine circulation, mixing, and residence times in the Salish Sea. J. Geophys. Res. Oceans 126:e2020JC016738. doi: 10.1029/2020JC016738

CrossRef Full Text | Google Scholar

McCabe, R. M., Hickey, B. M., Kudela, R. M., Lefebvre, K. A., Adams, N. G., Bill, B. D., et al. (2016). An unprecedente d coastwide toxic algal bloom linked to anomalous ocean conditions. Geophys. Res. Lett. 43, 366–10,376. doi: 10.1002/2016GL070023.

PubMed Abstract | CrossRef Full Text | Google Scholar

Mohamedali, T., Roberts, M., Sackmann, B. S., and Kolosseus A. (2011). Puget Sound Dissolved Oxygen Model: Nutrient Load Summary for 1999–2008. Publication No. 11-03-057. Olympia: Washington State Department of Ecology.

Google Scholar

Moore-Maley, B., and Allen, S. E. (2021). Wind driven upwelling and surface nutrient delivery in a semi-enclosed coastal sea. J. Ocean Sci. doi: 10.5194/os-2021-21

CrossRef Full Text | Google Scholar

NOAA (2016). NOAA Fisheries Mobilizes to Gauge Unprecedented West Coast Toxic Algal Bloom. Washington, DC: NOAA.

Google Scholar

Olson, E. M., Allen, S. E., Do, V., Dunphy, M., and Ianson, D. (2020). Assessment of nutrient supply by a tidal jet in the northern Strait of Georgia based on a biogeochemical model. J. Geophys. Res. Oceans 125:e2019JC015766. doi: 10.1029/2019JC015766

CrossRef Full Text | Google Scholar

Pearce, A. F., and Feng, M. (2013). The rise and fall of the ‘ marine heatwave’ off Western Australia during the summer of 2010/2011. J. Mar. Syst. 111–112, 139–156. doi: 10.1016/j.jmarsys.2012.10.009

CrossRef Full Text | Google Scholar

Peterson, W. T., Fisher, J. L., Strub, P. T., Du, X., Risien, C., Peterson, J., et al. (2017). The pelagic ecosystem in the Northern California current off oregon during the 2014–2016 warm anomalies within the context of the past 20 years. J. Geophys. Res. Oceans 122, 7267–7290. doi: 10.1002/2017JC012952

PubMed Abstract | CrossRef Full Text | Google Scholar

Psemp Marine Waters Workgroup (2020). Puget Sound Marine Waters: 2019 Overview. California, CA: PSEMP.

Google Scholar

Rattray, M. Jr. (1967). “Some aspects of the dynamics of circulation in fjords,” in Estuaries, ed. G. H. Lauff (Washington D.C: America Association for the Advancement of Science Publication).

Google Scholar

Scannell, H., Pershing, A., Alexander, A. M., Thomas, A. C., and Mills, K. E. (2016). Frequency of marine heatwaves in the North Atlantic and North Pacific since 1950. Geophys. Res. Lett. 43, 2069–2076. doi: 10.1002/2015GL067308

CrossRef Full Text | Google Scholar

Sutherland, D. A., MacCready, P., Banas, N. S., and Smedstad, L. F. (2011). A model study of the Salish Sea estuarine circulation. J. Phys. Oceanog. 41, 125–1143. doi: 10.1175/2011JPO4540.1

CrossRef Full Text | Google Scholar

Szpilka, C., Dresback, K., Kolar, R., and Massey, T. C. (2018). Improvements for the Eastern North Pacific ADCIRC tidal database (ENPAC15). J. Mar. Sci. Eng. 2018:131. doi: 10.3390/jmse6040131

CrossRef Full Text | Google Scholar

Ummenhofer, C. C., and Meehl, G. A. (2017). Extreme weather and climate events with ecological relevance: a review. Philos. Trans. R. Soc. B Biol. Sci. 372:20160135.

Google Scholar

Wang, S. Y., Hipps, L., Gillies, R. R., and Yoon, J.-H. (2014). Probable causes of the abnormal ridge accompanying the 2013–2014 Californiadrought: ENSO precursor and anthropogenic warming footprint. Geophys. Res. Lett. 41, 3220–3226. doi: 10.1002/2014GL059748

CrossRef Full Text | Google Scholar

Whitney, F. A. (2015). Anomalous winter winds decrease 2014 transition zone productivity in the NE Pacific. Geophys. Res. Lett. 42, 428–431. doi: 10.1002/2014GL062634

CrossRef Full Text | Google Scholar

Willmott, C. J. (1982). Some comments on the evaluation of model performance. Bull. Amer. Meteor. Soc. 63, 1309–1313. doi: 10.1175/1520-0477(1982)063<1309:SCOTEO>2.0.CO;2

CrossRef Full Text | Google Scholar