Drivers of Marine Heatwaves in the Northwest Atlantic: The Role of Air–Sea Interaction During Onset and Decline

Marine heatwaves (MHWs) are increasing in duration and intensity at a global scale and are projected to continue to increase due to the anthropogenic warming of the climate. Because MHWs may have drastic impacts on fisheries and other marine goods and services, there is a growing interest in understanding the predictability and developing practical predictions of these events. A necessary step toward prediction is to develop a better understanding of the drivers and processes responsible for the development of MHWs. Prior research has shown that air–sea heat flux and ocean advection across sharp thermal gradients are common physical processes governing these anomalous events. In this study we apply various statistical analyses and employ the self-organizing map (SOM) technique to determine specifically which of the many candidate physical processes, informed by a theoretical mixed-layer heat budget, have the most pronounced effect on the onset and/or decline of MHWs on the Northwest Atlantic continental shelf. It was found that latent heat flux is the most common driver of the onset of MHWs. Mixed layer depth (MLD) also strongly modulates the onset of MHWs. During the decay of MHWs, atmospheric forcing does not explain the evolution of the MHWs well, suggesting that oceanic processes are important in the decay of MHWs. The SOM analysis revealed three primary synoptic scale patterns during MHWs: low-pressure cyclonic Autumn-Winter systems, high-pressure anti-cyclonic Spring-Summer blocking, and mild but long-lasting Summer blocking. Our results show that nearly half of past MHWs on the Northwest Atlantic shelf are initiated by positive heat flux anomaly into the ocean, but less than one fifth of MHWs decay due to this process, suggesting that oceanic processes, e.g., advection and mixing are the primary driver for the decay of most MHWs.


INTRODUCTION
The Northwest (NW) Atlantic continental shelf stretches along the coast of North America from the Middle Atlantic Bight (MAB) in the south to the Labrador Shelf in the north. Over the continental shelf, a continuous flow brings cold and fresh high-latitude waters equatorward (e.g., Loder et al., 1998). Off the shelf break, the Gulf Stream transports warm salty water poleward and separates from the coast near Cape Hatteras. Between the Gulf Stream and the shelf, large Gulf Stream meanders and warm-core rings frequently form and influence the shelf environment (e.g., Garfield and Evans, 1987;Joyce et al., 1992;Ryan et al., 2001;Wei and Wang, 2009;Gawarkiewicz et al., 2012;Chen et al., 2014b;Zhang and Gawarkiewicz, 2015). As with most western boundary current systems, the Gulf Stream and the adjacent NW Atlantic shelf-slope region is an area of the ocean experiencing a rate of warming much higher than the global average (Wu et al., 2012;Forsyth et al., 2015;Pershing et al., 2015). The presence of the Jet Stream over the NW Atlantic also plays a role in the formation of anomalously warm bodies of water (e.g., Chen et al., 2014a). While the longterm change of this large-scale atmospheric feature is being debated (e.g., Francis and Vavrus, 2012;Barnes, 2013;Francis and Vavrus, 2015;Blackport and Screen, 2020), jet stream related activity such as high pressure blocking significantly impacts the weather and climate system (e.g., Santos et al., 2013) in the NW Atlantic.
The important fisheries that the NW Atlantic has historically supported have been in flux due, in part, to rapid changes in temperature. The cod fishery was one of the largest in the NW Atlantic until it's collapse by 1993 (Myers et al., 1997). This was initially attributed solely to overfishing but a more recent analysis supports the theory that rising temperatures also played an important role (Pershing et al., 2015). This gap in the marine economy was largely filled by the lobster and snow crab fisheries (Mills et al., 2013). The lobster fishery received a boost from a record-breaking warm ocean temperature extreme, known as a marine heatwave (MHW) in 2012 (Mills et al., 2013), but this same event had a negative impact on the local snow crab fishery (Zisserson and Cook, 2017). This rapid increase in warmth appears to have been most beneficial to lobsters just north of the United States -Canada border, leading to some tensions between the two countries (Mills et al., 2013). These tensions are likely to continue as the center of the lobster fishery is slowly progressing northward (Greenan et al., 2019), tracking the rapid increases in bottom water temperatures of the Gulf of Maine (Kavanaugh et al., 2017). This measured northward shift is punctuated by extreme temperature events, like the 2012 MHW, which has prompted the necessity to understand what is driving these events.
Marine heatwaves have been increasing in duration and intensity for decades (Oliver et al., 2018). Because these events are driven primarily by increases to the mean sea surface temperature (SST) in a region, and not increases in variability (Oliver, 2019), MHWs are projected to increase indefinitely into the future (Frölicher et al., 2018;Oliver et al., 2020) as the anthropogenic warming of the planet slowly ratchets up mean global SST (Pachauri et al., 2014). While the anthropogenic warming signal may be responsible for the overall increase in MHWs, it is not itself the single explanation for what drives the onset or decline of individual events. There are commonly two primary drivers of MHWs: atmospheric forcing (i.e., air-sea heat flux through the ocean surface) and ocean advection (i.e., horizontal currents crossing temperature gradients). There is a rapidly growing body of literature covering the drivers of individual events and many of the major MHWs have been well studied [see Holbrook et al. (2019) for a review]. The first MHW to garner global attention was the 2003 event in the Mediterranean that was shown by Olita et al. (2007) to have been caused by anomalously warm air over Europe in combination with a drastic reduction in wind stress, leading to a reduction of upward heat flux. The very intense 2011 MHW off the southwest coast of Australia was primarily forced by the abnormal advection of the Leeuwin Current onto the coast (Feng et al., 2013), with the remainder due to air-sea heat flux (Benthuysen et al., 2014). The NW Atlantic 2012 MHW that had such a large impact on the lobster fishery was primarily forced by heat flux from the atmosphere that was heightened by an anomalous northward movement of the Jet Stream (Chen et al., 2014a(Chen et al., , 2015. The 2014 Pacific "blob" was similarly caused by anomalous air-sea flux presumably associated with jet stream activity and ocean advective subsequently played a role in the evolution of this event (Bond et al., 2015). These earlier studies of multiple high profile MHWs have revealed the complex interplay of atmospheric and oceanic processes in driving individual MHWs. In the meantime, there is a need for an overall understanding of the drivers of all historical MHWs, especially their onsets and declines, respectively.
In this study, we focus our effort on understanding the past MHWs on the NW Atlantic continental shelf, a region that has been experiencing many recent changes including rapid warming (Forsyth et al., 2015;Pershing et al., 2015;Chen et al., 2020) and intense marine heatwaves (Chen et al., 2014b(Chen et al., , 2015. We relate a set of candidate physical processes to the onset and decline of these MHWs. Much of the philosophy for this approach comes from Chen et al. (2016), in which the authors were able to illustrate the parts of the heat budget that were most likely driving the anomalous heat content in the surface of the ocean for the NW Atlantic from the years 2003 to 2014. The analysis in this paper seeks to build on this methodology by analyzing additional variables that do not appear explicitly in the mixed-layer heat budget equation, as well as by examining the spatial and seasonal differences over the NW Atlantic shelf. The period of study is also extended from 1993 to 2018. There are four layers of focus in the methodology. The first is the simple comparison of the magnitude of the change in SST anomalies (SSTa) during the onset and decline of MHWs against the magnitude of change in heat flux variables. The second is to use root mean square error (RMSE) to determine the heat flux variables most closely related to the change in SSTa. The third is to use correlations between SSTa and a range of atmospheric and oceanic variables during MHWs to determine which of these non-heat flux variables relates most closely. The fourth is to use a self-organizing map (SOM) to visualize the most common synoptic air-sea states during MHWs. This last step is performed because we want to know not just which specific variables are the most important drivers, but also if MHWs arise due to synoptic patterns such as blocking highs or anomalous Gulf Stream meanders. A discussion is provided on the drivers of MHWs and the importance of the differences found between regions and seasons. This paper concludes with a summary of the primary drivers of MHWs in the NW Atlantic.

Study Area
This study used the geographic regions defined by Richaud et al. (2016) to divide the NW Atlantic shelf into sub-regions based on their climatological SST and sea surface salinity (SSS) ( Figure 1A). We modified the scheme in Richaud et al. (2016) by defining a Cabot Strait region on the shelf as distinct from the more isolated Gulf of St. Lawrence region. We also excluded the Labrador Shelf region as we are primarily concerned with the Atlantic coast south of the subpolar gyre. The pixels for the different data products (see section "Data") found within each region polygon were spatially averaged to create a single time series per region per variable.
For the SOM analysis [see section "Self-Organizing Map (SOM)"], which uses fully spatially resolved fields of ocean and atmosphere variables, the extent of the study area was determined by expanding out two degrees of longitude and latitude from the furthest edges of the polygons described above ( Figure 1A). This encompasses broad synoptic scale variables that may be driving MHWs in the study regions, but should not be so broad as to begin to account for regional teleconnections, which are currently beyond the scope of this project. We also excluded as much of the Labrador Sea as possible, and so only extended the northern edge of the study area by 0.5 degrees of latitude from the northernmost point of the Gulf of St. Lawrence polygon.

Surface Mixed-Layer Heat Budget
A mixed-layer heat budget can be used to relate the tendency of surface mixed-layer temperatures to a set of physical processes (e.g., Moisan and Niiler, 1998): The left-hand side of the equation shows the rate of change of the temperature of the mixed layer (T mix ) with time (t). The right-hand side of the equation consists of three terms: airsea heat flux into the mixed layer, convergence of heat in the mixed layer due to horizontal advection, and a residual term that includes the additional physical processes not explicitly accounted for here including vertical advection, entrainment, and mixing [see Oliver et al. (2020) for a review of the mixed-layer heat budget for MHW studies]. The air-sea heat flux term shows that temperature changes are driven by the net downward heat flux (Q net ) scaled by the density of seawater (ρ), the specific heat content of seawater (c p ) and the mixed layer depth (H). The variable Q net itself is composed of the sum of four individual heat flux variables: latent heat flux (Q lh ), sensible heat flux (Q sh ), longwave radiation (Q lw ), and shortwave radiation (Q sw ). Throughout, the Q sw variable in this study has been corrected for with the shortwave radiation that passes through the bottom of the mixed layer [i.e., Q sw (−H ) ]. This loss is calculated following the radiation decay profile from Paulson and Simpson (1977), using a water type of 1B representative of the NW Atlantic. The horizontal advection term shows that temperature changes are also driven by average horizontal velocities in the mixed layer (u mix ) acting across horizontal temperature gradients (∇ h T mix ).

FIGURE 1 | (A)
The coastal regions from Richaud et al. (2016) shown as colored polygons. The numbers within each polygon show the total count of marine heatwaves (MHWs) detected there during the 1993-2018 period. The total counts across all regions by season are shown with white labels above the legend. More information on the application of the heat budget equation for understanding MHWs may be found in Benthuysen et al. (2014), Chen et al. (2015), and Oliver et al. (2017). The rate of change of temperature due solely to the air-sea heat flux term is given by The right-hand side is a linear combination of the effect due to each of the air-sea heat flux components (Q x , where x = lw, sw, lh, or sh) and so the temperature change during a MHW due to one of these components is given by the following time integral where t 1 is the start time of a phase of an MHW (i.e., onset or decline; see section "Agreement Metrics"), and Q x is the given air-sea heat flux variable. In order to generate anomalies of T Qx we first generated anomalies of the combined time series Q x /H, and applied the equation above to those. We used values of ρ = 1024 kg/m 3 and c p = 4,000 J/(K kg). Using this method we can calculate the time series of SST change during a MHW due to the effect of each air-sea heat flux component, i.e., T Qsw , T Qlw , T Qsh , and T Qlh . All of the variables chosen below allow for the comparison of surface temperatures against various aspects of the other terms in the mixed-layer heat budget equation.

Data
We used the NOAA OISST v2.1 product to obtain daily fields of SST in our study region. The OISST v2.1 product is a remotely sensed infrared (AVHRR) SST retrieval optimally-interpolated onto a daily 1/4-degree global grid (Reynolds et al., 2007;Banzon et al., 2016;Banzon et al., 2020). It is a level four (L4) product, meaning that any gaps in the daily retrieval are filled via interpolation and with the aid of any nearby in situ data from ships, buoys, etc. The atmospheric variables in this study were taken from the ERA5 reanalysis product (Copernicus Climate Change Service [C3S], 2017). ERA5 provided hourly atmospheric fields on a 1/4degree grid which were then made into daily means. The variables obtained included longwave (Q lw ) and shortwave (Q sw ) radiation, latent (Q lh ) and sensible (Q sh ) heat flux, mean sea-level pressure (MSLP), total precipitation (P), total evaporation (E), total cloud cover (TCC), wind at 10 m (U 10 , V 10 ), and air temperature at 2 m (T air ). The four air-sea heat flux variables were summed together to give Q net ; air-sea heat flux is expressed as positive downward. The other three variables created were precipitation minus evaporation (P-E), wind speed (W 10 ) and direction ( 10 ). The daily means for the non-heat flux variables from ERA5 were centered on noon, meaning that the daily averages were taken from midnight to midnight for a given calendar day. However, since the air sea heat flux term represents a rate of change in temperature, and thus it's integral is comparable to temperature, the daily means in air-sea heat fluxes were calculated centered on midnight, i.e., leading the other variables by 12 h. The integral from midnight to midnight would thus result in daily values centered on noon, comparable to the definition of temperature. This allows the heat flux variables to be compared against SST as an integral of the driving of the development of SST over time.
The oceanic variables used in this study were taken from the GLORYS12V1 (hereafter GLORYS) product, which are provided by the E.U. Copernicus Marine Service Information. This is a global ocean reanalysis product with 50 levels on a daily 1/12-degree grid. The variables obtained were: sea surface height (SSH), sea surface salinity (SSS), surface currents (U, V), mixed-layer depth (MLD in text; H in equations), and bottom temperature (T b ). From these variables the surface current speed (W) and direction ( ) were calculated. When finding the average of these variables per region (Figure 1A), the full 1/12-degree resolution was used. In order to compare the GLORYS data with ERA5 and OISST for the SOM analysis [see section "Self-Organizing Map (SOM)"], these variables were averaged onto a shared 1/4-degree grid.
The 1/4-degree grids of the OISST and the ERA5/GLORYS differ slightly, with the pixels staggered by 1/8 degree. For consistency, the 1/4-degree ERA 5 and GLORYS products were regridded to match the OISST grid. This was not done for the 1/12-degree GLORYS variables when finding the pixels within each polygon ( Figure 1A) because of the much higher resolution. The full time-period considered in the analysis was 1993 -2018 as this was the longest period of overlap between these datasets.
Because we are interested in observed MHWs from the historic record we decided to use an observations-based product (NOAA OISST) rather than a reanalysis (e.g., GLORYS). We conducted the full analysis (see below) using both, and found that MHWs in GLORYS are generally longer and less intense, which is not surprising as numerical models tend to be less capable of capturing higher frequency and wavenumber signals. Because we are interested in investigating MHWs at a higher temporal resolution including both short and long events, we chose to use the daily OISST rather than the GLORYS reanalysis product.

Marine Heatwaves (MHWs)
We use the Hobday et al. (2016) definition for MHWs. This defines MHWs as a period of five or more days at or above the 90th percentile threshold at any given location, based on a seasonally varying climatology at a daily interval. The MHWs for each region in the study area were detected within the OISST dataset using the heatwaveR implementation of this MHW definition (Schlegel and Smit, 2018; Figure 1B). Note that large and intense events may occur simultaneously over multiple regions and may have long-lasting impacts on the background state. As a result, a large and intense event may be detected as multiple MHWs in different regions and/or as multiple MHWs close to each other in time. However, from a heat budget perspective, the onset and decay of these individual MHWs across space and time may not be necessarily driven by the same processes. Our methodology of targeting all MHWs that satisfy the definition allows for an analysis of any potential differences in the physical drivers of each individual event (see below) or in the synoptic scale forcing of these events across regions [see section "Self-Organizing Map (SOM)"].
The climatological period of 1993 to 2018 was chosen as it was the longest period shared across all data products in this study. This 26-year period differs from the recommended minimum length of 30 years (Hobday et al., 2016), but we do not expect this to bias our results significantly (Schlegel et al., 2019).
Once an MHW is defined it can be characterized by a set of metrics. We have focused on the following four metrics for this study: (1) the duration (D) of a MHW is the number of days from the start to the end of when the temperature anomaly is in excess of the 90th percentile threshold; (2) the maximum intensity (i max ) is the highest temperature ( • C) in excess of the seasonally expected daily value during the MHW, i.e., the maximum temperature anomaly; (3) the mean intensity (i mean ) is the average of the temperatures ( • C) above the seasonally expected daily values throughout the MHW; (4) the cumulative intensity (i cum ) is the time integral of the temperature anomalies above the seasonally expected daily values throughout the MHW ( • C days). The full suite of explanations for all of the MHW metrics may be found in Table 2 of Hobday et al. (2016).
For grouping purposes throughout the methodology it was necessary to be able to assign a single season of occurrence to each MHW, regardless of how long they may last. It was decided that the peak date (day on which i max occurs) of each MHW would dictate its season of occurrence. All seasons are 3 months long with Winter consisting of: January, February, and March.

Agreement Metrics
We compared each candidate forcing variable outlined in the data sub-section with MHWs by constructing a time series of the forcing variable and SSTa during the events. Given these time series we can compare the forcing variable and SSTa during an MHW in three ways: (1) compare the magnitude of change in T Qnet and SSTa, (2) calculate the RMSE between SSTa and each of the four heat flux variables within T Qnet , and (3) find the correlation between the non-heat flux variables and SSTa.
For all variables we calculated anomalies from a seasonal climatology. We used the same climatology calculation method outlined in Hobday et al. (2016) as was used with SST for MHW detection. Those climatologies were subtracted from the observed values in each variable to create the anomaly time series and all reference to variables hereafter is to these anomalies. The climatology period used matched that for SST, i.e., 1993SST, i.e., -2018 Some variables, like SSS, change daily, so comparing them to daily changes in SSTa is straight forward. Some variables, like MLD, may change initially but then persist at a stable level. Or other variables, like W 10 or MSLP may have a cumulative effect on SSTa. This makes comparisons to daily changes in SSTa less meaningful. To address this issue we created the following cumulative anomaly time series: TCC cum , MSLP cum , W 10cum , P-E cum , and MLD cum . This was done by taking the anomaly value on the start date of the particular phase of an MHW (i.e., onset or decline), and adding the following day cumulatively until the end of that phase. For the onset phase of a MHW this is from the start date to the peak of the event, and for the decline phase this is from the peak to the end of the event. The peak date of an event is the day on which the maximum temperature anomaly occurs. Note that these values are not designed to quantify the daily changes in SSTa, but rather are used to provide a fluctuating time series that can be directly correlated with SSTa. This is because these variables alone may be expected to be proportional to the temperature tendency, so these cumulative values are more likely to be proportional to the daily temperatures.

Magnitude of Change in T Qnet
The magnitude of change in T Qnet during the onset of MHWs was calculated by subtracting the value of T Qnet at the start of the event (which is always 0) from the value of T Qnet at the peak. The magnitude of change during the decline of an event was calculated by subtracting the value of T Qnet at the peak of the event from the value of T Qnet at the end. The magnitudes of change in SSTa during the onset and decline of MHWs were calculated likewise so that they could be used to divide the magnitude of change in T Qnet , thus providing the proportion of change in SSTa attributable to T Qnet . As an example, during the onset of a given MHW SSTa increased from 0 to 3 • C at its peak, and T Qnet increased from 0 to 4 • C. This means that the magnitude of change for SSTa is 3 • C, 4 • C for T Qnet , and the proportion of change attributable to T Qnet is 1.3 (4/3). This tells us that anomalous net air-sea heat flux likely drove the onset of this event, and because the ratio is >1 other processes were working to compensate for the SSTa increase. Should the proportion of change attributable to T Qnet during the onset and/or decline of a MHW be at or below 0.5, we may conclude that air-sea heat flux was not the primary driver of the onset/decline of that event. The magnitude of change for a phase of a MHW (i.e., onset or decline) was not calculated if it was less than 3 days in length as this would prevent the calculation of the following metrics of agreement.

RMSE of T a and T Qx
Root mean square error is the square root of the mean of the squared differences between two sets of data and is used here as a measure of how well SSTa and T Qx time series agree during the onset and decline of MHWs. This is expressed as: where T a is SSTa, T Qx is a given heat flux variable, and N is the number of days in the given MHW phase being compared. The RMSE during the onset of a MHW is calculated from the first day of the event to the peak, and for the decline it is calculated from the peak to the end of the event. The resultant RMSE value shows what the average difference per day between an SSTa and T Qx time series is in • C for a phase of a MHW. The difference in magnitudes of change between SSTa and T Qnet show us how much of the temperature change at the peak and end of the MHWs could be attributed to T Qnet , but the RMSE calculations give us a more precise measure of how well the T Qx variables track the daily changes to SSTa. Because we are interested in MHWs that have been shown to be driven primarily by heat flux, we will only calculate RMSE for the phase of events where more than 50% of the change in SSTa could be attributed to T Qnet .

Correlations Between T a and Forcing Variables
Pearson correlations were used to determine which of the airsea variables (non-heat flux) showed the strongest relationship to SSTa during MHWs. These correlations were run between SSTa time series from the start to peak (onset) and peak to end (decline) of each MHW for the full suite of air-sea forcing variables (see section "Data"), and should provide us with an idea of which variables may be behaving similarly to the rise and fall of SSTa during the onset and decline of a MHW. These various correlations may be easily explained as: where r x is the correlation value for variable x. Variable x can be any of the non-heat flux variables described above (see section "Data").
Correlations are a commonly used statistic for finding how well two variables co-vary with one another, but they do not show an indication as to if the magnitude of variation is similar. Meaning, two variables could correlate quite well if they both increase at roughly the same rate, even if they are of entirely different units of measurement and the values themselves are far apart from each other. This is why RMSE is used for the T Qx variables, because it allows us to say more precisely how much of the change in one variable is reflected in the other. Unfortunately, this cannot be done for variables with different units like MLD vs. SSTa.
If the correlations between a variable and SSTa for all of the MHWs in a region/season were random, or at least very noisy, we would expect the distribution of r x values to look something like a normal distribution, with most r x values near 0, and increasingly fewer as we move out to the tails of the distribution (−1 and 1). If, however, there is a consistently strong relationship (positive or negative), the distribution of the resultant r x values should be strongly skewed toward one of the tails.

Self-Organizing Map (SOM)
The most representative air-sea synoptic patterns during MHWs were determined using a self-organizing map (SOM). A SOM is a method of clustering that differs from more traditional techniques, such as K-means clustering, in that it is better suited for working with highly dimensional data (e.g., gridded SST data), and after it clusters the data it organizes the clusters (hereafter referred to as nodes) into a grid that best reduces the stress in neighboring nodes (Hewitson and Crane, 2002). This technique has been seeing increasing use in climate studies (e.g., Cavazos, 2000;Johnson, 2013;Gibson et al., 2017), and has been used with MHWs (e.g., Schlegel et al., 2017;Oliver et al., 2018). The SOM algorithm from the yasomi R package was used in this analysis (Rossi, 2012).
We use a subset of the variables outlined in Section "Data" in order to keep the stress in the model as low as possible while still providing a representation of the air-sea states. The variables fed to the SOM were the anomalies of: SST, MLD, U, V, T air , U 10 , V 10 , and Q net . Note that this provides the same number of air and sea variables so as to provide an even representation for the model when it performs its clustering. To ensure all variables are weighted equally in the SOM they were each centered around their mean value then scaled by their standard deviation. The data from the three products were regridded to the 1/4 degree OISST grid resolution for the entire study area. Due to the size and complexity of these data, it was necessary to use principal component initialization with the SOM. This is determined by running a principal component analysis on the data and using the first two principal components to initialize the SOM (Akinduko et al., 2016). This ensured a quick convergence to a consistent final result each time the SOM was run.
As with the metrics of agreement, the SOM was run on anomaly values. This required that daily climatologies for each variable for each pixel be created using the same seasonally varying climatology method with the same 1993-2018 base period. From these anomalies mean synoptic states for each variable during each of the MHWs detected in each region were created. To do this, only the days in the anomaly datasets of the chosen variables during each of the observed MHWs were taken and time-averaged to create the synoptic states during each event. For example, if a MHW lasted from 1993-02-09 to 1993-04-15, the daily anomalies over that period for the entire study area would be averaged together in order to create the single synoptic anomaly state for each variable given to the SOM to represent that MHW. Once each MHW was clustered into a node, the stack of synoptic states for those events were averaged to produce the representative synoptic state per variable for each node. This methodology follows that used by Schlegel et al. (2017) and Oliver et al. (2018).
Experiments were conducted with different SOM node counts in order to arrive at a 12 node grid (4 × 3) as this produced important differences in the middle two columns that a 3 × 3 grid excluded. We tried a larger (4 × 4) grid but the additional nodes did not produce additional information while further complicating the interpretation of the SOM; smaller grid sizes were unable to capture the variation in the low-pressure systems and high pressure blocking patterns. In order to ensure that a sufficient grid size was being used an analysis of similarity (ANOSIM) was run for each of these experiments (Johnson, 2013). It was found that all of these experiments produced significantly different nodes at p < 0.001.

Marine Heatwaves (MHWs)
The occurrence of MHWs differs visually between regions ( Figure 1B) but with a general shared pattern of increased occurrence in more recent years. In fact, there are no MHWs detected in the first year of the analysis (1993) in any region and only four MHWs detected in the second year (1994) across all regions. Regular annual occurrences of MHWs began in 1999 with the highest count occurring in 2012 consistent with that being a known major MHW year for the NW Atlantic. Overall, the highest count of MHWs occurs in the MAB (57) and the lowest count on the Scotian Shelf (41; Figure 1A). The Scotian Shelf, however, had the longest, most intense MHWs on average, with the shortest and least intense events occurring in the Gulf of St. Lawrence (Table 1).
Across all regions, there were almost twice as many MHWs detected in Summer (105) as there were in Spring (58) and Winter (58; Figure 1A). MHWs detected in the winter were on average the longest but least intense (Table 1). It must be noted that while the Winter events were the least intense on average, the variance in cumulative intensities was the highest, meaning that there was also a greater probability for large events compared to other seasons. This result is consistent with previous findings that temperature in the MAB during wintertime has the largest interannual variability and is subject to highly variable atmospheric and oceanic forcings (Connolly and Lentz, 2014;Chen et al., 2016). The Summer MHWs were the shortest but most intense. The Spring events were similar to the Summer events, and the Autumn events a mix between Spring and Winter.
Of the 291 MHWs detected across the regions in this study, 251 had an onset phase at least 3 days in length, and 241 MHWs had a decline phase of at least 3 days. This threshold is necessary to allow for the calculation of the following measures of agreement.

Measures of Agreement
Magnitude Analysis of the dominant drivers during the onset and decline of MHWs reveals remarkable contrast. It was found that more of the change in SSTa during the onset of MHWs could be attributed to T Qnet than for the decline (Figure 2). As SSTa became larger during the onset of MHWs, T Qnet often matched or exceeded the magnitude of increase in SSTa, suggesting that the anomalous airsea fluxes drove the onset. This scenario is less common during declines (Figure 2A). The median proportion of change in SSTa during the onset of MHWs attributable to T Qnet was 0.46, and for the decline of events this was −0.04 ( Figure 2B). This means TABLE 1 | The central tendency for the duration (D, days), mean intensity (i mean ; • C), maximum intensity (i max ; • C), and cumulative intensity (i cum ; • C days) for MHWs for the total study area, by region (for all seasons), and by season (for all regions). All values shown within each grid cell are the mean ± the standard deviation (SD) for the variable in that column. The region abbreviations are given in Figure 1. Note that the large variance in the duration of events is due to the presence of a few very long events.

Group
that overall air-sea heat flux plays an important role for the intensification of MHWs, but does not consistently contribute to the decline of the temperature anomalies. In just under half of the declines of the recorded MHWs (110 out of 243), anomalous air-sea flux continues to be positive (i.e., T Qnet > 0) while the actual temperature is decreasing ( SSTa < 0; Figure 2B). In such cases, oceanic processes (e.g., horizontal heat divergence) are responsible for the decline of the temperature anomalies. This is much less common during the onset of events where the vast majority have at least some increase in temperature attributable to positive changes in T Qnet (207 out of 251; Figure 2B).
Overall, 49% (122 out of 251) of the past MHWs on the Northwest Atlantic shelf had T Qnet contribute more than half of the increase in SSTa during their onset, but only 17% (42 out of 243) of MHWs had more than half of their declines in SSTa contributed to by T Qnet (Table 2). Again, this means that air-sea heat flux is much more important for driving MHWs than it is for their decay.
The partition of the primary drivers of MHW during onset vs. decline was similar throughout all of the regions (±9%) with the exception of the Scotian Shelf. The region with the highest percent of MHWs with an onset attributable to T Qnet was the Newfoundland Shelf (57%), and the lowest was the Scotian Shelf (36%). The region with the greatest percent of MHWs whose decline was attributed to air-sea heat flux was the Scotian Shelf (32%), and the lowest was the Gulf of St. Lawrence (10%).
For seasons, the range of MHW onsets and declines attributable to T Qnet was greater than for the regions (±12%). The season with the highest percent of MHWs attributable to T Qnet was Autumn (61%), and the lowest was Summer (40%). The season with the highest percent of MHWs declining due to T Qnet was a tie between Summer and Winter (20%), the season with the lowest was Spring (7%).
Knowing when the onset or decline of MHWs may be attributed primarily to air-sea heat flux provides us with the broad overview necessary to conclude the importance of this process, but we would also like to know the relative importance of the turbulent and radiative flux terms within T Qnet during the evolution of SSTa for each event. In the following subsection we examine the RMSE of T Qx and SSTa during the onset and decline of MHWs for which T Qnet was determined to be the primary driver.

Root Mean Square Error (RMSE)
For all seasons the median RMSE for each T Qx variable during the onset of MHWs is lower than for the decline, with the exception of T Qlw in Spring (Figure 3A). This is generally true for the regions as well, except for the Gulf of Maine (GM) and the Newfoundland Shelf (NFS) where the RMSE for declines are lower than for onsets ( Figure 3B). This means that the T Qx variables tend to relate more closely during the intensification of a temperature anomaly and less so during their decay. Note that we are only looking at MHWs that are primarily driven by T Qnet during onset and/or decline, yet there is a closer agreement in the daily development of SSTa during the onset of events. This means that not only are a much higher percentage of the onset of MHWs driven by T Qnet (see section "Magnitude" and Table 2), the higher frequency variability of the SSTa from start to peak is also better described by the air-sea heat flux terms. This is further support for our finding that air-sea fluxes are more responsible for the onset of MHWs than for their decline, which is likely influenced by ocean advection or mixing.
The range of RMSE values per season are lowest for onset during Winter, then Autumn, then Summer ( Figure 3A). For the decline phase this is reversed, with the lowest values on average found in Summer, then Autumn, then Winter. The RMSE values during spring for both onset and decline are much higher than for the other seasons. This implies that the formation of warm The columns show the grouping of the MHWs and then the percentage of the onset or decline of MHWs where the majority of the corresponding change in SSTa could be attributed to air-sea heat flux. Note that these statistics are only calculated on MHW phases (onset/decline) of 3 days in length or greater. The region abbreviations are given in Figure 1.
temperature anomalies during the colder months is much more reliant on anomalous T Qx and that the daily developments in airsea heat flux more closely match to the daily changes in SSTa. For Spring, however, this is not the case. Even though these RMSE values are limited to only MHWs that were shown to be driven primarily by T Qnet , the RMSE values for the decline of Spring events are notably larger. This implies higher variability in either SSTa or T Qx or both. There are also large differences in the RMSE values for the onset and decline of events in the different regions ( Figure 3B). Most notable are the very large RMSE values for the decline phase of events in the MAB and Cabot Straight (CBS), where both regions are subject to large advective heat fluxes due to a combination of strong horizontal temperature gradient and proximity to notable ocean currents, i.e., Gulf Stream and coupled Labrador current and St. Lawrence runoff. For the MHWs driven primarily by T Qnet , T Qlh most frequently had the lowest RMSE, and T Qsw the second lowest ( Table 3, rows 1-2). T Qlw least frequently had the lowest RMSE during the onset of events across all regions and seasons, and for the decline of events it was T Qsh . This means that for the MHWs driven primary by heat flux anomalies (determined in section "Magnitude"), latent heat flux anomaly most frequently stand out as the major process modulating the SSTa during both onset and decline. In contrast, sensible heat flux and long-wave radiation have less effect in the fluctuation of SSTa during MHWs. This general pattern of lowest and highest RMSE variables applies well across the regions and seasons as well. The warmer seasons more closely follow the overall pattern above that the two variables with the lowest RMSE during onset were T Qlh and T Qsw . The decline of events in the warmer seasons were also generally best matched to T Qlh , but the second best match for both seasons was T Qlw . Note, however, that there were only three MHWs in Spring whose declines were driven primarily by air-sea heat flux anomalies.  The columns show grouping of the MHWs (total, region, and season), followed by the time series phase (onset or decline) being shown in that row. The total count of MHWs is then shown followed by a column for the count of each T Qx variable that had the lowest RMSE when compared against SSTa. The highest count per row is shown in bold italics, with the second highest in bold. In the case of ties, no more than two values are highlighted. The percentage of the total count per row for each grid cell is shown in brackets for convenience. The region abbreviations are given in Figure 1.
The colder seasons of Autumn and Winter differed more in their drivers of onset and decline. The variable with the lowest RMSE during the onset and decline of Autumn events was T Qlh , with the second lowest RMSE for onset being T Qsh . Only one MHW in Autumn didn't have T Qlh as the most relevant variable during decline. The variables with the lowest RMSE for the onset of Winter MHWs were T Qsh and T Qsw , with T Qlh and T Qsw tied for having the lowest RMSE during decline. This shows that, unlike the broad similarity in forcing mechanisms across regions, there is a notable seasonal variation in dominant air-sea flux terms, with a more consistent importance of T Qsh for onsets and declines in the colder months.

Correlation
The oceanic and atmospheric states that lead to favorable conditions for surface warming due to air-sea heat flux can be varied. Of the many ocean-atmosphere variables tested as potential drivers for the onset and decline of MHWs, several stood out. The three atmospheric variables that most strongly correlated to SSTa were: T air (positive), P-E cum (negative), and MSLP cum (positive + negative). W 10cum is also worth mentioning for the strong negative relationship it had during the onset of MHWs. The top three oceanic variables were: MLD cum (negative), SSS (negative), and T b (positive). This was determined by examining the mean correlation values for each variable, and how non-normal their distributions were. If there were little meaningful relationship between a variable and SSTa we would expect the distribution of the resultant r x values to look like a normal distribution. For these top variables that is not the case (Figure 4).
Of the atmospheric variables, T air had on average the strongest (positive) correlation during the onset of MHWs (Figure 4A). Out of both the atmospheric and oceanic variables, MLD cum had on average the strongest (negative) correlation with the onset of MHWs ( Figure 4B). This implies that of all of the non-heat flux variables that one could use to predict MHW onset, MLD cum may be the most important. For the decline of MHWs, the two most strongly correlated variables are again T air and MLD cum . One important result to note here, however, is that the negative relationship MLD cum has with the onset of MHWs means that a shoaling of the mixed layer is likely driving the MHW, but the positive relationship MLD cum has with the decline of MHWs means that the mixed layer is continuing to remain shallow throughout the decay of these events. This is yet another point of support for the conclusion that oceanic processes are more responsible for the decline of MHWs.
There were differences in the strength of these correlations between the different regions and seasons, with the differences being more pronounced across seasons. Generally speaking, MHWs that occurred in Spring and Summer showed similar means and interquartile ranges to each other for the correlations with the top air/sea variables (Figure 5). Autumn and Winter MHWs tended to have similar interquartile ranges of correlations, too, but Autumn events could also occasionally be similar to Spring/Summer events (Figure 5). One important difference is the role of MLD cum in the onset and decline of Winter MHWs. This variable is consistently strongly correlated to MHW onset and decline, except in Winter when the distribution is much more normal (centered around 0), meaning there is no clear relationship (Figure 5B; MLD cum ). This is because the water column is generally well mixed in most regions during winter and the mixed layer depth is large with little variation. Similarly, many Autumn MHWs also had little to no correlation with MLD cum during their decline.

Self-Organizing Map (SOM)
The 291 MHWs across all regions were clustered into a 12 node (4 × 3) SOM, with each node displaying different ocean and atmospheric patterns associated with MHWs. Even though the region and season of occurrence of each MHW was not provided to the SOM, it still clustered MHWs in noticeable patterns expressed in these attributes. Of the many variables that were fed to the SOM, it was the atmospheric variables that came through with the clearest patterns.
Before addressing the output of the SOM it is important to note the differences in the MHWs clustered into each node. The central two nodes (F, G) had the highest counts of all clustered MHWs and these events also tended to have the longest lasting, most intense events from the study ( Table 4). The MHWs from  The whiskers show all other values falling within 1.5 this interquartile distance, and the black dots are outliers. The further from the red line the center of the boxplots are, the less normal, and therefore more meaningful the relationship with that variable and MHW onset or decline is. The variable abbreviations may be found in Sections "Data" and "Agreement Metrics." node J were also amongst the longest lasting and most intense. The shortest, least intense MHWs were found in three of the four corner nodes: D, I, and L.
The methodology in this study was intentionally structured with MHWs occurring simultaneously across regions so that we could investigate the similarity of their synoptic scale forcing with the SOM. Of the 291 MHWs in this study, there were 321 pairs of events that overlapped by at least 1 day or more. Note that one MHW could overlap with multiple other events. Of these 321 overlapping pairs, 190 of them were clustered into the same Nodes. Of the 131 pairs that were not clustered together, 74 of them had at least one short (<10 days) event that didn't overlap much with its paired MHW. Indeed, the median proportion of overlapping days for events not clustered together was 0.21, and 0.47 for clustered events. This shows us that when synoptic scale patterns are causing multiple MHWs simultaneously across regions the SOM clusters these together accordingly, with the exception of when one of the events is very short, or the events that don't overlap by many days. For the shorter events this is likely because the cause of these events is a more transient driver that isn't resolved as clearly in the mean synoptic state for the longer event. When the events don't overlap by many days they will generally be related to different drivers that are emerging in close temporal proximity, but are not the same phenomenon and so are not clustered together in the SOM.
In the following sub-sections we will focus on the regional/seasonal and atmospheric patterns from the SOM as these provide the most clear and meaningful results. The figures showing the patterns for the SST + surface currents (Supplementary Figure 1) and Q net + MSLP (Supplementary  Figure 2) can be found in the supplementary, as well as a figure showing the date of occurrence, duration, and maximum intensity of the MHWs clustered into each node (Supplementary Figure 3).

Regions and Seasons
The central nodes tended to have more MHWs clustered into them (n ≥ 33) than the outer nodes (n = 9-32; Figure 6). The nodes in the top right corner (Figure 6; C, D, H) are filled mostly with events centered around the MAB and only a couple of events on the NFS. The top left corner nodes (A, B, E) have mostly NFS and CBS events. The distribution of MHWs is less focused within the same single region throughout the center two columns of nodes, though some nodes do show a tendency toward regionality. There is no apparent regional pattern as one moves from the top to the bottom.
For seasonal patterns, we see that the MHWs occurring in the middle two nodes (F, G) and the center-right column (C, G, K) are primarily summer events (Figure 6). The nodes in the left column (A, E, I) are predominantly Autumn events with high proportions in the neighboring seasons. The nodes in the right column (D, H, L) are predominantly Autumn and Winter events, with some occurring in Spring as well. The bottom right node (L) has a very high proportion of events in Winter and none in Summer. 10.0 ± 7.0 1.9 ± 0.3 2.3 ± 0.4 18.6 ± 11.0 L 9 6.9 ± 3.9 1.0 ± 0.3 1.2 ± 0.4 7.1 ± 4.3 The columns show the node, the count of events therein, the duration (D, days), mean intensity (i mean ; • C), maximum intensity (i max ; • C), and cumulative intensity (i cum ; • C days).

Atmospheric Variables
For air temperature anomalies we see that the air overlying the ocean in the east of the study region gets warmer as we move left and down through the nodes (Figure 7). The air over the continent becomes warmer as we move to the top right corner. The air temperatures for the middle top four nodes (B, C, F, G) are moderately warm over nearly the full study area. The air over the Atlantic Ocean tends to become warmer as we move up the nodes. The bottom row shows stronger meridional (North-South) temperature gradients whereas the top row tends toward stronger zonal (East-West) temperature gradients. The nodes in the bottom right (G, H, K, L) generally show northerly air flow whereas the top row and left column are mostly southerly air flow. If we compare these regions of high air temperature anomalies (Figure 7) to the regions of highest MHW occurrence (Figure 6), we see that there is a strong relationship between the two. This supports the earlier observation that the best atmospheric predictor of MHWs is T air (Figure 4). This is also why we see either strong negative or positive correlations for MSLP (Figure 4), FIGURE 7 | The average synoptic states for anomalous air temperature, wind speed/direction (vectors), and MSLP (mean sea-level pressure) for each node in the SOM results (A-L). Solid contours show positive MSLP anomaly and dashed contours show negative MSLP anomaly. The background colors of red to blue indicate air temperature anomaly. Note that MSLP was not one of the variables fed to the SOM, but is shown here to better outline the patterns in the wind and air temperatures.
because both low-pressure and high-pressure systems may be present during MHWs.
In most of the top row and right column nodes (Figures 7B-D,H) we see high-pressure blocking anticyclonic wind anomalies around thermal dipoles of varying magnitude with the northward moving air being anomalously warm; the warmer (colder) the anomaly the faster north (south) the anomalous wind movement is. We tend to see the opposite (low pressure systems) in the left column and bottom row (A, E, I, J, K, L). The air temperature anomalies for these low pressure cyclones tend to be less than for the high pressure anticyclones, but the pressure gradients and wind speeds are much stronger. The center nodes (F, G) are relatively quiescent states, i.e., the "dog days of summer." Nodes I and L very prominently display similar Autumn/Winter low pressure systems (Figures 6, 7). Node I shows this system linked with MHWs as it nears the Scotian Shelf, and node L shows the system not linked with MHWs until it has nearly passed from the study area. While it may appear that the Node I pattern should precede the Node L pattern in time, the average difference in days between events in these nodes is nearly 2 years. There is only one event in Node L that occurs within 2 weeks of an event in Node I.
A common pattern that emerges between atmospheric states and the onset/decline of MHWs driven by T Qnet is that any system that bring warm air up from the South onto the coast are more likely to drive the onset of MHWs due to heat flux than any other pattern. This same pattern is likewise the least likely to have heat flux driving the decline of the events that they have created.
In these results this is predominantly caused by high pressure blocking systems accompanied by anti-cyclonic wind anomalies, but some low-pressure systems also cause southerly winds. Zonal systems that lead to relatively quiescent states have a much lower likelihood of being associated with the driving of the onset of MHWs by T Qnet , but a much higher likelihood that the decline of events will be due to T Qnet . The MHWs in these nodes are also, on average, the longest and most intense.

DISCUSSION
In this paper we investigated the relationship between air-sea heat flux and the onset and/or decline of MHWs. This was done through a four step methodology: (1) finding the proportion of change in SSTa attributable to T Qnet , (2) finding the RMSE of SSTa and T Qx for MHWs whose onset and/or decline was driven predominantly by T Qnet , (3) correlations between SSTa and ocean-atmosphere state variables, and (4) a SOM analysis of the mean synoptic states during MHWs. It was found that air-sea heat flux is responsible for the onset of roughly half of the MHWs in the chosen study area, but less than a fifth of the declines. It was also found that T Qlh and T Qsw matched much more closely to daily changes in SSTa and that these differences were not the same across regions/seasons. Of the many ocean-atmosphere state variables investigated, T air and MLD cum were found to have the strongest relationships with MHW onsets and declines; with MLD cum decreasing during MHW onset and T air decreasing during MHW decline. Lastly, the SOM analysis revealed that both low-pressure systems and high-pressure blocking cells could be present during the occurrence of MHWs.

Measures of Agreement
The relationship between T Qnet and SSTa has lead us to infer that MHW onset in the NW Atlantic is predominantly forced by air-sea heat flux, while MHW decay is likely due to ocean advective or mixing processes not captured in our analysis. Not only is air-sea heat flux driving the onset of nearly half of the MHWs in the NW Atlantic and less than one fifth of the declines, when MHW decline may be attributed to air-sea heat flux, the fit of the T Qx variables to declining SSTa is worse. This is more support for our inference that the decline of MHWs is more often related to oceanic processes, e.g., advective flux, which is not explicitly captured in our methodology. The variable with the best correlation with SSTa during the decline of a MHW is T air . This alone is likely not what is driving the SSTa anomaly, but is rather a secondary variable being controlled by the same advective/mixing process responsible for the decline in the SSTa, or is itself being driven by the changes in SST.
It was also seen that MLD cum tended to have a strong positive correlation with SSTa during the decline of the events, implying a continued shoaling of mixed layer during the decline phase (a negative correlation would mean that the MLD is deepening as the temperature anomaly went back down to zero). This means that most MHWs begin when the MLD is shoaling, or has already shoaled, but then the MLD does not begin to deepen again at the decline of the event. Advective or mixing processes must thus be responsible for the decline of the temperature anomaly. It may also be that the T Qx variables are not found to relate closely to SSTa at the decline of events specifically because of the lack of MLD agreement during the decline phase. This continued shoaling of the MLD during the end of MHWs implies that a shoaled MLD may be a more long-lasting features than the MHWs, and may actually be the driving force behind the occurrence of two or more short events in rapid succession. This observation was tested, but it was found that short events (<10 days duration) occurring within 2 weeks or less of each other were almost always accompanied by an anomalously deep MLD, not shallow. The columns show what the dominant atmospheric pattern, region, and season of occurrence is for each SOM node. A brief characterization of the MHWs is then provided. The last two columns give an indication to the amount of MHWs whose onsets or declines were driven by T Qnet relative to the overall averages of 49% for onset and 17% for decline. We use the definitions as follows: same = overall average ± 3%, more ≥overall average + 3%, much more ≥overall average + 10%, less ≤overall average − 3%, much less ≤overall average − 10%, none = 0%. To aid in the concise summary of the results a shorthand notation is used. For patterns: LPCS, low pressure cyclonic system; HPB1, high pressure blocking with anticyclonic wind anomalies; HPB2, high pressure blocking with little to no anomalous wind. Each node tends to be focused around one or two seasons and regions. The region abbreviations are given in Figure 1.

Self-Organizing Map (SOM)
We saw from the SOM that there were three primary types of patterns.
(1) Low pressure cyclonic systems that occurred mostly in Autumn and Winter.
(2) High pressure blocking systems with anticyclonic wind anomalies that occurred mostly in the Spring and Summer; and (3) high pressure blocking cells that would sit over the study area for long periods of time almost exclusively in Summer. A brief synopsis of the dominant patterns of each node has been tabularized for convenience ( Table 5).
It was also noted that when the onset of a MHWs was driven by T Qnet , it was much more likely that the predominant atmospheric pattern present was bringing warm southerly wind against the coast. This could be due to low pressure systems, but was more often due to high pressure blocking. This same pattern was also very infrequently related to MHWs whose decline was driven primarily by T Qnet . Note also that the quiescent atmospheric states associated with long lasting and intense MHWs were also much more likely to occur during MHWs whose declines were driven by T Qnet . This provides two very strong starting points for future research on the relationships between atmospheric drivers of air-sea heat flux and the formation of MHWs.

Methodology
The novel methodology created for this study has proven effective at quantifying and numerating the processes responsible for the onset and decline of MHWs, but there are some issues that demand further consideration. One such issue is how best to arrive at the definition of a single MHW when multiple smaller events occur within a short period of time. For consistency throughout this study we adhered to the method for MHW detection from Hobday et al. (2016), which states that events are counted as separate once the daily SST value drops below the 90th percentile threshold for more than 2 days. We acknowledge that this definition for gaps between events could be debated further, but here we stick to the definition to be consistent with earlier studies. Large and intense events may occur over broad spatiotemporal scales and may have long-lasting impacts on the background state. As a result, a large and intense event may be detected as multiple MHWs in different regions and/or as multiple MHWs close to each other in time. Even so, the onset and decay of these individual MHWs across space and time may not necessarily be driven by the same processes. In this study we decided to keep the smaller events separate from each other in order to provide the measures of agreement between datasets as rich a diversity of 'types' of MHWs as possible. Along this same avenue of thought we opted out of performing any additional smoothing on the SST time series. This could have potentially worked to join short rapidly occurring MHWs together without compromising on the pre-established methodology, but that in turn would have exaggerated the importance of these minor events by combining them into much longer single events. The effect of low-pass smoothing on MHW detection would also make for an interesting future study, though to some extent it is already known that the further smoothing of the seasonal climatology used for MHW detection leads to the detection of more intense events in Summer/Winter and less intense events in Spring/Autumn (Schlegel et al., 2019). Finally, it is known that preconditioning is an important consideration for the study of anomalous temperatures in the NW Atlantic (Chen et al., 2016). This methodology does not account for preconditioning specifically as this falls outside of its scope. Given the known importance of this phenomenon in the area it would be an excellent topic to investigate in more depth in a future study. As we have shown here with a follow up investigation into the potential for a shoaled MLD to function as preconditioning for short rapid MHWs, the relationship of this phenomenon with MHWs is not as clear as would be hypothesized and deserves further study.
A final and interesting consideration that was not accounted for in this methodology was the effect that anthropogenic warming will have on the SOM analysis. Consider that as the world is gradually warming, it becomes relatively easier for an atmospheric or oceanic driver to force a MHW because the SST values needed in order to be above a historic climatology are less than they would be in the cooler past. This means that a more intense synoptic scale pattern would need to be present to drive MHWs in the past than in the future. However, most of the SOM nodes contain MHWs that occur over the whole time range of 1993 -2018, so any possible warming signal appears not to have had a noticeable impact on the results. It would, however, be an interesting study to look at results from much longer (centennialscale) time series to see at what point anthropogenic warming has a significant impact on the types of synoptic scale forcing that can be related to the occurrence of MHWs.

CONCLUSION
In this study we used a suite of agreement metrics to develop an overall understanding of the drivers of the onset and decline of MHWs over the NW Atlantic continental shelf. It was found that nearly half (49%) of the onset of MHWs could be attributed to air-sea heat flux, but less than a fifth (17%) of the decline of MHWs could be attributed to the same process. Of the individual T Qx variables, T Qlh and T Qsw were most frequently the top two relevant drivers of the onset and decline of MHWs across all regions and seasons. T Qlw was on average the least important (highest RMSE value) for onsets, and T Qsh for declines across all regions and seasons. Over the entire study domain, T Qlh and T Qsw were still the most important drivers for onset in all seasons, but for the colder seasons of Winter and Autumn the onset of events often corresponded better to T Qsh . Anomalous air-sea heat fluxes are responsible for the onset and decline of more MHWs during the cold seasons (on average), and anomalous advection may be more responsible for MHWs during the warm seasons. This is presumably associated with the seasonal variability of advective fluxes, and is consistent with the results in Chen and Kwon (2018).
Air-sea heat flux alone may not be the only potential driver, or indicator of the development of MHWs. In order to have a better understanding of the driving processes, we further evaluated the correlations between SSTa and a suite of atmospheric and oceanic state variables. It was found that the top state variables to correlate with SSTa were T air and MLD cum . This shows that not only does the heat flux itself (i.e., Q x ) play a role, but importantly MLD shoaling is critical for the onset of most events. In comparison, the decline of MHWs was strongly positively correlated to T air . It must be noted, however, that the decline of most MHWs are not explainable via air-sea heat flux, which implies that the strong relationship that T air has with SSTa during the decline of MHWs should not be over-interpreted. It appears that persistent atmospheric forcing is necessary to maintain MHWs, so that when this forcing dissipates other processes (presumably oceanic processes) can take over and dissipate the accumulated heat in the mixed-layer associated with the MHW. The answer to what drives the decline of MHWs lies in a wider, more comprehensive future analysis that includes ocean advection and mixing.
Lastly, in order to identify common synoptic scale patterns during MHWs a SOM analysis was used. From this we determined that there were three primary atmospheric patterns during MHWs. The first was a high pressure blocking cell present over the region accompanied by anticyclonic wind anomalies, usually in the warmer months, causing anomalously high levels of heat flux through the ocean surface. These patterns were most (least) often associated with MHW onsets (declines) driven by T Qnet . The second pattern was a low-pressure cyclonic system present over the region, usually in the colder months, that transported heat into the region from elsewhere and provided strong winds to enhance air-sea heat fluxes. The final pattern was found most often in the warmer months, and was associated with the longest and most intense MHWs. This is the "dog days of summer" type of synoptic pattern in which the wind, air, and MSLP exhibit weakly anomalous behavior, but they persist for much longer than they normally would and create a slow and steady build up of heat in the surface of the ocean. These patterns were most frequently associated with MHW declines driven by T Qnet .

DATA AVAILABILITY STATEMENT
The datasets generated for this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://github.com/ robwschlegel/MHWflux.

AUTHOR CONTRIBUTIONS
RS led the study and prepared the majority of the text, figures, synthesized the comments, and uploaded the manuscript. EO conceptualized the SOM portion of the study. KC conceptualized the magnitude and RMSE portions of the study. EO and KC both provided guidance on the physical oceanographic aspects of the study. RS conceptualized the correlation portion of the study and introduced the idea to split the analyses into onset/decline phases. EO and KC provided several rounds of comments on the manuscript as it was developed. All authors provided responses to reviewers during the editorial process.