Combined oceanic and atmospheric forcing of the 2013/14 marine heatwave in the northeast Pacific

An unprecedented warm sea surface temperature (SST) anomaly event, namely, the Blob, occurred in the northeast Pacific during the winter (October–January) of 2013/2014, causing substantial economic and ecological impacts. Here, we explore the driving forces of the Blob from both atmospheric and oceanic perspectives and show that the Blob primarily resulted from weak wintertime cooling due to the reduced air-sea heat flux transfer from the ocean to the atmosphere and the reduced horizontal advection of cold water in the upper ocean. Both mechanisms were attributed to an anomalous high-pressure system over the study region. Specifically, the anomalous air-sea heat flux, which was dominated by turbulent heat flux anomalies, was mainly induced by the increased air temperature (i.e., with a contribution of approximately 70%) and the weakened wind speed associated with the high-pressure system. The reduced horizontal heat advection was mainly due to the weakened winds acting on the ocean temperature meridional gradient. Using a regional ocean numerical model with different experimental runs, we evaluated the contributions of air temperature and wind drivers to the Blob at both the surface and subsurface of the ocean. The Blob was absent when the model was forced by climatology-air-temperature. Both the SST and integrated ocean heat content (OHC, 0–150 m) decreased, and the mixed layer depth (MLD) was deeper than that in the control run forced by real atmospheric conditions. In the climatology-winds experiment, obvious warm anomalies still existed, which were similar to but weaker than the control run. The SST (OHC) and MLD values in the climatology-winds run were between those of the climatology-air-temperature run and the control run. Compared to former studies that attribute the formation of the Blob to an anomalous air-sea heat flux and horizontal advection mainly induced by reduced winds, our study demonstrates that anomalous warm air temperatures played a more important role in its formation.


INTRODUCTION
Extreme sea surface temperature (SST) warming, namely, marine heatwaves (MHWs), defined as discrete, prolonged anomalous warm-water events 1,2 , have become more frequent and intense in recent years 3 . MHW events have been recorded in the global oceans, including the Mediterranean Sea 4,5 , western Australia 6 , northwestern Atlantic 7,8 , Tasman Sea 9 , and northeastern (NE) Pacific 10 . In particular, a warm anomaly developed in the NE Pacific from October 2013, which is referred to as "the Blob" 10 . The Blob first appeared in the southern Gulf of Alaska during the boreal winter of 2013/2014 (October-January), spread along the west coast of North America and subsequently stretched south to Baja California 11 . Accompanying the warm ocean temperature anomalies, a delayed onset of upwelling was observed 12 , and a record-breaking harmful diatom bloom 13,14 occurred off the Californian coast, which persisted from May to October and resulted in the closure of several economically important fisheries 13 .
Multiple processes have been suggested to contribute to ocean temperature variations, including air-sea heat flux, horizontal advection, and vertical processes 15 . Specifically, air-sea heat fluxes, composed of shortwave and longwave radiative fluxes as well as latent and sensible heat fluxes, are important links between the ocean and atmosphere 16,17 . In the winter of 2013/2014, there existed an anomalous high-pressure system that weakened the background westerlies in the NE Pacific due to the weakened Aleutian Low 10,18 . The strongly positive sea level pressure (SLP) was responsible for the positive temperature anomalies in the NE Pacific, which suppressed the heat losses from the ocean to the atmosphere and induced less cold advection than usual in the upper ocean 10 . The anomalous heat transferred by horizontal advection was partly due to the reduced wind-forced (Ekman) transport from the north acting on the meridional upper ocean temperature gradient, and an additional contribution was caused by the eastward component of the background current acting on the upper ocean temperature anomalies with zonal gradients. Additionally, the weakened winds related to the high-pressure system also affected the vertical mixing that inhibited the cold subsurface water from entraining to the surface 10 . Reduced cooling during that winter resulted in more heat being concentrated in shallower mixed layer, resulting in a large SST anomaly that exceeded those of historical records from satellite observations and reanalysis datasets 19 .
In general, the generation and persistence of MHWs in the NE Pacific may also be modulated by large-scale modes of climate variability 15 . A list of such modes includes at least the North Pacific Oscillation (NPO) 20 , the tropical Northern Hemisphere (TNH) pattern 21 and the El Niño-Southern Oscillation (ENSO) 22 . Among them, the NPO and TNH patterns are responsible for the generation of MHWs, including the Blob 18,23 . After the atmospheric ridge forced the warm Blob, a neutral to weak El Niño pattern developed in late winter 19 . Even though the El Niño phenomenon during 2014/2015 winter was not yet fully developed, equatorial warming was sufficient to produce atmospheric teleconnections to the extratropics 19 . The atmospheric forcing was changed, which induced warming in the previous winter and largely drove the Aleutian Low southward 24 . After the Aleutian Low changed, the Blob evolved from the Gulf of Alaska warming pattern to the arc-shaped warming pattern, similar to the oceanic expression of the Pacific Decadal Oscillation (PDO) 19,24,25 .
Most MHWs have been investigated using SST products, which are available as global-scale, continuous, long-term records 10,19,26 . However, MHWs are not limited to the surface layer, and anomalous warming also occurs at the subsurface 6,27-29 . The unusually warm water of the Blob extended from the surface to an 80 m depth in early 2014 and then reached a depth of 140 m in the spring (March-May), which was revealed from Argo float data in the NE Pacific 11,27 . Surface satellite data indicated that the Blob disappeared in September 2016, while Jackson et al. 27 revealed that anomalous warm conditions continued to exist below the mixed layer through at least March 2018 in the coastal ocean by using Argo float and ship-based data. The subsurface warming related to the MHW could last longer than the surface warming and even reemerge in subsequent winters, affecting the mixed layer temperature 30 , as observed in the NE Pacific in 2016 27 . The persistence of subsurface warming and the possible reoccurrence of surface anomalies might be conducive to the occurrence of multiyear MHW events 28 , which have a significant impact on marine ecosystems, especially by affecting marine life in the subsurface layer 13 .
Previous studies have revealed that this MHW event was partly induced by the anomalous air-sea heat flux, which was dominated by the turbulent heat flux anomaly attributed partly to the wind speeds 10,19 because the turbulent heat flux was determined according to the wind speed and air-sea temperature/specific humidity difference 31 . Similar air-sea heat flux-type MHW events have been experienced worldwide, revealing that air temperature also plays an important role in the formation of MHWs. For example, in the boreal summer of 2003, the anomalous high atmospheric pressure related to almost no wind and extremely high surface air temperature 32 led to a significant warming event over the western European continent and the western Mediterranean Sea 32,33 . Olita et al. 4 revealed that for a similar air-sea heat flux-type MHW, the reduced heat losses from the ocean to the atmosphere were associated with decreased wind and a substantial increase in air temperature 5 . Manta et al. 34 found that an unprecedented combination of persistent low wind speed and an extremely high air temperature was probably the reason for the 2017 record MHW event on the Southwestern Atlantic shelf. These wind and air temperature anomalies induced significant changes in the air-sea heat flux, especially the turbulent heat flux. However, the investigation of the relative contributions of wind speed and air temperature to the formation of the Blob is insufficient. In this study, we applied a regional numerical ocean model and reanalysis data to describe the spatiotemporal patterns of anomalous environmental conditions in the NE Pacific and investigated their connection to the Blob, specifically to facilitate assessments of the various components of the air-sea processes, which are key to surface and subsurface warming. Understanding the key factors affecting the formation and structure of this MHW may provide insight into the prediction of MHWs in the future.

RESULTS
Observations and driving mechanisms of the Blob The SST anomalies were investigated by using the latest reanalysis dataset (ERA5) from the European Center for Medium-Range Weather Forecasts (ECMWF) (Fig. 1a). Prominent warm anomalies were observed in February 2014 in the NE Pacific, and the study region was designated as the area between 135°W-150°W, 40°N-50°N (Fig. 1a). The average SSTA within the study region reached approximately 2.0°C, and the maximum of the temperature anomalies in the study area exceeded 2.5°C. The corresponding spatially averaged time series (Fig. 1b) showed that the anomalous warming initiated during late 2013 (October-November), developed in early 2014 (January-February), and disappeared in late 2014 (October). The SSTA calculated from other reanalysis datasets, e.g., the National Center for Environmental Prediction (NCEP) and the Objectively Analyzed air-sea Fluxes (OAFlux), also captured similar patterns of anomalously high SSTs ( Supplementary Fig. 1).
Atmospheric forcing played a key role in driving the SST variability and was analyzed to explore the driving mechanisms of the Blob. The SLP anomalies over the region showed significantly positive values between 40°N and 60°N, with a center near 50°N (Fig. 2a). Compared with the corresponding climatology calculated from 1982 to 2019, these positive SLP anomalies tended to weaken the Aleutian Low and resulted in reduced surface winds. Indeed, the wind speed in the winter of 2013/2014 was the weakest of the past 38 years (Fig. 2c). Consistently, the windforced (Ekman) transport of cold water from high latitudes acting on the meridional ocean temperature gradient also decreased ( Supplementary Fig. 2a). The combined positive SLP anomalies and negative wind speed anomalies suggested that atmospheric forcing was an important driver of the Blob. During the winter with the anomalously warm SSTs, there was also an anomalously southeastward wind over the study area on the west side of the anticyclonic SLP anomalies, which transported a warm air mass from the low-latitude region. Therefore, the time series of the air temperature (  (Fig. 2c). In particular, the monthly averaged air temperature in January 2014 was 2.7°C warmer than that in the same period from 1982 to 2019, a phenomenon that is mainly attributed to the existence of anomalous SLP conditions (Fig. 2b). The daily time series revealed that the warm air temperature anomalies occurred 5 days earlier than the SST anomalies (not shown), which indeed indicated the driving role of atmospheric forcing.
The physical processes responsible for the Blob were studied through an analysis of the mixed layer temperature budget, reflecting the temperature changes due to air-sea heat flux, horizontal heat advection and vertical processes 35 . The upper mixed layer temperature budget was quantified in the study region (black box in Fig. 1a) to obtain a time series of each winter from 1982 to 2019 (Fig. 3). The change in the mixed layer temperature and each physical process contributing to this change were calculated from October to January of each year with reanalysis datasets, including the NCEP Global Ocean Assimilation System (GODAS) gridded product and ERA5 datasets. During winter, the average mixed layer temperature decreased by approximately 5.4°C from October to January, 2.9°C of which was attributed to the heat transferred by air-sea heat flux and 1.5°C of which was due to cooling from vertical entrainment; the contribution of horizontal oceanic heat advection was 0.9°C. The results illustrated that the anomalous warming was significant in the winter of 2013/2014, and the mixed layer temperature was approximately 48% (2.6°C) higher than the average over the period from 1982 to 2019 (Fig. 3). The mixed layer temperature budget analysis demonstrated that the positive ocean Value of the changes in the mixed layer temperature in the region between 40°N-50°N, 135°W-150°W (black box in Fig. 1) from October to January and the mixed layer heat budget terms that caused this temperature change from 1982 to 2019. The budget terms include the heat transferred by net air-sea heat flux (red), horizontal advection (green) and vertical entrainment (gray) terms. The values represent the degree of temperature change associated with each term.
temperature anomalies in the NE Pacific were mainly dominated by the positive air-sea heat flux anomalies, contributing approximately 56% of the warm anomaly. The anomalous horizontal advection term also played a role in the formation of the mixed layer warming with a contribution of approximately 29%. The residuals, e.g., the entrainment from below, contributed 15% of the warm anomaly. The data from the Simple Ocean Data Assimilation (SODA) reanalysis global gridded product also revealed the major role of the anomalous air-sea heat flux in the formation of the Blob ( Supplementary Fig. 3).
Considering the critical role of the anomalous air-sea heat flux in generating the Blob (Fig. 3), we further examined the major components of the flux that drove heat from the atmosphere to the ocean. The monthly mean air-sea heat flux in climatology has obvious seasonal variation, which is composed of shortwave and longwave radiative fluxes and latent and sensible heat fluxes (Fig.  4a), and the positive heat flux represents the heat into the upper ocean. The climatological monthly average of the air-sea heat flux in winter, which was the average of the dashed line from October to January, was −84.2 W m −2 , indicating heat losses from the ocean to the atmosphere in winter (Fig. 4a). In the winter of 2013/ 2014, the air-sea heat flux was −42.9 W m −2 , revealing a significant positive anomaly (41.3 W m −2 ) compared with the climatological value and implying a weaker cooling during this winter (Fig. 4). A detailed investigation revealed that this decrease in heat transfer from the ocean to the atmosphere was mainly dominated by the anomalous turbulent heat flux, especially the latent heat flux (Fig. 4b). The value of the anomalous turbulent heat flux was 35.3 W m −2 during the winter of 2013/2014, which was approximately 85% of the positive air-sea heat flux anomalies. The radiation heat flux (shortwave and longwave radiation flux) contributed less to the change in the net air-sea heat flux than the turbulent heat flux (Fig. 4b).
The variability in the turbulent heat flux was determined by multiple factors, e.g., the wind speed, air-sea temperature difference and air-sea specific humidity difference 36 . Because the specific humidity strongly depends on the local air temperature, the increased air temperature generally increases with an increase in the specific humidity 4 . To simplify the interpretation of the analysis, the latent and sensible heat fluxes were considered to be mainly affected by the air temperature and wind speed. The results revealed that the positive turbulent heat flux anomalies were mainly caused by increased air temperature (Fig. 5), with a contribution of 70%, and the weakened wind speed and nonlinear effects related to both factors contributed to the rest. Thus, both factors were important for driving the anomalous turbulent heat flux, and the air temperature played a more important role in latent heat flux and sensible heat flux anomalies, contributing approximately 63% and 86%, respectively ( Supplementary Fig. 4). Consistently, the positive turbulent heat flux anomalies caused by the anomalously warm air temperatures in the NCEP and OAFlux datasets accounted for approximately 68% and 66% of the turbulent heat flux anomalies, respectively, and the anomalously warm air temperature contributed approximately 59% (84%) and 52% (94%) to the latent (sensible) heat flux anomalies, respectively.  87 89 91 93 95 97 99 01 03 05 07 09 11 13 15 17 19 Year advection. Here, winds can change the ocean circulation through dynamic effects or modify the ocean-atmosphere thermal coupling through the wind speed in the bulk formula of turbulent heat flux. The air temperature, whose anomalies are significantly correlated with the specific humidity anomalies, is primarily responsible for determining turbulent heat fluxes. In this study, we employed a Regional Ocean Model System (ROMS) to investigate their individual contributions to the formation of the Blob, and wind and air temperature were the variables of interest. First, we used a control run (CTRL) as the baseline, which was forced by daily longwave and shortwave radiation, air temperature from NCEP, and daily winds from the NOAA multiple-satellite blended dataset. Next, we conducted three experimental runs to isolate the effects of winds and air temperature on the formation of the Blob. The first experiment (Case 1) was forced by climatological daily longwave radiation, shortwave radiation and air temperature, while the wind was similar to the CTRL. Compared with the CTRL, the second (Case 2) and third experiments (Case 3) were forced by the same forcing except that the realistic air temperature and winds were switched to the climatological daily air temperature and winds, respectively. The response of the ocean to the variables of interest can be obtained by comparing the difference between the simulations with the CTRL. In the CTRL, where both the wind and air temperature were forced with real observations, the SST showed a prominent positive anomaly centered at 45°N and 140°W in February 2014 (Fig. 6a). This feature resembled the pattern shown in the reanalysis data (Fig. 1a). Similar to the CTRL, the Case 3 simulation produced a clear SST anomaly pattern, but the corresponding magnitude was approximately 30% smaller than those in the CTRL (Fig. 6e). The signals of the significantly warming SST were absent in both Case 1 (Fig. 6b) and Case 2 (Fig. 6d), and the corresponding SST time series were characterized by clear seasonal variabilities in SST (Fig. 6c). The average SST during the winter of 2013/2014 was 11.3°C in both Case 1 and Case 2, which is lower than that in the CTRL (12.1°C). Thus, the air temperature played an important role in driving the positive SST anomalies during the winter of 2013/2014. In Case 3, the SST also showed clear seasonal variability, and the SST (11.7°C) was higher than that in Case 1 and Case 2 but lower than that in the CTRL (Fig. 6c), indicating that the winds also contributed to the formation of positive SST anomalies but had a smaller (27%) role than the air temperature.

Contribution of factors in modeling analyses
MHWs are not limited to surface temperature anomalies but also penetrate to the subsurface. Thus, the vertical extent of the anomalous warming water and its evolution in the study area were fully analyzed with the integrated ocean heat content (OHC) over the upper 150 m. The model outputs captured the anomalously warm temperature profile, which was consistent with the observational measurements. In late 2013, there were strong vertical temperature anomalies that were mainly confined to depths shallower than the mixed layer. The corresponding temperature anomaly was over 2.0°C (Fig. 7a), which was consistent with the surface warming. Later, this warm anomaly started to appear in the subsurface waters beneath the mixed layer. Associated with the propagation of the temperature anomaly to deeper water, this anomaly also expanded eastward. Consistently, anomalously warm surface waters were found in coastal regions along North America 13,27,37 . In addition, during the winter of 2013/2014, the mean values of the MLD increased from 38.8 m to 42.7 m in Case 1 and to 43.4 m in Case 2, increasing by approximately 10% and 12% from the CTRL, respectively (Fig. 7b, c, Table 1). In Case 3, there was a significantly positive vertical temperature anomaly in late 2013 similar to that in the CTRL, and the MLD decreased to 27.5 m, which was equivalent to 71% of the value in the CTRL (Fig. 7e, Table 1). Moreover, the integrated OHC of the upper ocean (0-150 m) in the CTRL was 3.2 × 10 8 J m −2 higher than the mean value of the 3-year (2010-2012) average. The integrated OHC (0-150 m) decreased during the three experimental runs (Fig. 6f), with positive anomalies of 0.9 × 10 8 J m −2 in Case 1, 1.2 × 10 8 J m −2 in Case 2 and 1.6 × 10 8 J m −2 in Case 3, which were approximately 28%, 38% and 50% of the value in the CTRL, respectively (Table 1). These results indicate the importance of anomalous atmospheric conditions in driving subsurface properties in the NE Pacific during this event.

DISCUSSION
The Blob in the NE Pacific, which was experienced during the winter of 2013/2014 (Fig. 1), was investigated with reanalysis products and a regional ocean numerical model to comprehensively determine its driving mechanism. The atmospheric forcing and its connection with the Blob in the NE Pacific are summarized with a conceptual model (Fig. 8). The Blob primarily resulted from the reduced Aleutian Low that caused an anomalous highpressure system. This high-pressure system not only weakened the intensity of surface winds but also increased the air temperature. As a result, the heat losses from the ocean to the atmosphere were reduced, and the horizontal advection of cold water and wind-driven ocean vertical mixing were depressed in the NE Pacific (Fig. 3). Thus, more heat was maintained over a shallower mixed layer depth, inducing ocean surface warming anomalies that exceeded the climatological state by~2.0°C. This feature was highly consistent with the observations, e.g., the SST was 1.0°C-4.0°C higher than average in the NE Pacific 13  Often, MHWs related to anomalous air-sea heat fluxes are accompanied by increased shortwave radiative heat fluxes resulting from less cloud cover and more insolation or less turbulent heat flux due to a warm surface air mass 35 . These processes, which can occur independently or simultaneously, are usually responsible for air-sea heat flux-driven MHWs. During the winter of 2013/2014, the shortwave radiation decreased slightly compared to the normal condition (Fig. 4b), but the turbulent heat fluxes, which are the primary mechanisms for transferring the absorbed solar radiation back to the atmosphere, showed positive anomalies due to the weakened Aleutian Low (Fig. 2) 18 . The positive anomalies of the turbulent heat flux were consistent with the mixed layer heat budget analysis result, which revealed that the anomalous air-sea heat flux played an important role in the Blob's formation (Fig. 3). Other MHW events caused by anomalous air-sea heat flux have been observed in other parts of the global ocean. For example, in the summer of 2003, a near lack of wind and a high air temperature associated with an anomalous highpressure system resulted in a significant MHW in the western Mediterranean Sea 32,33 . Another MHW in the Mediterranean Sea was associated with a significant increase in air temperature, reduction in wind intensity and decrease in all components of the air-sea heat flux from the ocean, characterized by an obvious SST warming by 3.0-4.0°C 4,5 . These events all indicated the importance of wind fields and air temperatures in driving the occurrence of anomalous ocean warming. In the winter of 2013/ 2014, the reduced wind speed associated with the high-pressure system was connected to the NPO (Fig. 2a) and suppressed the evaporation of sea surface water, releasing both water masses and energy to the atmosphere 39 . A substantial increase in air temperature was observed over the waters of interest (Fig. 2), accompanied by the obvious positive anomalies of air humidity in the winter of 2013/2014 ( Supplementary Fig. 2b), which led to changes in the air-sea temperature/humidity differences and further suppressed the heat losses from the ocean to the atmosphere. Therefore, the combination of the reduced wind speeds and increased air temperatures significantly reduced the transfer of turbulent heat flux from the ocean to the atmosphere, especially the air temperature, which played a more important role than the wind speed (Fig. 5).  Winds can not only modify ocean-atmosphere thermal coupling through wind speed in the bulk formula for turbulent heat flux but also change heat advection by ocean circulation through dynamic effects. In this event, the reduced southward cold-water transport associated with weakened background westerlies, which were due to the weakened Aleutian Low partly contributed to anomalous positive horizontal heat advection (Fig. 3). An additional contribution was made by the eastward component of the background current associated with the westerlies acting on the upper ocean temperature anomalies with a zonal gradient, which was significantly weaker than the wind-induced meridional transport (not shown). In addition, the negative wind stress curl anomalies and, therefore, the Ekman pumping anomalies might enhance ocean stratification, which was consistent with the relatively weak negative entrainment that transported less cold water from depths below the mixed layer to the surface layer ( Supplementary  Fig. 2a). Meanwhile, the weakening of the westerlies limited the energy transferred from the atmosphere to the ocean for mixing processes 40 and contributed to the relatively weak entrainment of cold water from below. The MLD in the winter of 2013/2014, induced by the weakened winds, was significantly shallower than the normal condition (decreased by 44%, Fig. 7). Thus, more heat contained in the shallow mixed layer increased the ocean stability and reduced the entrainment of cool water from below. Together with the abnormal air-sea heat flux, which was dominated by the turbulent heat flux induced by weakened wind speed and warm air temperature, there was less wintertime cooling in 2013/2014 than usual, resulting in the Blob (Fig. 4b). The increased SST further reduced the winds and established a positive ocean and atmosphere feedback, referred to as the winds-evaporation-SST feedback (WES) 41 . This positive feedback might have activated the meridional modes 42,43 , which are conducive to the development of El Niño in the central equatorial Pacific by their propagation and amplification of the warming anomalies from the subtropics to the central equatorial Pacific 24 . Once the El Niño had developed, the abnormal tropical convection excited atmospheric Rossby waves in the higher troposphere, which injected the changes into the extratropics atmospheric forcing through teleconnections during 2014/2015 winter; then, positive ocean temperature anomalies evolved from the Gulf of Alaska's warming pattern to an arcshaped warming pattern 24,44 .
The contribution of the anomalous atmospheric conditions driving the Blob, together with their impacts on the surface and subsurface water layers, was quantified with a regional ocean numerical model. By comparing the SST anomaly distribution and the values of the winter-averaged SST and integrated OHC (0-150 m) from the model simulations, we found that the abnormal winds and air temperature both played important roles in the Blob's formation, especially in terms of the increased air temperature. In the absence of abnormal winds, the intensity of the ocean temperature anomaly was weakened because the winds affected the air-sea heat flux through the wind speed in the bulk formula for the turbulent heat flux and changed the ocean heat advection through dynamic effects, consistent with a previous study 10,19 . Additionally, in this study, we emphasized the importance of the air temperature anomalies on the Blob's formation, which was the key factor in the positive turbulent heat flux anomalies during the winter of 2013/2014. Specifically, in the absence of abnormal air temperatures, there were almost no significant ocean temperature warming anomalies (Fig. 6). In addition, the CTRL model result revealed that the Blob was indeed a three-dimensional structure where ocean temperature warming persisted not only in the surface but also in the subsurface layer (Fig. 7), which was consistent with the observation that the surface and subsurface warming at 140 m were both observed roughly under the location of the Blob 27 . Previous studies indicated that the temperature anomalies from the surface to the subsurface in the NE Pacific were related to diabatic subduction 27 , detrainment, horizontal advection 11,45 and/or adiabatic isopycnal heave processes 28 . During the transition of the MLD from winter to spring in 2014, there existed significant subsurface warming below the mixed layer (Fig. 7), which was consistent with the Argo observations that indicated that diabatic subduction or horizontal mixing might be responsible for the penetration of temperature anomalies within the seasonal pycnocline 28 .
Previous studies have revealed that the onset and persistence of the Blob might have been modulated by climatic variability modes, such as the ENSO, NPO, and TNH. Among them, the NPO and TNH patterns might be responsible for the generation of the Blob 18,23 , indicating atmospheric teleconnection from tropical regions to the NE Pacific 18,46,47 . Additionally, the prolonged persistence of the Blob was partially related to the teleconnections back to the extratropics from the tropics through the atmosphere 24 , such as the weak El Niño phenomenon in the winter of 2014/2015, which was a key source of the long-lived event. In addition, the multiyear persistence of the Blob had unprecedented impacts on marine ecosystems and socioeconomically important fisheries. Accompanying the warm ocean temperature anomalies, low primary productivity was observed in January 2014 in the NE Pacific due to the reduced nutrient flux supplied to the mixed layer associated with the increased ocean stratification 48   ocean temperature significantly changed the species composition of the community 13 , and the shift in planktonic species led to the starvation and death of seabirds and marine mammals 13,49 . An unprecedented harmful algal bloom dominated by diatoms occurred from southern California to southeast Alaska in 2015, which led to fishery closures that resulted in millions of dollars in economic losses 13 . Due to the possibility of continuous ocean warming in the twenty-first century, the MHW frequency in the global ocean is expected to continue to increase, with serious implications for marine ecosystems 3,50 . Based on our results and those of other studies, it is particularly important to investigate the formation and vertical structure of MHWs and to assess the variabilities in their associated biological responses by combining subsurface observations, e.g., BGC-Argo, and coupled physicalbiogeochemical models. In this way, we can gain better insight into the potential effect of a given event and subsequently benefit fisheries and wildlife management in the ocean.
Overall, this multifaceted analysis delineates the driving forces of the Blob under anomalous atmospheric and oceanic conditions, including the changes in the SLP, air temperature, wind field and mixed layer. Our results highlight the importance of warm air temperatures in the formation of the Blob. Due to the anomalous environment, less wintertime cooling resulted from the reduced air-sea heat flux transferred from the ocean to the atmosphere and weakened the horizontal advection of cold water; thus, more heat was retained in the shallower water column and formed the Blob. Moreover, the model results confirmed that the anomalous air temperature played an important role in Blob's formation, and the Blob was indeed a three-dimensional structure where warming occurred simultaneously in the surface and subsurface layers. By improving the understanding of the formation of the Blob, this study provides a dynamical framework that may help to predict the occurrence of extreme events and mitigate the risk induced by future climate change.

Reanalysis data
We used the latest reanalysis products (ERA5) from the ECMWF from January 1982 to December 2020 and the ERA5 reanalysis for global climate and weather data 51 . This dataset includes the SST and other monthly variables, such as the 10 m surface winds (u, v, and wind speed), SLP, air temperature at 2 m, net shortwave radiation flux, net longwave radiation flux, net latent heat flux, and net sensible heat flux. The spatial resolution of the dataset is 0.25°× 0.25°, and it covers the period from 1979 to the present. Monthly data from NCEP and OAFlux were also obtained at the global level for comparison, with horizontal resolutions of 2.5°× 2.5°and 1°× 1°, respectively 39,52 . The subsurface temperature, current velocity and MLD were derived from the GODAS gridded product. This reanalysis dataset provides global coverage of the oceans at a monthly resolution from 1980 to the present at a 1/3°latitudinal and 1°longitudinal resolution with 40 vertical levels and a fine resolution in the upper ocean 53 . The MLD is defined as the depth where the buoyancy difference with respect to the surface level is equal to 0.03 cm s −2 53 . We also used ocean reanalysis data from the SODA global reanalysis product, which includes ocean temperature and ocean currents from 5 to 5400 m at 50 vertical levels with a horizontal resolution of 1/2°× 1/2°and a monthly temporal resolution from 1980 to the present 54 56 .

MHW definition and identification
We used monthly data instead of daily data to identify MHW events according to ref. 1 . When the SST exceeded a 90th percentile threshold based on long-term climatology and lasted for at least 2 months, it was considered an MHW event. The 90th percentile threshold was calculated for each month by using 30-year SST data (1982 to 2011), and the MATLAB "m_mhw" toolbox was used to analyze the MHWs 57 .
Ocean mixed layer heat budget Analyses of the heat budget within the mixed layer have been used to describe the processes of MHW formation, evolution, and decay 7,8,58,59 . This approach relates temperature changes in the mixed layer to physical processes, including horizontal advection, air-sea heat fluxes, and entrainment of water into the mixed layer from below. The temperature tendency in the mixed layer is calculated as follows 60 : where T is the vertically averaged mixed layer temperature, and ∂T ∂t is the change rate of the vertically averaged mixed layer temperature. Q is the air-sea heat flux into the ocean surface; ρ and C p are the reference density and specific heat of seawater (ρC p = 4.088 × 10 6 J°C −1 m −3 ), respectively; h is the mixed layer depth; u a is the vertically averaged ocean current velocity within the mixed layer; w h is the vertical velocity at the bottom of the ocean mixed layer, which is positive for upward flow; T h is the temperature at the bottom of the mixed layer; κ is the diffusivity constant, and z is the depth. The second term on the right-hand side of the equation represents the horizontal heat advection, and the third and fourth terms on the right-hand side of the equation are collectively referred to as the residuals, which are obtained by subtracting the first and second terms from the mixed layer temperature tendency. A similar approach has been applied in other parts of the global ocean, e.g., the Tasman Sea 9 , the Yellow and East China Seas 61 , and our study region 10,58 , which shows that this method can effectively identify the drivers of surface warming.
The values of the air-sea heat flux from the atmosphere to the ocean included the four following components: where Q net is the air-sea heat flux; SW and LW are the shortwave and longwave radiation flux, respectively; and LH and SH are the latent and sensible heat flux, respectively. The first two fluxes are radiative heat fluxes, and the latter two fluxes are turbulent heat fluxes.
Wind speed, air-sea temperature difference and air-sea specific humidity difference were the key factors controlling the variability in the turbulent heat flux. To investigate the key factors, we decomposed and calculated the latent and sensible heat flux anomalies using the methods published by Tanimoto et al. 31 and Fairall et al. 36 , which are calculated as follows: LH 0 ¼ ρ a LC E U 0 10 q s À q a ð Þþρ a LC E U 10 q 0 s À q 0 a À Á þ ρ a LC E U 0 10 q 0 s À q 0 a À Á where ρ a is the density of air. L is the latent heat of vaporization and a function of the SST, and it is expressed as L ¼ ð2:501 À 0:00237 SSTÞ 1:0 6 . C E and C h are the transfer coefficients of the latent and sensible heat fluxes, respectively. C p is the specific heat capacity of air, U 10 is the 10 m wind speed, q s is the saturation humidity at the SST, and q a is the specific humidity at a reference height of 2 m. T s is the SST, and T a is the air temperature at a reference height of 2 m. The three terms on the right-hand side of Eq. (3) are considered to represent the respective contributions from wind speed anomalies, air-sea temperature difference anomalies and the nonlinear effect of the two factors on the sensible heat flux anomalies. To avoid any complications during this study, we use the same terminology for the corresponding individual contributions to the anomalous latent heat flux because surface humidity anomalies generally depend on local air temperature anomalies 31 . This choice is reasonable because a good correlation (R = 0.98, above 95% confidence) was observed between the air temperature anomalies and the surface specific humidity anomalies.

Numerical model
To characterize the warming pattern and explore the process associated with warm anomalies in the NE Pacific, we used the ROMS model for the Pacific Ocean 62 . The ROMS is a free-surface, terrain-following, primitive ocean model equation 63,64 with 30 layers in the vertical direction and a horizontal resolution of 1/8°. This model was constructed for the Pacific (45°S-65°N, 99°E-70°W), and the initial and climatological temperature and salinity field data are from the World Ocean Atlas 2001. The model was forced by surface forcing fields, including the daily air-sea heat flux and freshwater obtained from the NCEP/National Center for Atmospheric Research (NCAR) reanalysis data 52 . In addition, the daily surface wind data were obtained from the NOAA multiple-satellite blended sea surface winds 65 . The model was integrated for 26 years from January 1991 to December 2016, and the outputs from 2010 to 2012 were used for climatological analysis. The sensitivity experiments were restarted and ran from January 2013 to December 2016. To evaluate the model's performance, reanalysis data from the study region were collected. As shown in Supplementary Fig. 5, the correlation coefficient between the modeled SST and the reanalysis data was 0.97, and the modeled OHC (0-150 m) matched the reanalysis data very well, with a greater correlation coefficient of 0.98; both were significant at the 95% confidence level. Overall, the modeled SST and OHC data performed well in the study region for the warming periods, allowing us to investigate the dynamics of the Blob's formation.

DATA AVAILABILITY
All data used in this study are available from publicly accessible depositories. The ERA5 atmosphere reanalysis data obtained from the ECMWF are available at cds.climate.copernicus.eu. The NCEP Atmosphere reanalysis data and the GODAS ocean reanalysis data provided by the NOAA PSL are available at: psl.noaa.gov/data/ gridded/data.ncep.reanalysis.html and psl.noaa.gov/data/gridded/data.godas.html, respectively. The OAFlux atmosphere reanalysis data are available at oaflux.whoi.edu/data-access/, the MHW identification code is available at github.com/ ZijieZhaoMMHW/m_mhw1.0, and the SODA global reanalysis dataset is available at https://www2.atmos.umd.edu/~ocean/.