Introduction

Marine heatwaves (MHWs) have become distinguishable features of northeast (NE) Pacific Ocean temperature variability that disrupt the productivity of marine ecosystems and their services (Smale et al., 2019). These prolonged, discrete, and anomalously warm water events (Hobday et al., 2016) are most recognizable at the sea surface and are influenced by anthropogenic warming (Laufkötter et al., 2020). The effects of long-term ocean warming have led to a near doubling in the average annual count of MHW days globally since the early twentieth century (Oliver et al., 2018). Although MHWs have occurred throughout the global ocean, the NE Pacific has recently emerged as a hot spot for extremely persistent and large-scale events that are forced by anomalous air-sea heat flux driven by remote forcing from the tropics (Di Lorenzo & Mantua, 2016;Holbrook et al., 2019), in addition to long-term warming from anthropogenic greenhouse forcing (Laufkötter et al., 2020). The most remarkable NE Pacific MHWs have occurred in 2013-2016 and 2019-2020 and are colloquially referred to as The Blob (Bond et al., 2015) and Blob2.0 (Amaya et al., 2020), respectively (Figure 1 and Figure S1 in the supporting information).

The magnitude of sea surface temperature (SST) anomalies associated with MHWs depends critically on the seasonal evolution of the mixed layer depth (MLD), which deepens in winter and shoals in summer. If winter mixed layer MHW anomalies are present in the early spring when the mixed layer shoals, they can become trapped in the subsurface during the summer through detrainment. These detrained temperature anomalies are then stored in the subsurface and can reemerge the following winter when the mixed layer deepens and reentrains them (Alexander & Deser, 1995;Alexander et al., 1999Alexander et al., , 2001)). Alternatively, in the presence of downward Ekman pumping from wind stress curl, for example, in the North Pacific subtropical gyre, detrained anomalies can subduct, where they are further isolated from the mixed layer (Qiu & Huang, 1995). Here, we explore the role of detrainment and subduction in the sequestration of MHW anomalies into the permanent pycnocline where they can persist for years.

The evolution of the 2013-2016 NE Pacific MHW was complex and shaped by multiple drivers. Warm SST anomalies first appeared in the southern Gulf of Alaska centered on 40°N and 150°W and subsequently developed along the coast and south into the Southern California Current System near 25°N. In the Gulf of Alaska, lower rates of turbulent heat loss during the winter of 2013-2014 from the ocean to atmosphere and a reduction in wind-generated stirring allowed the winter mixed layer to remain unseasonably warm and shallow (Bond et al., 2015). The MHW transitioned to the south owing to local positive downward shortwave radiation anomalies and a positive SST-cloud feedback over the Southern California Current System that reinforced surface warming near the coast in 2014 (Myers et al., 2018;Schmeisser et al., 2019;Zaba & Rudnick, 2016). Below the mixed layer, anomalously warm and salty water was detrained to denser and deeper isopycnals, reaching depths of 140 m beginning in 2014 (Jackson et al., 2018). These subsurface anomalies lingered through at least 2018, long after the initial onset of anomalous atmospheric forcing in late 2013.

A similar situation played out during the summer of 2019 when a resurgence of Blob-like surface conditions intensified in the NE Pacific. Weakened surface wind speeds, driven by atmospheric teleconnections associated with SST anomalies in the tropical Pacific, resulted in reduced evaporative heat loss from the ocean to atmosphere and limited wind-driven mixing, resulting in a MHW off the U.S. West Coast (Amaya et al., 2020). Increased shortwave radiation and a positive SST-cloud feedback helped to maintain the MHW over an exceptionally shallow summertime mixed layer (Amaya et al., 2020). Here, we show evidence for the role of salinity anomalies in increasing upper ocean stability and describe the propagation and persistence of the 2019-2020 NE Pacific MHW in the subsurface.

In this study, we examine the connections between surface MHWs and the subsurface structure of temperature, salinity, and density by analyzing objectively mapped monthly Argo data in the NE Pacific, comparing and contrasting the 2013-2016 and 2019-2020 MHWs. We characterize the spatiotemporal evolution of anomalous subsurface conditions and their connection to mixed layer properties from January 2004 through June 2020, and we quantify the change in water mass properties and ocean heat content anomalies within and below the mixed layer. Understanding the subsurface evolution and persistence of MHWs gives insight into the potential predictability and reemergence of these events in the future, where a trend toward shallower summertime MLDs is expected to increase the likelihood and intensity of MHWs in the North Pacific (D.J. Amaya, personal communication, October 2, 2020). The persistence and potential reoccurrence of MHWs could result in long-lasting impacts on the health of marine ecosystems, especially in the subsurface where the effects of warming on marine life (i.e., thermal stress) can persist for years (Cavole et al., 2016).

Data

We analyze monthly mean SST maps from the Optimum Interpolation SST version 2 (OISSTv2) data set on a 0.25°longitude by 0.25°latitude global grid from 1982 through present (Reynolds et al., 2002(Reynolds et al., , 2007)). These SST maps are generated from a blend of satellite (Advanced Very High Resolution Radiometer only), ship, buoy (both moored and drifting), and Argo float data. The satellite data are interpolated to fill gaps and are bias corrected with reference to buoys to account for platform differences. We use the OISSTv2 data set as it incorporates in situ observations, offers complete global coverage, and spans almost 40 years.

We also analyze monthly mean fields from January 2004 through June 2020 from the updated Roemmich-Gilson Argo Climatology (Roemmich & Gilson, 2009; hereafter RG09) to examine the vertical structure of temperature, salinity, and density anomalies associated with MHWs. Argo is a global network of autonomous profiling floats that continuously measures the temperature and salinity of the upper 2,000 m of the ocean. The Argo program began in 1999 and now consists of over 3,800 active floats and more than 2 million hydrographic profiles reported thanks to a coordinated effort from dozens of countries worldwide (Jayne et al., 2017). Archived and near-real-time float data are made publicly available (http://sio-argo.ucsd.edu/ RG_Climatology.html) and are incorporated into monthly maps on a 1°longitude by 1°latitude grid beginning in January 2004 when the global array had at least 1,000 floats and first approached sparse global coverage (RG09). These maps are made in 58 pressure layers with the shallowest centered on 2.5 dbar and the deepest on 1,975 dbar, with finer resolution near the surface (e.g., spaced 10 dbar apart from 10 to 170 dbar). The 2.5-dbar monthly temperature anomalies in RG09 closely track the monthly OISSTv2 anomalies in the NE Pacific, capturing large-scale spatial and temporal variability.

In addition to the mapped temperature and salinity versus pressure fields from RG09, we also analyze 19,697 quality-controlled Argo profiles in the NE Pacific (35.5-51.5°N, 135.5-154.5°W; box in Figure 1) to compute the MLD from January 2004 through June 2020 using the density algorithm of Holte and Talley (2009). The sampling frequency from Argo in the NE Pacific (35.5-51.5°N, 135.5-154.5°W) steadily increases from the early 2000s, achieving over 1,000 profiles per year starting in 2012 (Figure S2). These profiles were downloaded from one of the two Argo Global Data Assembly Centers (https://nrlgodae1.nrlmry.navy.mil/argo/ argo.html) in August 2020.

Analysis

We define MHWs locally when SST exceeds the monthly climatological 90 th percentile for at least a month using monthly data from January 2004 through June 2020. Our definition for MHWs is similar to that proposed in Hobday et al. (2016) with modifications in the length of the climatological period and in the minimum event duration. Owing to the prominence and persistence of the 2013-2016 and 2019-2020 MHWs, our definition highlights the same large-scale features described in previous studies using daily data (e.g., Fewings & Brown, 2019;Gentemann et al., 2017).

Before analyzing the RG09 data set, we fit temperature and salinity at each spatial point to a mean, trend, annual, and semiannual harmonics using least squares regression from January 2004 through June 2020.

We then remove the mean, annual, and semiannual harmonics (but not the trend) to generate anomalies. Following MHW conventions (e.g., Hobday et al., 2016), we choose to retain the warming trend in the analysis using a fixed climatology computed over the entire record. Furthermore, the trend would not be accurately estimated over such a short period and would be extremely biased by the 2013-2016 and 2019-2020 MHWs at one end of the time series. Finally, detrending would effectively remove part of the strong MHW signal that we observe toward the latter end of the record. We therefore retain it. Next, we smooth the anomalies and the regression coefficients with a 5-month Hanning filter and then a 6°latitude × 6°1 0.1029/2020GL090548

Geophysical Research Letters

longitude locally weighted scatter plot smoothing (LOESS; Cleveland & Devlin, 1988) filter to reduce mesoscale signals that are retained in the RG09 maps. We then reconstruct the total smoothed in situ temperature and practical salinity maps using the smoothed anomalies and smoothed model coefficients. We apply the thermodynamic equation of seawater (Intergovernmental Oceanographic Commission, SCOR, & IAPSO, 2010) to compute the absolute salinity (S A ) and conservative temperature (Θ) at each space and time grid point. Using S A and Θ, we also compute the potential density anomaly (σ θ ) with reference to 0 dbar, expressed as a particular potential density minus 1,000 kg m -3 . The potential density represents the density a fluid parcel would acquire if it were brought adiabatically to the sea surface, thus eliminating the density dependence on pressure. We also map the RG09 fields of S A , Θ, and pressure (P) to a vertical density coordinate, σ θ . We compute anomalies in S A , Θ, and P in σ θ coordinates, as well as S A , Θ, and σ θ in P coordinates, by removing the monthly means of these quantities across the entire 198-month time series at each spatial point and for each vertical coordinate system (σ θ and P) to get the anomalies. We describe changes in S A , Θ, and P on an isopycnal (25.4 kg m -3 ) that may outcrop during winter. The properties of outcropped isopycnals are easily modified through air-sea interactions that may drive surface MHWs. Once isopycnals subduct below the mixed layer, their properties can be modified through mixing and/or lateral advection, which is usually less effective than direct air-sea heat and freshwater exchange.

We examine the ocean heat content anomaly (Q′) within the mixed layer (10-90 dbar), thermocline (100-180 dbar), and just below the thermocline (200-280 dbar). These layers of equal thickness are chosen based on the vertical profiles of subsurface temperature in the NE Pacific (Figure 4b). They typify the surface, pycnocline, and interior ocean in the region, allowing for the distinction of the changes in Q′ with depth. We

, where g = 9.8 ms -2 is the acceleration due to gravity, c p = 3991.8680 J kg -1 K -1 is the standard specific heat of seawater when using Θ, Θ′ is the conservative temperature anomaly, and ∫ dp is the integral over each of these three 80-dbar thick layers.

We apply the Holte and Talley (2009) density algorithm to 19,697 Argo float profiles in the NE Pacific (35.5-51.5°N, 135.5-154.5°W; box in Figure 1) to estimate monthly MLDs from January 2004 through June 2020. This method searches for the depth at which the density increases by 0.03 kg m -3 relative to a near-surface reference level.

We quantify the bulk stratification of the upper ocean using the Brunt-Väisälä frequency squared N 2 ¼ -

Here, dρ dz is the change in potential density with reference to 0 dbar between 2.5 and 200 dbar. Larger values of N 2 correspond to greater upper ocean stratification-a more stable water column. We compute anomalies in N 2 , again with respect to monthly long-term means, to quantify the change in the stratification of the upper ocean due to MHW variations in both Θ and S A .

To further examine the relationships among Θ, S A , and σ θ , we examine Θ -S A diagrams with contours of constant density and spice to show changes in water mass properties between different MHW years in the NE Pacific. Spice quantifies Θ -S A variations along isopycnals (Munk, 1981), where warm/salty anomalies are spicy and cool/fresh anomalies are minty. We compute spice following McDougall and Krzysik (2015) using a potential density with reference to 0 dbar. Isopycnal variations in spiciness can be used to describe MHW impacts on isopycnal water mass properties in density units.

Results

Anomalies in Θ -S A on isopycnals can be tracked following the surface evolution of SST anomalies during MHWs and can either be warm/salty (spicy) or cool/fresh (minty), such that the density of that isopycnal does not change (Movie S1). The winter-intensified 2013-2016 MHW had spicy anomalies on 25.4 kg m -3 , which lagged the spatiotemporal evolution of SST anomalies within the MHW (Movie S1, hatching in Figure 1). For example, surface MHW conditions appeared nearshore by late 2014 and began to fade as early as 2015, whereas subsurface spice anomalies did not appear along the coast until winter 2015 and persisted into 2016 (Movie S1). By comparison, summer Θ -S A anomalies in 2019 lacked the transitional nature of the 2013-2016 MHW, yet they were much more widespread. Minty anomalies on 25.4 kg m -3 encompassed 10.1029/2020GL090548

Geophysical Research Letters

nearly the entire Gulf of Alaska from late summer 2018 through summer 2020, while spicy anomalies lingered off the coast between Baja California and Hawai'i (Figure 1 and Movie S1).

Positive stratification (N 2 ) anomalies occurred for both the 2013-2016 and 2019-2020 MHWs; however, they were much greater in 2019 (Figure 1 and Movie S1). Warm and fresh near-surface anomalies in 2019 decreased density and increased the stratification (Figure 2), whereas in 2013-2016 the near-surface density reduction from a warm anomaly was partially offset by a salty anomaly. The increase in pressure along 25.4 kg m -3 beginning in 2018 reflects an increase in stratification even before the onset of the 2019 MHW (Figure 2). The 2019 large and positive stratification anomaly likely inhibited the surface MHW from penetrating as deeply as the 2013-2016 MHW and furthermore may have enhanced the surface buildup of heat.

Prior to 2013, two other noteworthy MHWs occurred in the NE Pacific in 2004-2005and 2008-2009 (Figure 2) (Figure 2). Warm subsurface Θ anomalies during these MHWs extended to depths beyond 100 dbar, and anomalies at 25.4 kg m -3 were spicy, similar to that of the 2013-2016 event (Figure 2). Warm and salty anomalies reduced subsurface density and increased the stratification of the surface layer. The 2004-2005 MHW was more stratified than the 2008-2009 event owing to the larger surface density anomaly (Figures 2e, 4b, and 4c).

The simultaneous change in temperature from 0-200 dbar in 2008-2009 could have resulted from isopycnal heave, as indicated by the downward deflection of 26.3 kg m -3 (Figure 2a). Heave can occur in response to Ekman pumping due to wind stress curl that depresses the main thermocline (Bindoff & McDougall, 1994) or from other dynamic features such as large-scale Rossby waves (Xie et al., 2016) or eddies (Pegliasco et al., 2015). Positive pressure anomalies on 26 kg m -3 indicate a deepening of the thermocline in 2008-2009 at approximately 130 dbar (Figure 2f ). These vertical isopycnal motions are nearly adiabatic.

As seen from the conservation of water mass properties on the isopycnal (Figures 2b and 2d), there is little exchange of heat or salinity with the surrounding environment. As a result, warm and fresh anomalies in 2008-2009 that occurred along the 150-200 isobars were negligible on 26.3 kg m -3 , which ranges from 150 to 200 dbar (Figure 2). (lower dashed), and 26.3 kg m -3 (lower solid) isopycnal surfaces vary with pressure (a, c, e), but are constant when plotted against density (b, d, f ).

10.1029/2020GL090548 3a). The warmest wintertime temperatures occurred in 2013/2014 when Θ -S A variations were confined to lighter isopycnals (<26 kg m -3 ).

Winter spice anomalies in 2013/2014 likely mixed to denser isopycnals in the permanent halocline by summer, as can be seen along 25.6 kg m -3 during the summers of 2014 and 2015 (Figure 3b). By summer 2016, spice anomalies within the permanent halocline returned to near normal; however, the seasonal thermocline remained anomalously spicy. Spice anomalies during the summer 2019 MHW were minty compared to average. Minty conditions in June-July-August of 2019 were greatest within the seasonal thermocline above 25.5 kg m -3 (Figure 3b). As a consequence, the near-surface Θ -S A properties were much lighter compared to 2014-2016, both in winter and summer seasons. Minty conditions persisted into the winter of 2019/2020.

A connection between the evolution of surface and subsurface anomalies was a recurring theme during recent 2013-2016 and 2019-2020 NE Pacific MHWs and is visible in both Figures 2 and 4. To quantify the time lags associated with the penetrations of surface anomalies into the subsurface, we compute the lagged cross correlation for Θ and S A on isobars and isopycnals with values at 2.5 dbar and 25.7 kg m -3 , respectively. Significant positive correlations between surface and subsurface Θ -S A anomalies increase with positive lag and density between 25.7 and 27 kg m -3 . For example, the maximum cross correlation on 26.3 kg m -3 occurs at 6 months positive lag (Figure S3). On the other hand, subsurface Θ anomalies (between 150 and 220 dbar) are most strongly correlated with the surface conditions for positive lags of 1-2 years, while subsurface S A correlations peak at 6-12 months positive lags (Figures S3 and S4).

The downward progression of surface Θ and S A anomalies suggests that the NE Pacific Ocean is capable of maintaining long-term memory of surface MHWs. One measure of memory is the heat content anomaly, Q′, evaluated here over equal thickness subsurface layers. The largest Q′ values occur within the seasonally varying mixed layer (10-90 dbar) where temperature fluctuations are the strongest (Figure 4). The largest positive anomalies are present during the 2013-2016 MHW. After a period of strong cooling, Q′ steadily increased beginning in 2018 through present. Prior to 2013 there were two smaller MHWs that occurred in 2004-2005 and 2008-2009 that also had small gains of heat content. Evaluating Q′ over layers spanning the pycnocline (100-180 dbar) and interior (200-280 dbar) reveals the persistence of Θ anomalies below the surface temperature variability. Once Θ -S A anomalies get into the subsurface, their properties are nearly conserved even after the surface cools (Figure 4).

Geophysical Research Letters

An increase in upper ocean heat content can affect the stability of the upper ocean. The depth of the mixed layer also shoals, which can be seen during the winters of 2013/2014 and 2014/2015 (Figure 4). The increase in stratification reduces entrainment of cool water from below and can exacerbate warming by reducing the thickness of the surface layer that accepts heat from the atmosphere, making the surface ocean easier to warm. The upper ocean stratification anomaly was noticeably higher (large N 2 anomaly values) in 2014-2015, with the largest values occurring in 2019 (Figure 4c). The very high values in 2019-2020 arise from the anomalously fresh near-surface conditions during that MHW. 10.1029/2020GL090548 Geophysical Research Letters SCANNELL ET AL.

Discussion

This study examines 21 st Century MHWs in the NE Pacific based on gridded SST data and also the evolution of subsurface Θ -S A anomalies from Argo on both isobars and isopycnals during the 2013-2016 and 2019-2020 NE Pacific MHWs. Upper ocean salinity was anomalously fresh in the Gulf of Alaska during the 2019-2020 MHW, which greatly increased the buoyancy of the surface layer. Indeed, there was a net freshwater input from precipitation as can be seen in the 2018 precipitation anomaly in the Gulf of Alaska (Yu et al., 2019) that likely contributed to the decrease in surface salinity (Reagan et al., 2019). The resulting increase in stratification during 2019-2020 likely contributed to the decrease in the depth (and density) to which water property anomalies from this event were detrained and in places subducted. The confinement of warm anomalies to the near-surface likely enhanced the MHW's intensity.

There are several dynamical pathways by which surface MHW anomalies in the NE Pacific could reach the subsurface, by means of detrainment, diabatic subduction (Jackson et al., 2018), lateral advection (Chao et al., 2017;Zaba et al., 2020), and/or adiabatic isopycnal heave. Subduction occurs in subtropical regions after temperature anomalies within the deep wintertime mixed layer detrain as a result of the mixed layer retreating in late spring. During the 2014 and 2015 spring transition of the MLD, subsurface warming occurred along both isopycnals and isobars below the mixed layer, suggesting that diabatic vertical or horizontal mixing could play a role in the penetration of MHW anomalies within the seasonal pycnocline. Indeed, Zaba et al. (2020) attribute positive subsurface heat content anomalies within the California Undercurrent to an increase in poleward heat transport from the tropics in September 2015. Alternatively, subsurface warming that occurs primarily on isobars and not on isopycnals was likely the result of isopycnal heave, defined as the downward deflection of a potential density surface. We speculate that heave is most likely responsible for the near-simultaneous appearance of anomalies below 150 dbar, for example, during the 2008-2009 MHW; however, the exact mechanisms of heave (e.g., from Ekman pumping due to wind stress curl) are not investigated here.

Once surface MHW anomalies are detrained out of the deep wintertime mixed layer, they may propagate downward. The lag associated with the vertical propagation of surface anomalies causes the subsurface heat content to remain anomalously high even after surface conditions return to normal. This persistence of subsurface heat and the possible seasonal reemergence of surface anomalies could in fact help supercharge the occurrence of multiyear events. As future warming trends favor a more stratified upper ocean (Li et al., 2020), we expect that detrainment out of the mixed layer may become less effective in storing MHW anomalies in the subsurface and therefore further amplify surface warming. This possibility is concerning owing to the impacts that accumulated heat stress and stratification have on pelagic marine ecosystems and primary production (Cavole et al., 2016;Jacox et al., 2016;Smale et al., 2019).

Mixed layer heat budgets are frequently used to diagnose the drivers of surface warming associated with MHWs; however, the influence of salinity and subsurface water mass properties are often overlooked (Holbrook et al., 2020). Using data from the global Argo array, this study motivates complementary analyses on the role of salinity and subsurface Θ -S A anomalies to better understand the ocean's role in the persistence and evolution of long-lived events. Further investigation into the drivers of salinity anomalies and their role in the development of NE Pacific MHWs would appear to be a fruitful avenue of future research. Analysis of the full 4-D heat budget using high-resolution numerical models could be undertaken to investigate the local mechanisms of subsurface warming.