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

Effects and impacts of the Northeast Pacific marine heatwave of 2014–2016 on the inner coastal estuarine waters of the Salish Sea were examined using a combination of monitoring data and an established three-dimensional hydrodynamic and biogeochemical model of the region. The anomalous high temperatures reached the U.S. Pacific Northwest continental shelf toward the end of 2014 and primarily entered the Salish Sea waters through an existing strong estuarine exchange. Elevated temperatures up to + 2.3°C were observed at the monitoring stations throughout 2015 and 2016 relative to 2013 before dissipating in 2017. The hydrodynamic and biogeochemical responses to this circulating high-temperature event were examined using the Salish Sea Model over a 5-year window from 2013 to 2017. Responses of conventional water-quality indicator variables, such as temperature and salinity, nutrients and phytoplankton, zooplankton, dissolved oxygen, and pH, were evaluated relative to a baseline without the marine heatwave forcing. The simulation results relative to 2014 show 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). However, sensitivity tests show that this increase was a direct result of higher freshwater and associated nutrient loads accompanied by stronger estuarine exchange with the Pacific Ocean rather than warming due to the heatwave. Strong vertical circulation and mixing provided mitigation with only ≈+0.6°C domain-wide annual average temperature increase within Salish Sea, and served as a physical buffer to keep waters cooler relative to the continental shelf during the marine heatwave.


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.
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 higherthan-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 . 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)   (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 waterquality indicator variables and parameters of concern such as temperature (T) and salinity (S) [

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  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 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 4 , nitrite + nitrate (NO 2 + NO 3 ), labile and refractory dissolved organic nitrogen and particulate organic nitrogen, PO 4 , 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 2 , and pH . 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 2 + NO 3 , 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.
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.
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: where T = temperature ( • C), T opt = optimal temperature for plankton production ( • C), KTg 1 = effect of temperature below T opt on production ( • C −2 ), and KTg 2 = effect  of temperature above T opt on production ( • C −2 ). These parameters (KTg 1 , KTg 2 , 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(Keister et al., , 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.

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. 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 (downwardupward), (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.  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.

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 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 2 ) 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 5 m 3 /s (Sutherland et al., 2011;Khangaonkar et al., 2017Khangaonkar et al., , 2018Olson et al., 2020), MacCready et al., 2021). It is a function of mean water depth (∝ H 3 ) and depth-averaged density (salinity) gradient (∝ ∂S ∂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 subbasins. 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 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.
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.

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 Ahmed et al., 2019). 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 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.

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(Khangaonkar et al., , 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/Dataresources/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.