Importance of ocean dynamics in the onset and persistence of the 2013-15 and 2019-20 northeast Pacific marine heatwaves

Domain: earth_science
Published: 2025-11-11
Editorial outcome: major_revision

Abstract
Large-scale marine heatwaves in the mid-latitude northeast Pacific have garnered significant attention due to their vast spatial extent, intensity, prolonged duration, and detrimental impacts on marine ecosystems and fisheries. Contrary to the conventional understanding that surface heat fluxes drive their formation, here we show that the two severe and impactful multi-season marine heatwave events in 2013-15 and 2019-20 were caused by the large-scale northward displacement of warm subtropical waters into the typically colder subarctic region. These oceanic changes are mainly explained by wind-driven circulation changes and Sverdrup balance adjustment. The marine heatwave decay phase corresponds with anomalous northwesterly winds which transport cold, dry air, enhancing latent heat loss and leading to ocean surface cooling. The physical driver is linked to the Tropical/Northern Hemisphere (TNH) teleconnection pattern, which aligns with previous studies. Collectively, the characteristic interannual timescale of the oceanic dynamic response and TNH explain the multi-season persistence of these extreme events.

Introduction
Over the past decade, two large-scale marine heatwaves in the northeast Pacific—colloquially known as the ‘Blob’—have attracted considerable scientific attention due to their unprecedented duration, spatial extent, and ecological impacts1–10. ‘Blob 1.0’ (2013–2015) emerged in the boreal winter (September–February)1, while ‘Blob 2.0’ (2019–2020) intensified in the boreal summer (May–August)6. Although interannual climate variability has been proposed as the key driver2,11, the reasons for their distinct onset month6,12 and multi-year persistence2,11,13 remain debated, highlighting gaps in our understanding of regional air-sea interactions.

Traditional mixed-layer heat budget analyses attribute the formation of these ‘Blob’ events primarily to reduced ocean heat loss to the atmosphere, with secondary contributions from oceanic advection and subsurface downward heat fluxes1,6,14–19. Subsequent numerical and observational studies have explored the influences of atmospheric forcing on surface heat flux and upper-ocean mixing6,18,19, highlighting mixed-layer depth variability16 and subsurface freshening20 as important mechanisms. However, we have identified that certain assumptions within these heat budget frameworks may obscure the important causal mechanisms.

Sea surface temperature anomalies (SSTA) in the northeast Pacific region exhibit significant interannual variability from 1980 to 2019 (Fig. 2c in ref. 2; Fig. 1b in ref. 6; Fig. 1a in ref. 12). This implies while marine heatwaves are typically defined as extreme high temperatures exceeding a threshold for at least five consecutive days21, understanding these northeast Pacific marine heatwaves requires focusing on interannual SSTA rather than raw, unfiltered data (Supplementary Note 1). On the other hand, atmospheric changes occur at higher frequencies, meaning that temperature tendencies driven by surface heat flux dominate over slower oceanic processes. Therefore, analyzing unfiltered SSTA without accounting for the timescales most relevant to extreme events can potentially bias toward the atmospheric processes.

We note that classic mixed-layer heat budget equations may inadequately assess oceanic contributions. For example, advection terms often represent heat redistribution within a domain rather than bulk water parcel movement22. Additionally, combining lateral and vertical advection into entrainment can obscure key oceanic advection effects (see “Methods”). Oceanic processes are interdependent: geostrophy influences horizontal heat transport and modulates water column thickness, affecting stratification and vertical mixing. These nonlinear interactions complicate mechanistic interpretations of linearly decomposed heat budget terms. While the standard surface mixed-layer heat budget framework (temperature tendency = surface heat flux + temperature advection + residual) is widely used to understand local temperature change, the individual components of the residual term may contain critical insights when examining larger-scaled dynamics.

Notably, the spatiotemporal scales of the 2013–15 and 2019–20 northeast Pacific marine heatwaves align with changes in the large-scale ocean dynamics. Although initial sea level anomalies were proposed as precursors23, their underlying dynamics remain unclear. The North Pacific Current (NPC), a confluence of the North Pacific subarctic and subtropical gyres, traverses regions with SSTA exceeding +2 °C24, marking the major temperature front in the northeast Pacific. We hypothesize that large-scale ocean circulation variations—particularly meridional excursions of the NPC—play a key role.

Our analysis demonstrates that large-scale interannual SSTA variations in the northeast Pacific—specifically, the dominant component of unfiltered SSTA—is primarily governed by large-scale ocean dynamics, a critical mechanism previously underestimated in the literature. Using a mixed-layer heat budget and a box model, we decompose oceanic contributions into lateral heat advection and a residual term (mainly vertical heat flux). The initiation of these two major northeast Pacific marine heatwave events, arises from anomalous wind stress curl driving northward subtropical warm water advection, whereas its decay is mechanistically linked to northwesterly wind anomalies that induce cold air advection and amplify latent heat loss. Collectively, these processes identify the Tropical/Northern Hemisphere (TNH) teleconnection pattern as the dominant driver of northeast Pacific marine heatwave dynamics25.

Results
Gyre boundary variation in the northeast Pacific
The North Pacific subtropical gyre, characterized by warm and saline anticyclonic (clockwise) circulation, converges with the cool and fresh cyclonic (anticlockwise) subarctic gyre, resulting in a northeastward flow. The eastern segment of this flow is recognized as the North Pacific Current, whose axis marks both the gyre boundary and a basin scale temperature front (Fig. 1a).

The gyre boundary (see Methods) exhibited a southward displacement during the periods of 1994–2001 and 2010–2012, while a northward shifted from 2002–2007 and after 2013, coinciding with the occurrences of the ‘Blob 1.0’ and ‘Blob 2.0’ events (Fig. 1b). At the beginning of 2013, the gyre boundary was situated near its mean position (150 °W–130 °W, Fig. 1c, blue solid line), subsequently displacing northward from 2014 to 2016. Following a brief southward displacement in 2017, the boundary reached a northern latitude of approximately 54°N during the summer of 2019, a position that was unprecedented within the satellite altimeter record from 1993 to 2022. This shift coincided with the documented intensification of ‘Blob 2.0’.

The interannual displacement of the gyre boundary from its mean latitudinal position dominates its annual cycle (Fig. 1c, red line), with a poleward trend of +1.33° latitude per decade (Fig. 1c, blue solid line). The daily SSTA (Fig. 1c, gray shading) is consistent with interannual variations of the gyre boundary. Moreover, the zero-lag correlation coefficient was found to be +0.7 (p < 0.05) between the 1-year filtered, detrended monthly SSTA and the latitude displacement anomaly of the gyre boundary.

In order to conduct a more thorough examination of the potential impact of gyre boundary displacements on the SSTA in the northeast Pacific, which are linked to variations in the basin-scale ocean circulation, we created composite maps depicting sea level and SST for the years when the gyre boundary reached its northernmost positions (1993, 2002, 2005, and 2019; henceforth referred to as the ‘positive phase’) and its southernmost positions (1998, 1999, 2011, and 2017; henceforth referred to as the ‘negative phase’) (Fig. 2).

During the positive phase, the subtropical gyre located east of 175°W expands northward, leading to a southwest-northeast oriented tilt of the gyre boundary (Fig. 2a). In contrast, during the negative phase, the gyre boundary shifts southward and presents an almost zonal orientation (Fig. 2d). The composite maps of SST exhibit a pattern that closely resembles that of sea level (Fig. 2b, e). Given that the subtropical gyre is characterized by higher sea level and SST, while the subarctic gyre is associated with lower values, the northward displacement of the subtropical gyre into the region typically occupied by the subarctic gyre results in positive anomalies in both sea level and SST. Consequently, the phase difference of sea level (Fig. 2g) shows a similar spatial pattern to that of SST (Fig. 2h).

The orientations of the gyre boundary during the ‘Blob 1.0’ event (Fig. 2c, black line) and ‘Blob 2.0’ event (Fig. 2f, black line) exhibit characteristics like those observed during the positive phase (Fig. 2a). This correspondence also extends to spatial patterns of SST (Comparison of the color shading in Fig. 2c, f with that in Fig. 2b). In addition, the correlation coefficient (Fig. 2i) between the latitude displacement of the gyre boundary (Fig. 1c, blue dashed line) and the SSTA (also detrended) aligns with the spatial pattern of the SST difference between the positive and negative phases (Fig. 2h). Notably, a high correlation value (>0.6, p < 0.05, Fig. 2i, black thick line) is observed in regions traversed by the gyre boundary.

Collectively, these findings imply that the poleward shifts of the gyre boundary substantially influence SSTA in the mid latitude northeast Pacific through lateral advection. This interpretation challenges the existing paradigm, highlighting the need to analyze the mixed layer heat budget in greater detail to understand the reason for this inconsistency.

Mixed-layer heat budget
As stated in the Introduction, it is imperative to select an appropriate timescale to undertake mixed layer heat budget diagnostic analyses (Fig. 3a). Our analysis shows the interannual component of the mixed layer temperature anomalies (MLTA) accounted for 80% of the variations observed in the unfiltered MLTA, while the ratio of their respective daily tendencies represented only 4% (see Methods). This discrepancy is attributed to the prolonged duration and consistent sign of the interannual component of the MLTA tendency. Specifically, if a 90-day low-pass filtered MLTA time series is selected for the analysis, the tendency associated with variations spanning from 90 to 365 days may obscure the tendency linked to interannual variations, potentially resulting in a misinterpretation of the prevailing physical processes (Comparison of the magnitude between Fig. 3b, c).

Collectively, understanding the physical mechanisms that govern these large-scale northeast Pacific marine heatwaves require explaining their interannual SSTA variations. Our framework reveals significant discrepancies with previous studies. For example, there was a warm anomaly in the unfiltered MLTA (Fig. 3a, red shading) from 1 June 2013 to 15 September 2013, and a transient cool anomaly from 15 September 2013 to 23 October 2013, which again turned to a warm anomaly after November 2013. The cool anomaly lasted for about 40 days and was induced by the negative surface heat flux anomaly (Fig. 3b, blue line, black double-head arrow). The previous studies did not take this short cool anomaly into account and focused on the heat budget in November 2013. Since the interannual component of the MLTA is dominant, its neutral point (the time point when the cool anomaly shifts to the warm anomaly) should be regarded as a more reliable reference point for analyzing the underlying physical mechanisms.

This approach shifts the reference point of SSTA increase for ‘Blob 1.0’ from November 2013 to March 2013, and identifies May 2018 as the reference for ‘Blob 2.0’ (Fig. 3a, black line, magenta arrow), contrasting with the previous diagnosis of May 20196. This adjustment affects the interpretation of the underlying physical drivers, particularly given the differing atmospheric conditions between March–May and November2,19.

Our interannual-scale mixed layer heat budget analysis6,26,27 partitions MLTA tendencies into atmospheric (surface heat flux) and oceanic (MLTA tendency minus net surface heat flux term) contributions (Fig. 3b, c). Our results show persistent positive oceanic contributions (Fig. 3c, black line in dashed rectangles) and negative atmospheric contributions during the formation of these marine heatwaves, indicating that oceanic processes predominantly drive the onset while surface fluxes govern decay.

To address potential limitations in using mixed layer heat budget analyses in representing the oceanic dynamics (see Methods), we employed a box model (Fig. 3d) that decomposes the oceanic term into lateral advection and a residual (primarily vertical processes). The result demonstrates the dominant role of lateral advective heat flux in the formation of these two large-scale northeast Pacific marine heatwaves, with minimal surface flux influence (Fig. 3d, black dashed rectangle), supporting our advection hypothesis. These findings raise fundamental questions: What factors cause the interannual variations in ocean circulation and net heat flux?

Wind stress curl anomalies drive gyre boundary poleward displacements
Variations in the northeast Pacific Ocean large-scale circulation are predominantly controlled by the wind stress curl28. To analyze this relationship, we generated composite wind stress curl maps by using ERA5 data during extreme northward or southward displacements of the gyre boundary, corresponding to its positive and negative phases (Fig. 4a, b). The zero wind stress curl contour closely coincides with the gyre boundary positions (magenta lines). Positive wind stress curl anomalies in the southwest and negative anomalies in the northeast of the region, delineated along a diagonal axis from 125°W-40°N to 160°W-55°N (Fig. 4c), correspond with northward displacements of the gyre. Conversely, the negative phase features a weaker and spatially reversed anomaly pattern (Fig. 4d). As expected, the annual mean wind stress curl anomalies in 2014 (Fig. 4e) and 2019 (Fig. 4f) closely mirror the pattern associated with northward gyre boundary shifts, suggesting changes in the wind stress curl were a critical driver of these large-scale northeast Pacific marine heatwave events.

The North Pacific subtropical (subarctic) gyre responds to the large-scale negative (positive) wind stress curl through Sverdrup dynamics24: βMy = curl(τ), where β represents the latitudinal gradient of Coriolis parameter, curl(τ) the vertical component of wind stress curl, and My the meridional mass transport. Negative wind stress curl anomalies enhance the subtropical anticyclonic circulation while weakening the cyclonic subarctic gyre. Although the Sverdrup relation represents the steady-state solution, our analysis does not concentrate on the transition of ocean circulation between its positive and negative phases (Fig. 2a, d), which pertains to time-varying solutions. Rather, the Sverdrup relation here explains the ocean circulation during these two phases.

Positive-phase wind stress curl anomalies induce a northeast-tilted gyre boundary displacement east of 175°W, contrasting with zonal boundary orientations during negative phases. Sverdrup-derived mass transport stream function zero contours (blue lines) align closely with the observed gyre boundaries (magenta lines) in both phases (Fig. 4).

Sea level alone inadequately represents the geostrophic stream function due to steric height contributions (sea level changes caused by thermal expansion and saline contraction). Composite analyses of the 2000 m-referenced geostrophic stream function and temperature reveal northward (39.1°N) and southward (36.4°N) surface displacements of the gyre boundary during the positive and negative phases respectively, creating a 2.7° latitudinal displacement (Fig. 5a, b). Upper 200 m stream function anomalies indicate subtropical gyre expansion and subarctic gyre contraction during the positive phases (Fig. 5d). The exceptional 2019 northward shift to 43.9°N reflects intensified anticyclonic anomalies across 30-60°N (Fig. 5f).

The spatial expansion of the subtropical gyre results in two key consequences: the warm subtropical water moves poleward, and the local thermocline deepens, causing warming within the latitude band of 30–60°N. Additionally, the deepening of the thermocline accounts for the presence of the warm anomaly beneath the mixed layer16,29 and indicates that advection plays a more important role in the subsurface warming.

Phase-dependent wind vector anomalies (Fig. 4c, d), though not directly incorporated into turbulent heat flux calculations, influence air mass sources over the northeast Pacific region. Given the role of surface heat flux in the decay of these two events, these wind pattern variations motivate investigation of the atmospheric conditions that caused the heat flux anomalies.

Anomalous northwesterly wind explains marine heatwave decay
The net heat flux term in the mixed layer heat budget (Fig. 3c, blue line) exhibits similar variations with the surface net heat flux (Fig. 6a, black line) during the negative peaks of SSTA during the decay phases of the 2013–15 and 2019–20 northeast Pacific marine heatwaves (Fig. 6, vertical magenta lines). Building on prior comprehensive analysis16 of mixed layer depth impacts on surface heat flux components, we focus on clarifying the role of the surface heat flux in the dynamics of these two marine heatwave events.

Latent heat flux dominates the net surface heat flux (Fig. 6a, blue line), with secondary contributions from sensible heat flux and longwave radiation (Fig. 6a, red and green lines). The primary driver of the latent heat flux anomalies—calculated as wind speed multiplied by the specific humidity difference between saturation and 2m air—was decomposed into wind speed variations (Fig. 6b, blue line) and humidity difference variations (Fig. 6b, red line).

Although oceanic heat advection played a critical role in the evolution of these two marine heatwave events, weakened surface heat loss occurred in 2013 and around 2018 (Fig. 6a). Latent heat flux decomposition reveals that specific humidity differences predominantly governed variations in early 2013 and late 2017 (Fig. 6b), while wind speed reductions (Fig. 6e) gained comparable importance during late 2013 and early 2018 (Fig. 6b). Throughout the 2016–2017 and 2020–2021 northeast Pacific marine heatwave decay periods, specific humidity differences remained the principal factor.

Saturation-specific humidity (Fig. 6c, red shading) and 2m air humidity (Fig. 6c, blue shading) display consistent variations, with their difference (Fig. 6c, black line) aligning closely with the humidity difference component in latent flux decomposition (Fig. 6b, red line). However, spatial patterns of wind vector anomalies (Fig. 4c, d) and directional shifts (Fig. 6d) suggest wind direction variations critically influenced both wind speed (Fig. 6e) and humidity differences.

The initiation of these two large-scale northeast Pacific marine heatwave events coincided with anomalous southeasterlies (Fig. 4c), characterized by positive anomaly of wind direction (Fig. 6e), driving warm air advection toward lower-SST regions. This induced air mass compression (reducing wind speeds, Fig. 6e), while warm, moist southern air traversing low SSTA region generated negative humidity difference anomalies (Fig. 6c), collectively suppressing ocean heat loss. On the other hand, the decay of these two marine heatwaves featured a negative anomaly of the wind direction, marked by anomalous northwesterlies (Fig. 4d), with subsequent wind speed increases attributable to cool air advection toward higher-SST areas, causing air mass expansion. These cold winds enhanced the ocean surface latent heat flux in the northeast Pacific, accelerating SSTA reduction through intensified ocean heat loss.

Schematic view and physical driver
During the onset phase of the 2013–15 and 2019–20 northeast Pacific marine heatwaves, wind stress curl anomalies drove a northward migration of the gyre boundary that served as the primary contributor to increased SSTA. In contrast, surface turbulent heat flux anomalies generated by wind direction anomalies emerged as the dominant driver of decreased SSTA during the subsequent decay phase. The distinct temporal responses of SSTA to heat advection versus diabatic cooling processes reveal a nonlinear relationship between the local upper ocean thermal response and the wind forcing (Fig. 7).

Specifically, these variations in the temperature tendency result from the interaction between oceanic advection and the surface heat flux. During the positive SSTA tendency (Trough→0→Ridge), southeasterly wind anomalies advect warm humid air into the northeast Pacific region. The condition of ‘warm winds over low-SSTA’ suppress ocean heat loss, while negative wind stress curl anomalies enable subtropical warm water intrusion, jointly elevating SSTA via advection and reduced latent heat flux. Once the SSTA exceeds the value of zero (0), the air-sea interactions transit to ‘warm winds over high-SSTA’, where advection is the dominant contributor to the amplification of the SSTA as latent flux stabilizes at the 0 anomaly (Fig. 3c). Conversely, northwesterly wind anomalies drove cold air southward in correspondence with the negative SSTA tendency (Ridge→0→Trough). Specifically, ‘Cold winds over high-SSTA’ enhanced turbulent heat loss, while positive wind stress curl anomalies caused warm water divergence out of the northeast Pacific region. Meanwhile, the eastward North Pacific Current continuously drove the inflow of warm water, causing the oceanic heat advection to remain positive. When the gyre boundary shifted past its mean position, ‘cold winds over low-SSTA’ restored the positive latent heat flux, as southward gyre displacements amplified the advective cooling.

Given the Sverdrup relation (βMy = curl(τ)) describes the linear response of the ocean circulation to the wind stress curl (Fig. 4a, b), we found that the wind stress curl over the ‘Blob’ region led gyre boundary latitude anomalies (Fig. 1c, blue dashed line) by 4 months at the maximum correlation coefficient (Fig. 8a). As theoretically anticipated, a negative correlation (peak r < −0.6, p < 0.05) manifested in the Alaskan gyre region (Fig. 8a, green box) north of the ‘Blob’ region, indicating that the northward displacement of the gyre boundary was predominantly governed by the subarctic gyre contraction. Importantly, we found positive correlations near 12°N (135°E-160°W) and negative correlations in the equatorial Pacific (0–3°N, 180–120°W), suggesting that wind fluctuations over the northeast Pacific region were most likely influenced by broader-scale atmospheric changes.

Within the Alaskan gyre region, the wind stress curl corresponds strongly and significantly with the sea level pressure anomalies (−0.81, p < 0.05, SLPAs, Fig. 8b). Lag-correlation analysis reveals that SLPA perturbations originate in the tropical eastern Pacific with 5-month lead before the Alaskan gyre intensified (Fig. 8c--h). The zero-lag correlation pattern exhibits positive SLPA values over the Gulf of Alaska and from the Gulf of Mexico northeastward across the western North Atlantic, contrasted by negative anomalies across eastern Canada. This spatial configuration bears strong resemblance to the positive phase of the Tropical/Northern Hemisphere teleconnection pattern (TNH, Fig. 8h)25, which modulates the longitudinal extension of the Pacific jet stream and influences the intensity/position of the climatological Hudson Bay Low. Our findings align with previous conclusions13 that persistent positive TNH phases facilitate prolonged ‘Blob’ marine heatwave events through sustained atmospheric forcing.

Methods
Sealevel and sea surface temperature data
The daily sea level data were obtained from the Copernicus Marine Service, with the ProductID of SEALEVEL_GLO_PHY_L4_MY_008_047. This 0.25° gridded product has global coverage. The data used here are from 1 January 1993 to 7 June 2023.

The U.S. National Oceanic and Atmospheric Administration Daily Optimum Interpolation Sea Surface Temperature (OISST V2.1) data35 are used to show the sea surface temperature in Figs. 1 and 2. These SST data are optimally interpolated onto a 0.25° × 0.25° grid. The SSTA were obtained by removing the climatological daily mean values, calculated between January 1993 and December 2021, from the SST data.

Wind and ocean reanalysis data
We also analyzed 10 m wind vectors above the sea surface based on data from the European Center for Medium Weather Forecasting (ECMWF) reanalysis, with the product name of ‘ECMWF Reanalysis v5 hourly data on single levels from 1940 to present’ (ERA5). The data used in the study are 6-hourly from 00:00, with a spatial resolution of 0.25° × 0.25°. The wind stresses are calculated following ref. 36.

The daily ocean reanalysis data used in this study are from the Bluelink ReANalysis 2020 (BRAN2020)37. These data are 0.1° × 0.1° spatial resolution extending from 75°S-75°N. The data cover the period from January 1993–September 2023.

Calculating the gyre boundary
The monthly sea level data within the region 180°W-116°W, 32°N-55°N are used to calculate the location of the gyre boundary. The longest +0.47 m contour is identified as the gyre boundary and North Pacific Current flow axis. For each month, the latitude of the gyre boundary was averaged from the original 0.25° resolution in longitude onto a 1° grid. Then, the climatological monthly mean value was removed to obtain the latitude anomaly of the gyre boundary from its mean position (color shading in Fig. 1b).

Calculating the geostrophic stream function
The geostrophic stream function as a function of latitude and depth (the sectional view)—within the longitude band 145°W–135°W and latitude band 30°N–61°N where the west coast of Alaska is located—was calculated using daily seawater temperatures and salinities from the Bluelink reanalysis 2020 (BRAN2020) and validated by comparing with the result from the IAP gridded dataset (Supplementary Note 2).

For any day, the matrix of the temperature and salinity data is in the form of longitude × latitude × depth. According to the zonal resolution of the ocean reanalysis data, the temperature and salinity data were divided into meridional sections, each section in the form of latitude × depth. Then the zonal geostrophic velocity was calculated, based on a reference level of no motion of 2000 m, or the deepest grid level where water depth was shallower than 2000 m. Then the sectional geostrophic velocity was zonally averaged. The geostrophic stream function is the horizontal integral of the geostrophic velocity from the northernmost to the southernmost extent of the section.

Heat budget calculation
The heat budget calculations. The mixed layer heat budget equation applied here is in a simple form of:
∂MLTA / ∂t = Qnet / (ρCp h) + residual (1)

where Qnet is the net heat flux (the sum of shortwave radiation, longwave radiation, latent and sensible heat flux), Qh is the penetrative shortwave radiation (based on the relationship of:
Qh = Qsw * (0.58e^(-h/0.35) + 0.42e^(-h/23)) (2)

where Qsw is the shortwave radiation, h is the time varying mixed layer depth, and ρ and Cp are the seawater density and specific heat capacity, respectively. Here, ρCp = 4.088 × 10^6 J °C^-1 m^-3, following ref. 16, and the residual term comprising the ocean dynamic contributions calculated as the difference between ∂MLTA/∂t and Qnet/(ρCph).

The temporally and spatially varying mixed layer depth, from the model output, is defined as the first depth where its buoyancy is 0.0003 m/s^2 different from the surface. The analysis is calculated over the region of 150°W-135°W, 40°N–50°N, and is then averaged by considering the area of each grid cell (areal mean value).

The climatological daily mean value is removed by assuming the equation is balanced for the annual cycle.

The box model quantifies the total temperature tendency (∂BOXTA/∂t) through six interfaces of a predefined cubic control volume in the upper ocean. The model assumes the heat flux into the box is positive and is expressed as:
∂BOXTA / ∂t = Qnet / (ρCp H) - (1 / (ρCp H)) * ∬(uL TL dz) dl + ∬(wTB dxdy) + mixing (3)

where H is the thickness of the box (in our case it is the total thickness of the upper 10 grids, with a value of about 61 m), uL is the horizontal velocity normal to the boundaries of the computational domain, by assuming inflow is positive, T is the seawater temperature at the boundaries. The accumulated horizontal advection through the four lateral interfaces (west, east, north, south) is calculated as:
∫(uL TL dz) dl = ∬(uw Tw dydz) + ∬(ue Te dydz) + ∬(us Ts dxdy) + ∬(un Tn dxdy) (4)

with uw/ve and ue/vn representing zonal (east−west) and meridional (north−south) velocities at the corresponding interfaces, Tw, Te, Ts, Tn as temperatures. The vertical advection term is derived from wTB (vertical velocity × temperature at the bottom interface). Each term in the box model is represented as a net heat flux (velocity × temperature, e.g., (1/H) * ∫(uL TL dz) dl and (1/H) * ∬(TB dxdy)) at the boundaries, rather than an effective heat flux (e.g., velocity × temperature gradient, e.g., -u * ∇T and (1/H) * T * (∂h/∂t + u * ∇h + w)).

The high correlation coefficients for the total temperature tendencies (+0.93, p < 0.05, between the red line in Fig. 3d and red shading in Fig. 3c), net heat flux terms (+0.84, p < 0.05, between the blue lines in Fig. 3d and Fig. 3c), and the respective residual terms that represent the oceanic contribution (+0.75, p < 0.05) suggest that the decomposition of the box model reflects the same physical processes as that of the classic mixed layer heat budget.

Examining the entrainment term. We analyzed the complete form of the classic mixed-layer heat budget equation to show why this equation cannot fully capture the advection process.

The complete equation of the classic mixed layer heat budget is:
∂MLTA / ∂t = Qnet / (ρCph) - u * ∇T + Mixing - (1/h) * (T - Th) * (∂h/∂t + u * ∇h + w) (5)

The second term on the right-hand side represents advection, while the fourth term represents entrainment. If we take the entrainment term apart, we can get,
(1/h) * (T - Th) * (∂h/∂t + u * ∇h + w) = (1/h) * T * w + (1/h) * Th * w - (1/h) * T * (∂h/∂t + u * ∇h) (6)

where T is the surface temperature and Th is the temperature at the mixed layer depth. The term (1/h) * T * w represents the vertical advection, which is induced by the divergence or convergence of the horizontal velocity.

Moreover, the entrainment term in the classic mixed layer heat budget is in the form of the effective flux, whereby the temperature difference (T - Th) multiplies vertical velocity (the term of ∂h/∂t + u * ∇h + w has a unit of m/s). Given that the temporal variation of the mixed layer depth (∂h/∂t) may also be affected by the vertical motion, the entrainment term likely contains greater contribution from the lateral advection. Therefore, the advection term in this equation may underestimate the contribution of the advection process.

Another simple way to understand this issue is by applying the classic mixed layer heat budget equation to a cubic control volume. Consider a cubic control volume with fixed dimensions (width, length, and thickness), similar to a basic grid cell in a numerical model. Since the volume is time-invariant, the ∂h/∂t and ∇h terms vanish. This leads to:
∂BOXTA / ∂t - Qnet / (ρCpH) = -u * ∇T - (1/H) * (T - Th) * w + Mixing (7)

Expanding the second term in the right-hand side apart and rearranging, we have:
∂BOXTA / ∂t - Qnet / (ρCpH) = -u * ∇T - (1/H) * T * w + (1/H) * Th * w + Mixing (8)

where -u * ∇T - (1/H) * T * w is the lateral advection and (1/H) * Th * w represents the vertical heat advection.

A simple scale analysis could help to evaluate the relative importance of each term. In the ‘Blob’ region, the average horizontal velocity at the sea surface is 0.05 m/s, the maximum horizontal temperature gradient is 9 * 10^-6 °C/m (1 °C/110 km * 0.001 km/m), the vertical velocity in the upper ocean is 10^-4 m/s, the annual mean mixed layer depth is 50 m, SST = 12 °C, Th = 12 – 0.3 °C = 11.7 °C (the mixed layer depth is commonly defined as the depth at which temperature first becomes 0.3 °C different from the surface), and suppose T = SST + Th = 11.85 °C.

On the right-hand side of the above equation, the first term (-u * ∇T) is on the order of 4.5 * 10^-7 °C/s, the second term (-(1/H) * T * w) is on the order of 2.37 * 10^-5 °C/s, and the third term ((1/H) * Th * w) is on the order of 2.34 * 10^-5 °C/s.

Therefore, both the lateral advection (-(1/H) * T * w) and the vertical advection at the bottom of the mixed layer ((1/H) * Th * w) play important roles in the heat budget equation. However, if we write them in the form of the effective flux (-(1/H) * (T - Th) * w, temperature difference multiplying vertical velocity), the two terms cancel out and generate a small term that is on the order of -3 * 10^-7 °C/s (-2.37 * 10^-5 °C/s + 2.34 * 10^-5 °C/s), which is comparable to the term of -u * ∇T. This analysis shows that if the predominance of a particular term is assessed based on its magnitude, both the vertical and lateral processes may be undervalued when the entrainment term is represented as the effective flux.

Decomposition of the latent heat flux. The latent heat flux can be represented as:
QLH = ρL C * (qs - qa) * W (9)

where W is the wind speed, qs is the saturation specific humidity, and qa is the 2 m air humidity. The latent heat flux variations can be decomposed into the contributions of the major factors:
(qs - qa) * W = (qs - qa) * W + W * (qs - qa) + residual (10)

Where tilde denotes the anomaly while overline denotes the annual mean. The climatological mean was subtracted from each term obtained from the above equation. Subsequently, a 365-day low-pass filter was implemented.

Interannual component contribution of the MLTA to the total MLTA. The interannual component contribution of the MLTA to the total MLTA is defined as the standard deviation of the interannual component of the MLTA (black line in Fig. 3a) divided by the standard deviation of the total MLTA (red shading in Fig. 3a). The interannual component contribution of the MLTA tendency to the total MLTA tendency is the standard deviation of the interannual component MLTA tendency (red shading in Fig. 3c) divided by the standard deviation of the total MLTA tendency (tendency of red shading in Fig. 3a).

Data availability
We have used publicly available data only; new data were not generated by this study. ADT: https://data.marine.copernicus.eu/product/SEALEVEL_GLO_PHY_CLIMATE_L4_MY_008_057/description. OISST: https://www.ncei.noaa.gov/products/optimum-interpolation-sst. ERA5: https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-single-levels?tab=overview. BRAN2020 Ocean reanalysis: https://research.csiro.au/bluelink/bran2020-data-released/

Code availability
The code used to analyze these data and generate the results presented in the study can be obtained from Supplementary Dataset 1–3.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This manuscript reports on an investigation of a novel hypothesis that a dynamic ocean response to wind stress curl forcing played a fundamental role in two recent multi-year NE Pacific marine heatwaves. The authors conclude that interannual variations in gyre scale frontal dynamics comprehensively explain the spatial characteristics and multi-year persistence of the two events. The article is interesting for offering up this idea and some provocative interpretations of the results of their analyses. However, I do not believe that the work detailed here supports such grand conclusions. As noted in my specific comments below, simply showing that these events featured meridional excursions of the subarctic frontal boundary falls short of demonstrating cause and effect. Moreover, there is really no analysis or theory offered for linking the proposed mechanisms to the observed time-scales involved in the multi-year marine heatwaves of interest. If the event duration times were determined by the atmospheric forcing, for instance, one would need to explain the observed evolution of the atmospheric forcing field (something that DiLorenzo and Mantua focused on for the 2014-15 period, and Amaya et al. focused on for summer 2019). This work has focused more on a time-varying but spatially invariant interpretation of the atmospheric forcing - for example, creating composite maps (Figs 2 and 3) for sea level, SST, and the surface wind stress fields during what are termed “positive” and “negative” phase years in their detrended and 1-yr filtered gyre boundary index (which has its maximum value in 1993, a year that is not used in the composites).\n\nThis work has raised some interesting questions for me that I think are worth addressing. These include:\n• The monthly mean gyre boundary index shows a substantial trend over the 1993-2022 period - is there a corresponding trend in the surface wind stress field (for instance, in the zero-curl contour latitude in the NE Pacific)?\n• How would the gyre strengths and positions vary in response to the kind of atmospheric forcing variations described by Di Lorenzo and Mantua (2016)? Would they include the expansion of the subtropical gyre in 2013 through winter 2014, or a slowing of both the subarctic and subtropical gyres in response to the atmospheric forcing in the “North Pacific Oscillation” pattern that DiLorenzo and collaborators have linked with the North Pacific Gyre Oscillation?\n• Would the fall 2014/winter of 2015 feature an expansion of the subarctic gyre but contraction of the subtropical gyre? These SSH variations and their links to the gyre circulations and atmospheric forcing are basically those described by Cummins and Lagerlof (2004), your reference 20.\n• Could you better evaluate the relative contributions of basin-scale ocean dynamics versus surface flux forcing through simulations that involved forced ocean responses in a slab-ocean (mixed layer only) versus a dynamically-responding ocean model?\nIn summary, I believe there is merit in the research directions this work reports on. However, I do not believe that the research reported here supports the conclusions that are offered. I recommend this manuscript be returned to the authors with a request for major revisions.\nNate Mantua\nNOAA/NMFS Southwest Fisheries Science Center\nSanta Cruz, CA\nSpecific comments:\n40: I would add “and fisheries” after “marine ecosystems”\n43-44: revise to “meridional expansion (contraction” of the subtropical (subarctic) gyre” … current text doesn’t include subarctic\n45-47: this analysis does not isolate the contribution to gyre scale frontal dynamics from other processes, nor does it connect those dynamics to the time scale of the events, so it falls short of “comprehensively explain(ing) the spatial characteristics and multi-year persistence of these … heatwaves”\n60: reasons for multi-year persistence have been investigated and explained by others (e.g., DiLorenzo and Mantua 2016)\n78-80: Isn’t a dominant tendency term required to understand the origins of a broad-scale warming of the upper ocean that happens within a single season, hence the focus on atmosphere-ocean interactions?\n83-86: This may be true in an ideal case, but this real-world case study does not involve a surface mixed layer having top and bottom interfaces that are adiabatic.\n120-125: Causation may go either way, or there may be a common forcing for the SSTa and latitude of the gyre boundary. I also think that the trend over the 1994-2023 period of record is striking. Is there is a similar trend in atmospheric forcing? Why doesn’t the daily SSTA time series share the longer-term trend with the gyre boundary index?\n131-134: Why do you conclude that the subarctic gyre west of 150W expands southward? I can’t see the southward expansion of the subarctic gyre west of 150 from the SLH or SST difference maps (panels e and f) – do you mean the negative SSH differences between 40-50N and 180-170W?. I don’t see this in Fig 1b either (I would expect to see a dipole centered on ~150W, but that is only a feature of a few years like 2014-2016). The main changes in frontal position are east of 140W from 2014-2016, and east of 165W from 2018-2020 (based on Fig 1b).\nFigures 2 and 3: I think it would be best to use the same projection in the maps shown in Figures 2 and 3.\n168-170: Interesting that the zero curl contour at 180 is about 5 deg latitude farther south in your positive phase than negative phase composite, while it is at about the same latitude for both phases on the west coast of N. America ~125W; This seems at odds with the SSH composites in Fig 2\n176-179: Fig 3a shows your “positive phase” having the gyre boundary east of 150W shifting northward, not southward, in response to a negative wind stress curl anomaly in the Gulf of Alaska (north of 50N).\n180-181: This is confusing. Perhaps you could include the actual wind stress curl anomalies for the composites of positive and negative phases and this would be easier for readers to see and understand than trying to assess anomalies from the composites shown here.\nFigure 3 a-b-c and Fig 5 a-b: Wind vectors are really small and difficult to see, even at 200% magnification. Please make them bigger.\n211-212: Thermocline deepening is inconsistent with previously published findings (like Amaya et al 2020 that show a shoaling mixed layer depth in summer 2019)?\n241-249: I don’t find this new definition of event onset to be very useful or convincing. Is there something special about climate dynamics in the month of May that makes you think this event onset criteria leads to a new and better understanding for NEP marine heatwave events? Would you say that Blob 3.0 started in January 2022 by your criteria? Did the 2013-14 event end in 2015 or 2016 by your criteria?\n276-277: Did the marine heatwaves “result from” or “include” poleward expansion (contraction) of the subtropical (subarctic) gyre?\n281-282: This analysis does not address the issue of “persistence”, but instead focuses on “onset”. What sets the time scale for these events? Is it the persistence and pattern of atmospheric forcing (along the lines of Di Lorenzo and Mantua’s analysis of the 2013-2015 period)? Or something else?\n295-297: If the basin-scale wind-stress curl is the key driver for the gyre variations, then multi-year persistence of these marine heatwaves is due to persistence in the wind-stress curl forcing. So this then begs the question – what provides persistence in the basin-scale atmospheric forcing of the gyre circulation? As noted, the ocean circulation (and mixed-layer heat content) response integrates atmospheric forcing so it has greater variance at longer time scales than the atmospheric forcing, yet the evolution of the atmospheric forcing has already been shown to play a key role in understanding the evolution and persistence of the 2014-2016 (DiLorenzo and Mantua) and 2019-2020 (Amaya et al) events.\n324-328: Why would a ML heat budget fail to capture warm water advection into the “Blob” region? Are the boxes drawn with boundaries that lie outside that of the warm water advection process? Wouldn’t frontal shift-induced SSTA be maximized along the axis of the SST front, and zonally confined to the SST front?\n342-343: Is the 30 deg latitude extent of the SSTA warm pattern consistent with the 200km (2 deg latitude) gyre boundary excursion? What about the scale of atmospheric forcing anomalies associated with surface heat flux anomalies?", "reviewer_2": "This study tried to explore the potential role of gyre scale ocean dynamics in iconic northeast Pacific marine heatwaves, focusing on the interannual or low-frequency timescales. I think this viewpoint is very important and potentially crucial for understanding the climate variability in the Northeast Pacific. But I think the current presentation is quite preliminary in many aspects. Please see my detailed comments below.\nMajor comments:\nIn your mixed-layer heat budget equation, the H is time varying mixed layer depth. Please specify how you calculate the surface heat flux term, that is Qnet/(ρCpH). In particular, how does the mixed layer depth vary in this term? Moreover, I cannot tell exactly how do the gyre variations work on the MLTA in the blob region, because authors merged the oceanic processes into one term. The physical picture is not clear to me.\nIn many previous studies, as you also shown, the oceanic processes may not be that important for the warm blobs. How to interpret your results?\nIn this study, authors did not pay special attention to the differences of mechanisms between seasons.\nWhat is the frontal dynamics? I think authors did not use too much frontal theory in this study.\nFigure 2: I think there should be a more objective method to select northernmost and southernmost years. Then, significance test should be added to confirm the robustness of your results. Currently, the sample numbers are small. Moreover, authors can add the climatological +0.47 line to show the changes. Same issue for Fig. 3. Currently, I cannot understand well the texts in Line 135-137.\nLine 215-218: Previous studies claimed the shallower mixed layer depth during the Northeast Pacific warm blobs/marine heatwaves (Amaya et al. 2020 NC; Shi et al. 2022 GRL). Could authors give more explanations on the deepening of thermocline here.\nThe neutral point appears a bit arbitrary. Authors tried to raise a new method to identify the onset of the blob, but they should compare in details with the previous results. Why this method is more reasonable that others? Moreover, did Amaya et al. (2020) propose the onset of blob 2.0 in May 2018? I think they only focused on the summer of 2019, rather than 2018.\nMinor comments:\nThe word “Iconic” is a bit hard to clearly understand. Could authors change into another word?\nAbstract, Line 43: anomalous positive or negative wind stress curl drives the expansion of the subtropical gyre. Please specify here.\nBlob 1.0 and 2.0: It is not appropriate to appear in the abstract without any background here. Readers may not know their similarities and differences.\nAbstract, Line 52: It is better to specify the sea surface temperature, rather than temperature.\nLine 61: Chen et al. (2021) classified two types of warm blobs and reported different mechanisms of the warm blobs in different seasons.\nChen, Z., Shi, J., & Li, C. (2021). Two types of warm blobs in the Northeast Pacific and their potential effect on the El Niño. International Journal of Climatology, 41(4), 2810–2827.\nLine 64: For Blob 1.0, due to its multi-year persistence over 2013-2015, the growth months may be different in these years.\nThere are some errors in the reference list, such as the missing page information, etc.\nLine 85-86: Why the advection is zeros here? Could you please draw a simple diagram to illustrate? Moreover, why do you give this ideal case? I feel confused to understand. How does advection outside the computational domain explain the temperature tendency in the domain?\nLine 98-99: Please cite some related papers to confirm this.\nLine 391-392: I am wondering whether +0.47 is reasonable for the identification of flow axis in each year/month.\nFig. 1a: climatological mean? Annual-mean? Please specify.\nFig. 1b: Could authors add significance test for this figure? The while shading to the rightmost area indicates there is no +0.47 value, right?\nFig. 1c: the interannual component should be frequency higher-than 9 years. Why do authors use the 1-year filter here?\nFig. 3: Authors did not make it clear the wind vectors and curl are anomalies or original values.\nLine 178: confusing for “In the mean”\nLine 182 and 184: west of or east of? Moreover, I cannot distinguish the changes without climatological position superimposed. Also for the right column in terms of stream function. Moreover, why these panels are averaged over a narrow band as indicated in Line 193?\nLine 203-204: Please give the figures, maybe in the supplement.\nLine 276-277: Confusing sentence for “the increased mixed layer heat loss to the air is the response”.\nLine 337: what did “residual” refer to here?\nLine 350-351: What variable for the spatial extent of approximately 30 ° in latitude?\nFig. 5: What is the difference between left and right figures in each panel? In panel c, why does the warm anomaly correspond to cyclonic wind anomalies? PV dynamics is not included in the main texts.", "reviewer_3": "This study analyzed the causes of the \"blobs\" that occurred in the Northeastern Pacific through a new perspective on large-scale ocean dynamics. The study argued that the positive-negative wind stress curl anomalies in the Northeastern Pacific drove the formation of the blob by shifting the gyre boundary between the subtropical and subarctic gyres. The study has investigated the cause of blobs through a simple heat budget analysis with respect to the gyre boundary variations. Through the negative wind stress curl anomaly, the subarctic gyre contracted, and the subtropical gyre expanded, leading to oceanic heat gain.\nThe value of this study lies in its examination of blob formation through the oscillations of subtropical and subarctic gyres, rather than focusing solely on the climate variability perspective of air-sea interactions, which was proposed in previous studies. This new viewpoint could provide insights for future research on marine heatwaves and ecosystems in other regions, making it a relevant contribution to the scope of the journal.\nHowever, there are some limitations in the analysis and interpretation of results, and I believe addressing these will significantly enhance the quality of the paper. Here are some major comments:\n1. The heat budget analysis is overly simplified. The study divides surface net heat flux and the residual term, attributing the residual term to the ocean's influence. However, the equation presented in the paper assumes that all shortwave radiation enters and exits the mixed layer, while some of it is likely transferred to the subsurface layer beneath the mixed layer. This could result in an overestimation of surface net heat flux and an underestimation of the ocean’s influence, which might align with the paper's intent, but accurate calculations and proposals of values would be more appropriate. Additionally, the entrainment term could be helpful to understand the effect of the upwelling & mixed layer depth changes on the blob formation due to the negative wind stress curl anomaly. Relevant references include:\nReferences:\nQiu, C., Kawamura, H., Mao, H., & Wu, J. (2014). Mechanisms of the disappearance of sea surface temperature fronts in the subtropical North Pacific Ocean. Journal of Geophysical Research: Oceans, 119(7), 4389-4398.\nMurata, K., Kido, S., & Tozuka, T. (2020). Role of reemergence in the central North Pacific revealed by a mixed layer heat budget analysis. Geophysical Research Letters, 47(13), e2020GL088194.\n2. The positive-negative wind stress curl anomalies weaken (strengthen) the subarctic (subtropical) gyre. According to the Ekman theory, this implies a reduction (intensification) in divergence (convergence) within the surface mixed layer, leading to heat accumulation at the surface and a weakening (strengthening) of upwelling (downwelling), reducing heat loss. Therefore, in the gyre’s interior region, the blob formation might be more due to the weakening (strengthening) of the subarctic (subtropical) gyres affecting mixed layer depth and upwelling, (downwelling) rather than advection of the gyre boundary especially. The analysis in the paper does not sufficiently support the idea that heat is transferred from the subtropical gyre to the subarctic gyre. Additional analysis showing the spatiotemporal transmission of heat from the subtropical gyre across the front to the subarctic gyre is needed. This is the reason for that the entrainment term is included in the heat budget analysis.\nAdditionally, the schematic diagram in Figure 5 suggests that warm water columns move northward via potential vorticity (PV) conservation from the subtropical gyre to the subarctic gyre across the front. However, the negative wind stress curl occurs in the subarctic gyre, not the subtropical gyre, so the figure does not align with the content of the paper. In PV conservation, an increase in H does not necessarily lead to northward movement to higher latitudes but could instead increase relative vorticity. And it is also questionable whether the water column would be able to cross the front and move northward to conserve PV. Therefore, I believe Figure 5 needs significant revision, and if the content of Figure 5 is correct, sufficient discussion should be included in the manuscript.\n3. The positive phase period should include not only the year the blob formed but also the years when the blob existed. Therefore, the years 2015 and 2020 should also be included. Additionally, does the negative wind stress curl anomaly occur in both 2015 & 2020 and does the marine heatwaves appear in 2005?\n4. In this paper, it is argued that Blob 1.0 started in May 2013, based on 365-day lowpass filtered data. Marine heatwaves are phenomena typically observed on a daily scale, so it is questionable whether identifying them using low-pass filtered data is appropriate. Additionally, an explanation is needed as to why a cold anomaly lasting 40 days was included in the Blob period, and what advantages this inclusion brings.\nMinor comments\n- Line 138. ‘the gyre boundary shifts southward (eastward),’ Please verify the sentence.\n- Line 182-184, ‘resulting in the gyre boundary east of 150W ~ , resulting in the gyre boundary west of 150W shifting northward.’ Please verify the sentence."}, "author_response": "Thank you very much for your thoughtful review. We have now revised our manuscript based on your comments and those of the other two reviewers. Our detailed responses are provided below, with your original comments in black font and our replies in blue font.\nIn the revised manuscript, we have endeavored to contextualize our analysis and proposed mechanisms within the framework of existing theories. Overall, the heat budget analysis indicates that the advection term is responsible for the onset of the Blobs, while the increased surface net heat flux from the ocean to the atmosphere accounts for its decay.\nWe employed the Sverdrup relation to elucidate the oceanic dynamic response to the anomalous wind stress curl. The Sverdrup relation serves as a foundational theory in ocean general circulation, delineating the relationship between meridional mass transport and the thickness of the water column. Although this relationship is effectively time-independent, the composite maps are effective for interpreting the quasi-steady-state (or equilibrium) responses to differences in the wind stress curl, and has been validated in some previous studies (Jiang et al., 2006; Hill et al., 2008; Thomas 2014).\nSpecifically, during the positive phase, the negative anomaly in wind stress curl over the subarctic gyre results in a weakening of the gyre, leading to its contraction and facilitating the northward intrusion of warm subtropical waters. This explanation contrasts with the previous understanding that reduced ocean heat loss to the atmosphere played a major role in the formation of the Blobs.\nWe find the anomalous latent heat flux is primarily responsible for the decay of the Blobs. It results from the increase in latent heat flux loss from the ocean to atmosphere, linked to the ocean thermodynamic response to the northwesterly surface wind anomaly that transports cold air.\nDue to significant changes in the manuscript compared to the previous submission, we have summarized the revisions here. The Introduction has been completely rewritten to provide a more comprehensive explanation of our concerns about the classic mixed layer heat budget analysis applied previously. The 'Mixed layer heat budget' section has been restructured to follow the 'Gyre boundary variation in the 'Blob' region', highlighting the importance of advection in the heat budget. Furthermore, we have added two new sections: one examining the variations in net heat flux and another addressing the physical drivers of the Blobs. The final section has also been revised to incorporate our latest findings."}, {"round": 2, "reviewer_comments": {"reviewer_1": "Reviewer #1 (Remarks to the Author):\nThis manuscript reports on an investigation of a novel hypothesis that a dynamic ocean response to\nwind stress curl forcing played a fundamental role in two recent multi-year NE Pacific marine\nheatwaves. The authors conclude that interannual variations in gyre scale frontal dynamics\ncomprehensively explain the spatial characteristics and multi-year persistence of the two events. The\narticle is interesting for offering up this idea and some provocative interpretations of the results of their\nanalyses. However, I do not believe that the work detailed here supports such grand conclusions. As\nnoted in my specific comments below, simply showing that these events featured meridional\nexcursions of the subarctic frontal boundary falls short of demonstrating cause and effect. Moreover,\nthere is really no analysis or theory offered for linking the proposed mechanisms to the observed time-\nscales involved in the multi-year marine heatwaves of interest. If the event duration times were\ndetermined by the atmospheric forcing, for instance, one would need to explain the observed evolution\nof the atmospheric forcing field (something that DiLorenzo and Mantua focused on for the 2014-\n15 period, and Amaya et al. focused on for summer 2019). This work has focused more on a time-\nvarying but spatially invariant interpretation of the atmospheric forcing - for example, creating\ncomposite maps (Figs 2 and 3) for sea level, SST, and the surface wind stress fields during what are\ntermed “positive” and “negative” phase years in their detrended and 1-yr filtered gyre boundary index\n(which has its maximum value in 1993, a year that is not used in the composites).\nThank you very much for your thoughtful review. We have now revised our manuscript based on\nyour comments and those of the other two reviewers. Our detailed responses are provided below, with\nyour original comments in black font and our replies in blue font.\nIn the revised manuscript, we have endeavored to contextualize our analysis and proposed\nmechanisms within the framework of existing theories. Overall, the heat budget analysis indicates that\nthe advection term is responsible for the onset of the Blobs, while the increased surface net heat flux\nfrom the ocean to the atmosphere accounts for its decay.\nWe employed the Sverdrup relation to elucidate the oceanic dynamic response to the anomalous\nwind stress curl. The Sverdrup relation serves as a foundational theory in ocean general circulation,\ndelineating the relationship between meridional mass transport and the thickness of the water column.\nAlthough this relationship is effectively time-independent, the composite maps are effective for\ninterpreting the quasi-steady-state (or equilibrium) responses to differences in the wind stress curl, and\nhas been validated in some previous studies (Jiang et al., 2006; Hill et al., 2008; Thomas 2014).\nSpecifically, during the positive phase, the negative anomaly in wind stress curl over the subarctic gyre\nresults in a weakening of the gyre, leading to its contraction and facilitating the northward intrusion of\nwarm subtropical waters. This explanation contrasts with the previous understanding that reduced\nocean heat loss to the atmosphere played a major role in the formation of the Blobs.\nWe find the anomalous latent heat flux is primarily responsible for the decay of the Blobs. It results\nfrom the increase in latent heat flux loss from the ocean to atmosphere, linked to the ocean\nthermodynamic response to the northwesterly surface wind anomaly that transports cold air.\nDue to significant changes in the manuscript compared to the previous submission, we have\nsummarized the revisions here. The Introduction has been completely rewritten to provide a more\ncomprehensive explanation of our concerns about the classic mixed layer heat budget analysis applied\npreviously. The 'Mixed layer heat budget' section has been restructured to follow the 'Gyre boundary\nvariation in the 'Blob' region', highlighting the importance of advection in the heat budget. Furthermore,\nwe have added two new sections: one examining the variations in net heat flux and another addressing\nthe physical drivers of the Blobs. The final section has also been revised to incorporate our latest\nfindings.\nReferences\nJiang, H., Wang, H., Zhu, J. & Tan, B. Relationship between real meridional volume transport and\nSverdrup transport in the North Subtropical Pacific. Chinese Science Bulletin 51, 1757–1760\n(2006).\nHill, K. L., Rintoul, S. R., Coleman, R. & Ridgway, K. R. Wind forced low frequency variability of\nthe East Australia Current. Geophysical Research Letters 35, (2008).\nThomas, M. D., De Boer, A. M., Johnson, H. L. & Stevens, D. P. Spatial and temporal scales of\nSverdrup balance*. Journal of Physical Oceanography 44, 2644–2660 (2014).\nThis work has raised some interesting questions for me that I think are worth addressing. These include:\n• The monthly mean gyre boundary index shows a substantial trend over the 1993-2022 period - is\nthere a corresponding trend in the surface wind stress field (for instance, in the zero-curl contour\nlatitude in the NE Pacific)?\nThank you for your comment. Yes, there is a corresponding trend in the atmospheric field.\nIn the Alaskan gyre region where the correlation coefficient between the gyre boundary and wind\nstress curl reaches a maximum (140-152 °W, 50-57 °N), we found a noticeable linear decrease in the\nabsolute values of the surface wind stress curl with -3.9×10-9 N/m3 per decade, and in the sea level\npressure anomalies with a value of 79 Pa per decade (Figure R1-1; we also added this figure as Fig.\n7b in our revised manuscript). This suggests in the region where the wind stress curl anomaly has the\ngreatest influence on the gyre boundary, there exists a corresponding trend in the surface wind stress\nfield.\nHowever, the zero-curl contours in the NE Pacific are winding and with many smaller scale\nfeatures, even using the annual mean data (Figure R1-2, black contours). Consequently, we selected\nthe longest zero-curl contour (Figure R1-2, black thick contours) and averaged its latitude to a\nresolution of 1 degree in longitude. Subsequently, we computed the linear trend of the zero-curl contour\nand identified a northward displacement, characterized by a trend of 1 degree per decade (Figure R1-\n3).\nFigure R1-1: Areal mean wind stress curl (red line) and sea level pressure (blue line) over the region\n(green boxes) where the correlation coefficient between the gyre boundary latitude and wind stress\ncurl anomaly reaches maximum.\nFigure R1-2: The yearly mean wind stress curl. Black contours denote the zero-curl contour and black\nthick contours denote the wind stress curl selected to calculate the zero-curl contour latitude.\nFigure R1-3: The linear trend of the zero-curl contour latitude. The red dots represent the point where\nthe p-value is less than 0.05.\n• How would the gyre strengths and positions vary in response to the kind of atmospheric forcing\nvariations described by Di Lorenzo and Mantua (2016)? Would they include the expansion of the\nsubtropical gyre in 2013 through winter 2014, or a slowing of both the subarctic and subtropical gyres\nin response to the atmospheric forcing in the “North Pacific Oscillation” pattern that DiLorenzo and\ncollaborators have linked with the North Pacific Gyre Oscillation?\nThank you for your comment. Di Lorenzo and Mantua (2016) did not provide a distinct pattern of\nwind stress curl; rather, they showed the seasonal evolution of sea level pressure anomalies in their\nFigure 1. Our analysis of the time series of gyre boundary latitude indicates that interannual variations\nare significantly greater than seasonal variations (Figure 1c, comparison between the red and blue\nlines). This implies the ocean dynamic response to the atmospheric driving field favors the interannual\nvariation, rather than intra-annual variation.\nIt has been previously argued that the ocean's response to atmospheric dynamics in the Northeast\nPacific has a damping time scale of 1-3 years (Cummins & Lagerloef, 2004; Cummins & Freeland,\n2007; Cummins & Masson, 2018). In this context, it would seem less appropriate to discuss changes\nin ocean circulation based on the atmospheric forcing pattern averaged over a specific month or a few\nmonths unless it designates the state of slow variations on the timescale over one-year. If your concern\nis that the pattern may present the interannual variation, please refer to our response to the following\ncomment.\nDi Lorenzo and Mantua (2016) reported, “We have shown that both the 2014 and 2015 SSTa\npatterns can be reconstructed independently using direct atmospheric forcing functions. This suggests\nthat ocean internal processes may not have played a key role in persisting and evolving the record\namplitude of the SSTa through 2015.” According to their analysis of the AR1 model, it was concluded\npreviously that ocean dynamics might not have been a significant factor. The previous research does\nnot account for the expansion of the subtropical gyre from 2013 to winter 2014, nor the slowing of\nboth the subarctic and subtropical gyres as far as we can establish here.\nReferences\nCummins, P. F. & Lagerloef, G. S. E. Wind-driven interannual variability over the northeast Pacific\nOcean. Deep Sea Research Part I Oceanographic Research Papers 51, 2105–2121 (2004).\nCummins, P. F. & Freeland, H. J. Variability of the North Pacific Current and its bifurcation. Progress\nin Oceanography 75, 253–265 (2007).\nCummins, P. F. & Masson, D. Low-frequency isopycnal variability in the Alaska Gyre from Argo.\nProgress in Oceanography 168, 310–324 (2018).\n• Would the fall 2014/winter of 2015 feature an expansion of the subarctic gyre but contraction of the\nsubtropical gyre? These SSH variations and their links to the gyre circulations and atmospheric forcing\nare basically those described by Cummins and Lagerlof (2004), your reference 20.\nThank you for your comment.\nIf we consider the atmospheric forcing pattern from fall 2014 to winter 2015 as a representation\nof interannual variation, it can influence ocean circulation changes. Previous studies have explored\nhow wind patterns affect circulation in this area (Cummins & Masson, 2018; Hristova et al., 2019),\nso\nwe will only provide a brief overview. According to ocean general circulation theory, specifically here\nthe Sverdrup relation, the North Pacific subtropical gyre is forced by negative wind stress curl, while\nthe subarctic gyre is influenced by positive wind stress curl. For the interannual variation, a reduction\nin the negative wind stress curl over the subtropical gyre is expected to result in either a weakening or\ncontraction of the gyre. Conversely, an increase in the positive wind stress curl over the subarctic gyre\nis anticipated to lead to either a strengthening or expansion of that gyre. This relationship also holds\ntrue in the reverse scenario. Furthermore, negative sea level pressure (SLP) anomalies typically align\nwith positive wind stress curl, and vice versa.\nThe uniformly negative SLP anomaly over the northeast Pacific (atmospheric forcing pattern from\nfall 2014 to winter 2015) leads to a positive wind stress curl anomaly, which increases the positive\nwind stress curl in the subarctic gyre and decreases the negative wind stress curl in the subtropical\ngyre. This results in the expansion of the subarctic gyre and the contraction of the subtropical gyre,\neffectively shifting the gyre boundary southward. Our findings indicate that there was indeed a\nsouthward shift in the gyre boundary during fall 2014/winter 2015 (Fig. 1c, blue solid line).\nReferences\nCummins, P. F. & Masson, D. Low-frequency isopycnal variability in the Alaska Gyre from Argo.\nProgress in Oceanography 168, 310–324 (2018).\nHristova, H. G., Ladd, C. & Stabeno, P. J. Variability and Trends of the Alaska Gyre from Argo and\nSatellite Altimetry. Journal of Geophysical Research Oceans 124, 5870–5887 (2019).\n• Could you better evaluate the relative contributions of basin-scale ocean dynamics versus surface\nflux forcing through simulations that involved forced ocean responses in a slab-ocean (mixed layer\nonly) versus a dynamically-responding ocean model?\nThank you for your comment and your suggestion, which we acknowledge would be a very nice\nseparate modelling study. However, we consider this to be beyond the scope of the current manuscript.\nIn lieu of this suggested modelling approach, we have instead substantially improved our diagnostic\nanalysis to strengthen the reliability of our conclusions. In the 'Mixed layer heat budget' section, we\nused a box model to decompose the temperature tendency into the terms representing the net surface\nheat flux, the net lateral advection solely, and the residual (Figure 3), to more completely explain the\nadvection contribution. The result shows the lateral advection term accounts for the positive\ntemperature tendency during the onset of the Blob and supports our hypothesis that the warm\nsubtropical water displacing into the typical subarctic region is the causal mechanism.\nWe included a more in-depth examination of the impact of net surface heat flux during the Blobs\nin the “Anomalous wind direction explains net heat fluxes” section, and found the decrease in SSTA\nduring the decay phase of the Blobs is primarily attributed to an increase in latent heat flux loss from\nthe ocean to the atmosphere, linked to the ocean thermodynamic response to the northwesterly surface\nwind anomaly that transports cold air into the Blob region.\nIn summary, I believe there is merit in the research directions this work reports on. However, I do not\nbelieve that the research reported here supports the conclusions that are offered. I recommend this\nmanuscript be returned to the authors with a request for major revisions.\nNate Mantua\nNOAA/NMFS Southwest Fisheries Science Center\nSanta Cruz, CA\nSpecific comments:\n40: I would add “and fisheries” after “marine ecosystems”\nThank you for your comment. We followed your advice by adding “and fisheries” after “marine\necosystems”.\n43-44: revise to “meridional expansion (contraction” of the subtropical (subarctic) gyre” … current\ntext doesn’t include subarctic\nThank you for your comment. This sentence is not revised because the northward excursion of the\ngyre boundary is the representation of the expansion of the subtropical gyre and contraction of the\nsubarctic gyre.\nThe contraction of the subarctic gyre is reported in the section of “Wind stress curl anomalies\ndrive gyre boundary poleward displacements”. Fig 5. (d) in the revision shows direct evidence of this\nvariation. The red shading over the upper 100m indicates an enhancement of the anticyclonic rotation\nover the northeast Pacific. If we add the background value, this anomaly value designates a contraction\nin the subarctic gyre and an expansion of the subtropical gyre.\nIn the revision, we have reduced the emphasis on changes in gyre boundaries and enhanced the\nfocus on the expansion and contraction of the two gyres.\n45-47: this analysis does not isolate the contribution to gyre scale frontal dynamics from other\nprocesses, nor does it connect those dynamics to the time scale of the events, so it falls short of\n“comprehensively explain(ing) the spatial characteristics and multi-year persistence of these …\nheatwaves”\nWe have deleted this sentence and instead followed your suggestions conducting more analysis on\nthe net heat flux and physical drivers, to make our analysis more comprehensive and conclusions more\ncompelling.\n60: reasons for multi-year persistence have been investigated and explained by others (e.g., DiLorenzo\nand Mantua 2016)\nThank you for your comments. We revised \"the reason for their different onset months and multi-\nyear persistence remains unclear\" to \" the reason for their different onset months and multi-year\npersistence remains the subject of further consideration and exploration.\" After \"multi-year\npersistence,\" we cited DiLorenzo and Mantua (2016).\n78-80: Isn’t a dominant tendency term required to understand the origins of a broad-scale warming of\nthe upper ocean that happens within a single season, hence the focus on atmosphere-ocean interactions?\nThank you for your comment. We don’t think so. Our original submission did not clarify this issue\nclearly, and we have now revised this section in the third paragraph of the Introduction in the revised\nmanuscript.\nIn particular, the heatwaves are linked to positive extremes in SSTA. Since the ocean and\natmosphere influence SSTA at different frequency bands or time scales during air-sea interactions, it\nis essential to first determine which frequency band or time scale is most influential in the changes in\nSSTA, and then analyze the dominant tendency term at that specific frequency band or time scale.\nCalculating temperature tendencies and identifying the main factors using SSTA filtered on an\nunsuitable timescale (or unfiltered) may lead to misunderstanding. This is because atmospheric\nchanges occur at a higher frequency than oceanic changes. The temperature tendencies resulting from\natmospheric processes are associated with high-value segments, while those from oceanic processes\nrelate to low-value segments. Analyzing unfiltered SSTA without considering the timescales that\nprimarily affect positive SSTA may skew the results towards higher values, which are coincidentally\nlinked to the atmospheric dynamics.\n83-86: This may be true in an ideal case, but this real-world case study does not involve a surface\nmixed layer having top and bottom interfaces that are adiabatic.\nThank you for your comment. We deleted this ideal case for it leads to confusion.\n120-125: Causation may go either way, or there may be a common forcing for the SSTa and latitude\nof the gyre boundary.\nThank you for your comment. The similarity between the latitude of the gyre boundary and SSTA\ncannot fully establish a causal relationship. However, the latitude of the gyre boundary corresponds to\nthe advection term in the mixed layer heat budget, and the mixed layer heat budget can help us\ndetermine that SSTA is indeed influenced by lateral advection.\nI also think that the trend over the 1994-2023 period of record is striking. Is there is a similar trend in\natmospheric forcing?\nYes, there is a corresponding trend in the atmospheric field. In the region where the correlation\ncoefficient between the gyre boundary and wind stress curl reaches a maximum (140-152 °W, 50-\n57 °N), we found a noticeable linear increase in the absolute values of the surface wind stress curl with\n-3.9×10-9 N/m3 per decade, and in the sea level pressure anomalies with a value of 79 Pa per decade\n(Figure R1-1; we also added it as Fig. 7b in our revised manuscript).\nWhy doesn’t the daily SSTA time series share the longer-term trend with the gyre boundary index?\nWe found that there is a linear trend primarily influenced by net heat flux in the net heat flux term\nanomaly, indicating a steady increase in ocean heat loss (Figure R1-4). Consequently, it is proposed\nthat the anomaly in temperature transport due to the northward expansion of the subtropical circulation\nis counterbalanced by the net heat flux anomaly, leading to a lack of a linear warming trend in local\nSSTA.\nFigure R1-4: Decomposition of the net surface heat flux term. The data used to draw the figure are\ñ𝑄 Q̃ 𝑄̅ℎ̃\nderived according to the equation: ( ) = − +𝑟𝑒𝑠𝑖𝑑𝑢𝑎𝑙 (Shi et al., 2022). A tilde denotes the\nh ℎ̅ ℎ̅2\nclimatological anomaly while an overline denotes the climatological daily mean. The climatological\ndaily mean was removed from the result and then a 365-day low-pass filter was applied. The data used\nin the following figure are from BRAN2020.\nReference\nShi, J. et al. Role of Mixed Layer Depth in the Location and Development of the Northeast Pacific\nWarm Blobs. Geophysical Research Letters 49, (2022).\n131-134: Why do you conclude that the subarctic gyre west of 150W expands southward? I can’t see\nthe southward expansion of the subarctic gyre west of 150 from the SLH or SST difference maps\n(panels e and f) – do you mean the negative SSH differences between 40-50N and 180-170W?. I don’t\nsee this in Fig 1b either (I would expect to see a dipole centered on ~150W, but that is only a feature\nof a few years like 2014-2016). The main changes in frontal position are east of 140W from 2014-\n2016, and east of 165W from 2018-2020 (based on Fig 1b).\nWe appreciate that you pointed out this point. We revised the text to ‘During the positive phase,\nthe subtropical gyre located east of 175 °W expands northward, leading to a southwest-northeast\noriented tilt of the gyre boundary (Fig. 2a). In contrast, during the negative phase, the gyre boundary\nshifts southward and assumes an almost zonal orientation (Fig. 2d)’.\nFigures 2 and 3: I think it would be best to use the same projection in the maps shown in Figures 2 and\n3.\nThank you for your comment. We now use the same map projection in the maps shown in Figures\n2 and 3 of the revised manuscript.\n168-170: Interesting that the zero curl contour at 180 is about 5 deg latitude farther south in your\npositive phase than negative phase composite, while it is at about the same latitude for both phases on\nthe west coast of N. America ~125W; This seems at odds with the SSH composites in Fig 2\nThank you for your comment. The zero curl contour and the latitude of the gyre boundary serve\nas indicators of changes in wind stress curl and ocean circulation. However, these two parameters do\nnot need to be exactly aligned in latitude, but rather in their orientation. The relationship between the\nocean circulation and wind stress curl can be confirmed by comparing Fig. 4c with Fig. 5a. The\nnegative anomaly in wind stress curl (shown in Fig. 4c) leads to a contraction of the subarctic gyre\n(indicated by blue shading in Fig. 5a and red shading in the subarctic region of Fig. 5d) and an\nexpansion of the subtropical gyre (represented by red shading in Fig. 5a and the subtropical region in\nFig. 5d). This variation occurs between latitudes 30-60 °N, with the gyre boundary being the area\nwhere these changes are most pronounced.\n176-179: Fig 3a shows your “positive phase” having the gyre boundary east of 150W shifting\nnorthward, not southward, in response to a negative wind stress curl anomaly in the Gulf of Alaska\n(north of 50N).\nWe appreciate that you pointed out this point. We revised ‘southward’ to ‘northward.’\n180-181: This is confusing. Perhaps you could include the actual wind stress curl anomalies for the\ncomposites of positive and negative phases and this would be easier for readers to see and understand\nthan trying to assess anomalies from the composites shown here.\nThank you for your comment. In revision, we added both the absolute wind stress curl and wind\nstress curl anomaly for the positive and negative phase, and wind stress curl anomaly in 2014 and 2019\nin Fig. 4.\nFigure 3 a-b-c and Fig 5 a-b: Wind vectors are really small and difficult to see, even at 200%\nmagnification. Please make them bigger.\nThank you for your comment. In revision, we made all the wind vectors bigger.\n211-212: Thermocline deepening is inconsistent with previously published findings (like Amaya et al\n2020 that show a shoaling mixed layer depth in summer 2019)?\nThank you for your comment. This is not contradictory; we are discussing the variations in the\nmain thermocline, primarily caused by geostrophic processes, while the key factors affecting the mixed\nlayer depth may differ. The mixed layer is influenced by changes in wind-induced mixing, surface net\nheat flux, stratification of the bottom of the mixed layer driven by changes of temperature and salinity,\nas well as advection processes.\nAlthough Amaya et al. (2020) noted a shallower mixed layer depth in the summer of 2019, Bond\net al. (2014) also observed a deeper mixed layer in winter. Additionally, Amaya's study area is located\nin the southeast of the Blob region. Their Figure 1c indicates that the shoaling of the mixed layer\nidentified in the area of 34°–47°N, 213°–232°W appears to be a northern extension of the deepening\nmixed layer found at lower latitudes (south of 34°N).\nScreenshot of Figure 2c. in Amaya, D. J. et al. Are Long-Term changes in mixed layer depth\ninfluencing North Pacific marine heatwaves? Bull. Amer. Meteor. Soc. 102, S59–S66 (2021).\n241-249: I don’t find this new definition of event onset to be very useful or convincing. Is there\nsomething special about climate dynamics in the month of May that makes you think this event onset\ncriteria leads to a new and better understanding for NEP marine heatwave events?\nThank you for your comment. In the revision, by carefully examining the 365-day low-passed\nSSTA, the neutral point for Blob 1.0 is identified as March 19, 2013, while the neutral point for Blob\n2.0 is determined to be May 6, 2018. Revisions have been made in accordance with these recent\nfindings.\nOur analysis of the emergence of Blob 1.0 seeks to address the question posed by Amaya et al.\n(2020) regarding the differences in timing between Blob 1.0 and Blob 2.0. They noted, \"An important\ndistinction between the recent temperature anomalies (referred to as Blob 2.0) and those from\n2013/2014 (referred to as Blob 1.0) is that Blob 2.0 primarily intensified during the summer months.\nIn contrast, Blob 1.0 developed in the winter due to a persistent atmospheric ridge in the Northeast\nPacific, which weakened the climatological Aleutian Low and associated surface winds.\"\nWhile defining the onset as October is reasonable from an event perspective, it leads to entirely\ndifferent conclusions when analyzing the mechanisms. A comparison of Fig. 3b and 3c reveals that the\nclassic mixed layer heat budget, after applying a 90-day low-pass filter, indicates that both ocean\nprocesses and surface net heat flux contribute to warming at the time marked by the black double-\nheaded arrow. In contrast, the results from a 365-day low-pass filter show that during the Blobs,\nocean processes consistently have a positive contribution, while surface net heat flux consistently has a\nnegative contribution (the direction of the flux is from ocean to atmosphere). This raises the question\nof which explanation to trust.\nThe primary fluctuations in SSTA in the Blob region occur on an interannual scale, suggesting that\nunderstanding the positive phase of the interannual component, rather than higher frequency variations\nof SSTA, can explain most of the positive anomalies. Therefore, we argue that focusing on the\ninterannual timescale, particularly by analyzing SSTA after a 365-day low-pass filter, is more suitable\nfor finding the casual mechanism.\nApparently, the climate dynamics differs in March and October. We did not examine this because\nthis is outside the scope of this manuscript.\nFigure R1-5 (Fig. 3 in the revision). The mixed layer heat budget. | (a) Mixed layer temperature\nanomalies (MLTA, ℃, red shading); also shown are the MLTA low-pass filtered by a 90-day cutoff\n(blue line) and MLTA low-pass filtered by a 365-day cutoff (black line). (b) Mixed layer temperature\ntendencies (℃/day) as a function of time, including net surface heat flux and integrated ocean process\ncontributions from the tendency budget, derived from daily resolution data, with the annual cycle\nremoved and 90-day lowpass filtered. (c) Same as (b) but 365-day lowpass filtered. (d) The heat budget\nfor a box model, and 365-day lowpass filtered. The magenta arrows denote the neutral point. The black\narrow denotes the cooling event in 2017. The black boxes in (c) denote the periods of ‘Blobs’. The\nblack boxes in (d) denote the period of the positive MLTA tendency that involved the onset of ‘Blobs’.\nWould you say that Blob 3.0 started in January 2022 by your criteria?\nWe are uncertain about this. According to our definition, Blob 3.0 appears to be an extension of\nBlob 2.0, as the interannual component of SSTA has not dropped below zero, indicating that there is\nno neutral point separating Blob 2.0 from Blob 3.0.\nDid the 2013-14 event end in 2015 or 2016 by your criteria?\nIt varies. The discrepancy in the timing (onset or decay) of the Blobs arises from the different foci\nof the disciplinary communities. The physical oceanographers’ concern is to explain how the positive\nSSTA was formed while the ecologists aim to understand how this above threshold SSTA affects the\nmarine ecosystem. Apparently, the timing of SSTA surpassing or dropping below the threshold works\nmore effectively for evaluating the biological impact of the temperature extreme, while the timing of\nSSTA exceeding or falling below the zero value works for examining the underlying physical\nmechanisms. These two definitions are not contradictory, as they reflect distinct perspectives from\ndisciplinary communities.\nThe SSTA dropped below the threshold during the winter of 2015, and it went below zero in the\nwinter of 2016. From a mechanistic perspective, we recommended the winter of 2016 as the endpoint,\nwhereas communities like ocean ecologists, who are more concerned with the impact of SSTA on the\necological environment, the winter of 2015 may be a better choice.\nBy and large, the neutral point of SSTA (from positive to negative phase, or from negative to\npositive phase) is a better choice when analyzing the physical mechanisms behind the marine heatwave\nevents.\n276-277: Did the marine heatwaves “result from” or “include” poleward expansion (contraction) of\nthe subtropical (subarctic) gyre?\nThe marine heatwaves are a consequence of the poleward movement (or shrinkage) of the\nsubtropical (or subarctic) gyre. The findings from the mixed layer heat budget attest to this. The mixed\nlayer heat budget analysis performed on the primary frequency band of SSTA indicates that oceanic\nadvection contributes significantly to the fluctuations in the SSTA.\n281-282: This analysis does not address the issue of “persistence”, but instead focuses on\n“onset”. What sets the time scale for these events? Is it the persistence and pattern of atmospheric\nforcing (along the lines of Di Lorenzo and Mantua’s analysis of the 2013-2015 period)? Or something\nelse?\nSee combined response to this question and the one below in our following response.\n295-297: If the basin-scale wind-stress curl is the key driver for the gyre variations, then multi-year\npersistence of these marine heatwaves is due to persistence in the wind-stress curl forcing. So this then\nbegs the question – what provides persistence in the basin-scale atmospheric forcing of the gyre\ncirculation? As noted, the ocean circulation (and mixed-layer heat content) response integrates\natmospheric forcing so it has greater variance at longer time scales than the atmospheric forcing, yet\nthe evolution of the atmospheric forcing has already been shown to play a key role in understanding\nthe evolution and persistence of the 2014-2016 (DiLorenzo and Mantua) and 2019-2020 (Amaya et al)\nevents.\nWe appreciate your motivating comment.\nThe responses to your two comments were combined into a single reply since they both address\nthe persistence of the Blob.\nWhat sets the time scale for these events? The physics behind the ocean dynamic response to the\nwind stress curl is essentially the Sverdrup relation, where the gyre-scale ocean dynamics favor the\nlow-frequency variation in the atmosphere. This is supported by the fact that the interannual changes\nin the gyre boundary are more significant than its seasonal variations (Fig. 1c).\nDi Lorenzo and Mantua (2016) noted in their Abstract, “… teleconnections between the North\nPacific and the weak 2014/2015 El Niño linked the atmospheric forcing patterns of this event. These\nteleconnection dynamics from the extratropics to the tropics during winter 2013/14, and then back to\nthe extratropics during winter 2014/15, are a key source of multi-year persistence of the North Pacific\natmosphere.” In their Fig. 4, the onset of the Blob is explained by the (atmospheric) North Pacific\nOscillation. However, from the perspective of the Sverdrup relation, changes driven by the NPO should\ninduce either equally strong or equally weak variations in the subtropical and subarctic gyre, rather\nthan be compensatory in their relationship. We conclude, therefore, that the NPO cannot be the\ndominant physical driver.\nInstead, the question becomes which atmospheric pattern drives the local atmospheric variation?\nWe analyzed the lag-correlation between the sea level pressure anomaly in the area where the wind\nstress curl anomaly has the highest correlation with the gyre boundary latitude and at each grid point\n(see Fig. 7) to identify the physical driver of the Blob. The results indicate that the Tropical Northern\nHemisphere pattern (as compared in Fig. 7h and the figure found at\nhttps://www.cpc.ncep.noaa.gov/data/teledoc/tnh_map.shtml) is likely the driving force behind the\nBlob. Interestingly, our results align with the conclusions of Liang et al. (2017), who suggested that\nthe Tropical Northern Hemisphere Pattern is closely linked to the formation of the Blob based on their\nstatistical analysis.\nMoreover, Mo and Livezey (1986) noted in Section 6 E, “The TNH pattern is principally\nconnected to interannual and longer time scales.” This offers additional support to our interpretation.\nOverall, the persistence of the Blob can be attributed to the ocean response to atmospheric\nconditions at the interannual timescale, with local interannual variations in the atmosphere\ncorresponding to the Tropical Northern Hemisphere Pattern (TNH), and the TNH features with\ninterannual variation. We added this explanation in our revised manuscript.\nReferences\nLiang, Y.-C., Yu, J.-Y., Saltzman, E. S. & Wang, F. Linking the tropical Northern hemisphere pattern\nto the Pacific warm blob and Atlantic cold blob. J. Clim. 30, 9041–9057 (2017).\nMo, K. C. & Livezey, R. E. Tropical-Extratropical Geopotential Height Teleconnections during the\nNorthern Hemisphere Winter. Monthly Weather Review 114, 2488–2515 (1986).\n324-328: Why would a ML heat budget fail to capture warm water advection into the “Blob”\nregion? Are the boxes drawn with boundaries that lie outside that of the warm water advection process?\nThank you for your comment. We have now hopefully better explained this in the Methods section\nof our revised manuscript.\nWe analyzed the complete form of the classic mixed layer heat budget equation to show why this\nequation cannot fully capture the advection process.\nThe complete equation of the classic mixed layer heat budget is:\n𝜕𝑀𝐿𝑇𝐴 𝑄 −𝑄 𝑇̅ −𝑇 𝜕𝐻\n= 𝑛𝑒𝑡 ℎ −𝒖̅ ∙𝛻𝑇̅ +𝑀𝑖𝑥𝑖𝑛𝑔−( −𝐻 )( +𝒖 ∙𝛻𝐻+𝑤 )\n𝜕𝑡 𝜌𝐶 𝐻 𝐻 𝜕𝑡 −𝐻 −𝐻\np\nThe second term on the right-hand side represents advection, while the fourth term represents\nentrainment. If we take the entrainment term apart, we can get,\n𝑇̅ −𝑇 𝜕𝐻 𝑇̅ 𝑇 𝑇̅ −𝑇 𝜕𝐻\n−𝐻 −𝐻 −𝐻\n−( )( +𝒖 ∙𝛻𝐻+𝑤 ) = − 𝑤 + 𝑤 −( )( +𝒖 ∙𝛻𝐻)\n𝐻 𝜕𝑡 −𝐻 −𝐻 𝐻 −𝐻 𝐻 −𝐻 𝐻 𝜕𝑡 −𝐻\nAccording to the continuity equation, the vertical velocity, 𝑤 can be written as 𝐻(𝛻 ∙𝒖 ),\n−𝐻\nwhich suggests the −\n𝑇̅\n𝑤 term can be written as −𝑇̅𝛻 ∙𝒖. This term represents the temperature\n−𝐻\nH\ntendency induced by the divergence or convergence of the horizontal velocity.\nMoreover, the entrainment term in the classic mixed layer heat budget is in the form of the\neffective flux, whereby the temperature difference (𝑇̅ −𝑇 ) multiplies vertical velocity (the term of\n−𝐻\n𝜕𝐻\n+𝒖 ∙𝛻𝐻 +𝑤 has a unit of m/s). Given that the temporal variation of the mixed layer depth\n−𝐻 −𝐻\n𝜕𝑡\n𝜕𝐻\n( ) may also be affected by the vertical motion, the entrainment term likely contains more\n𝜕𝑡\ncontribution from the lateral advection. Therefore, the advection term in this equation may\nunderestimate the contribution of the advection process.\nAnother simple way to understand this issue is, when we applied the classic mixed layer heat\n𝜕𝐻\nbudget equation to a box (the bottom of the box is flat and does not vary with time), the and 𝛻𝐻\n𝜕𝑡\nterms vanished. By moving the net heat flux term to the left-hand side, we get:\n𝜕𝑀𝐿𝑇𝐴 𝑄 −𝑄 𝑇̅ −𝑇\n− 𝑛𝑒𝑡 ℎ = −𝒖̅ ∙𝛻𝑇̅ −( −𝐻 )𝑤 +𝑀𝑖𝑥𝑖𝑛𝑔\n𝜕𝑡 𝜌𝐶 𝐻 𝐻 −𝐻\np\nExpanding the second term in the right-hand side and rearranging, we have:\n𝜕𝑀𝐿𝑇𝐴 𝑄 −𝑄 𝑇̅ 𝑇\n− 𝑛𝑒𝑡 ℎ = −𝒖̅ ∙𝛻𝑇̅ − 𝑤 + −𝐻 𝑤 +𝑀𝑖𝑥𝑖𝑛𝑔\n𝜕𝑡 𝜌𝐶 𝐻 𝐻 −𝐻 𝐻 −𝐻\np\nwhere −𝒖̅ ∙𝛻𝑇̅ − 𝑇̅ 𝑤 is the lateral advection and 𝑇−𝐻𝑤 represents the vertical heat\n−𝐻 −𝐻\nH H\nadvection.\nA simple scale analysis could help to evaluate the relative importance of each term. In the Blob\nregion, the average horizontal velocity at the sea surface is 0.05 m/s, the maximum horizontal\ntemperature gradient is 9*10-6 ℃/m (1℃/110km* 0.001 km/m), the vertical velocity in the upper ocean\nis 10-4 m/s, the annual mean mixed layer depth is 50 m, SST = 12 ℃, 𝑇 = 12 ℃-0.3℃=11.7℃ (0.3℃\n−𝐻\nis the common criteria for defining the mixed layer depth), and suppose 𝑇̅ = 𝑆𝑆𝑇+𝑇−𝐻 = 11.85℃.\n2\nOn the right-hand side of the above equation, the first term (𝒖̅ ∙𝛻𝑇̅) is on the order of 4.5*10-\n7 ℃/s, the second term ( 𝑇̅ 𝑤 ) is on the order of 2.37*10-5 ℃/s, and the third term ( 𝑇−𝐻𝑤 ) is on the\n−𝐻 −𝐻\nH H\norder of 2.34*10-5 ℃/s.\n𝑇̅\nTherefore, both the lateral advection ( 𝑤 ) and the vertical advection at the bottom of the mixed\n−𝐻\nH\nlayer (\n𝑇−𝐻𝑤\n) plays an important role in the heat budget equation. However, if we write them in the\n−𝐻\nH\nform of the effective flux (−\n𝑇̅\n𝑤 +\n𝑇−𝐻𝑤\n, temperature difference multiplying vertical velocity),\n−𝐻 −𝐻\nH H\nthe two terms cancel out and generate a small term that is on the order of -3*10-7 ℃/s (-2.37*10-\n5 ℃/s+2.34*10-5 ℃/s), which is comparable to the term of 𝒖̅ ∙𝛻𝑇̅. This analysis suggests that if the\npredominance of a particular term is assessed based on its magnitude, both the vertical and lateral\nprocesses may be undervalued when the entrainment term is represented as the effective flux.\nWouldn’t frontal shift-induced SSTA be maximized along the axis of the SST front, and zonally\nconfined to the SST front?\nThe variations in SSTA closely align with the ocean circulation (as shown in the plan view: Fig.\n2g and h, and sectional view: Fig. 5d with color shading and black line). The gyre boundary\ndisplacements serve as a crucial indicator of changes in ocean circulation. This boundary aligns with\nthe SST front, meaning it intersects the area of maximum SSTA (Fig. 2h). However, both the sea level\nand geostrophic stream function suggest that the changes are not limited to the frontal zone; rather,\nthey extend spatially from 30 to 60 degrees North. Consequently, the SSTA caused by the shift in the\nfront is not restricted to the SST front in a zonal manner.\n342-343: Is the 30 deg latitude extent of the SSTA warm pattern consistent with the 200km (2 deg\nlatitude) gyre boundary excursion?\nThank you for reminding us. We examined the thermal displacement and found it increased with\nlatitude (Fig. R1-4). Therefore, we revised this statement to “In this case, the northward thermal\ndisplacement caused by the gyre boundary excursion increases with latitude from less than 100\nkilometres at about 35 °N and over 200 kilometres at about 52 °N.”\nFigure R1-6: Comparison of isotherms of SST during the positive phase (black line) and negative\nphase (blue line).\nWhat about the scale of atmospheric forcing anomalies associated with surface heat flux anomalies?\nFocusing on yearly changes in the surface heat flux anomaly (Fig. R1-7), in 2013, the spatial\ndistribution of the annual average surface heat flux anomaly was concentrated in the northeast Pacific\nregion, whereas in other years, the spatial pattern was less uniform.\nFigure R1-7: The surface heat flux anomaly referring to an annual mean over 1993 to 2022. A\npositive value indicates that the ocean is absorbing heat.\nReviewer #2 (Remarks to the Author):\nThis study tried to explore the potential role of gyre scale ocean dynamics in iconic northeast Pacific\nmarine heatwaves, focusing on the interannual or low-frequency timescales. I think this viewpoint is\nvery important and potentially crucial for understanding the climate variability in the Northeast Pacific.\nBut I think the current presentation is quite preliminary in many aspects. Please see my detailed\ncomments below.\nWe appreciate your acknowledgment of the importance of this work. Below, you will find our\ndetailed responses to your comments, with your comments in black font and our replies in blue font.\nDue to significant changes in the manuscript compared to the previous submission, we want to\nsummarize the main points in the revision. The “Introduction” has been completely rewritten to\nprovide a more comprehensive explanation of our concerns about the conventional mixed layer heat\nbudget analysis. The section headed 'Mixed layer heat budget' has been restructured to follow the 'Gyre\nboundary variation in the 'Blob' region', highlighting the importance of advection in the heat budget.\nFurthermore, we have added two new sections: one examining the variations in net heat flux and the\nother addressing the physical drivers of the Blobs. The final section has also been revised to incorporate\nour new findings.\nMajor comments:\nIn your mixed-layer heat budget equation, the H is time varying mixed layer depth. Please specify how\nyou calculate the surface heat flux term, that is Qnet/(ρCpH). In particular, how does the mixed layer\ndepth vary in this term? Moreover, I cannot tell exactly how do the gyre variations work on the MLTA\nin the blob region, because authors merged the oceanic processes into one term. The physical picture\nis not clear to me. Moreover, I cannot tell exactly how do the gyre variations work on the MLTA in the\nblob region, because authors merged the oceanic processes into one term.\nThank you for your comment. In our revision, the surface heat flux term is calculated according\nto\n𝑄𝑛𝑒𝑡−𝑄ℎ,\nwhere 𝑄 is the surface net heat flux, calculated by the sum of the latent and sensible\n𝑛𝑒𝑡\n𝜌𝐶𝑝𝐻\nheat flux, the longwave and short wave radication, 𝑄 is the penetrative shortwave radiation, H is the\nh\ntime-varying mixed layer depth, and other factors are calculated as follows, 𝜌𝐶 = 4.088×\np\nh h\n106 J℃−1m−3 and 𝑄 ℎ = 𝑄 𝑠𝑤 (0.58𝑒 − 0.35 +0.42𝑒− 23) (Shi et al., 2022). All the data are daily\nresolved.\nSince the relationship between mixed layer depth and net surface heat flux is not our main focus\nand has been thoroughly addressed by Shi et al. (2022), and cited in our manuscript, we did not\nexamine its variations. The importance of advection is highlighted by its significant role in the positive\nMixed Layer Temperature Anomaly (MLTA) during the onset of the Blobs, as shown in Fig. 3c. We\nemployed a box model to analyze the temperature tendency. The rationale for not utilizing the\ntraditional mixed layer heat budget is explained in response to the next comment. The lateral advection\nterm in our box model effectively captures the contributions from the oceanic advection process (Fig.\n3d). This advection corresponds to the Sverdrup transport, which is fundamental to the wind-driven\nocean general circulation. Fig. 4 demonstrates that slow variations in the large-scale ocean circulation\ncan be explained by this theory, while Fig. 5 supports this by analyzing the geostrophic stream function\nin a sectional view. This reason was part of our previous submission, although the figures mentioned\nare arranged differently in the revision.\nAdditionally, Fig. 3c indicates that the net heat flux may contribute to the decay of the Blob.\nConsequently, we investigated the surface heat flux and found that the latent heat flux is the most\nsignificant factor (Fig. 6a). Further analysis revealed that the unusual wind direction affects the 2-\nmeter air humidity, which in turn influences the latent heat flux. The persistence of the Blob is briefly\ndiscussed through the lag-correlation between sea level pressure anomalies in the region where the\nwind stress curl anomaly shows the largest correlation coefficient with the gyre boundary latitude and\nat each geographical point. The result suggests the persistence can be attributed to a Tropical Northern\nHemisphere Pattern-like teleconnection pattern (Fig. 7).\nReferences\nShi, J. et al. Role of mixed layer depth in the location and development of the Northeast Pacific warm\nblobs. Geophys. Res. Lett. 49, e2022GL098849 (2022).\nIn many previous studies, as you also shown, the oceanic processes may not be that important for the\nwarm blobs. How to interpret your results?\nThank you for your comment. In previous studies, oceanic process contributions were found to be\nnegligible, with conclusions derived from diagnostic analysis of the mixed layer heat budget. However,\nwe find that the classic mixed layer heat budget misleads these interpretations.\nIn the present study, we analyzed the more complete mixed layer heat budget equation, to show\nwhy the advection term in this equation cannot fully capture the advection process.\nThe complete equation of the classic mixed layer heat budget is:\n𝜕𝑀𝐿𝑇𝐴 𝑄 −𝑄 𝑇̅ −𝑇 𝜕𝐻\n= 𝑛𝑒𝑡 ℎ −𝒖̅ ∙𝛻𝑇̅ +𝑀𝑖𝑥𝑖𝑛𝑔−( −𝐻 )( +𝒖 ∙𝛻𝐻+𝑤 )\n𝜕𝑡 𝜌𝐶 𝐻 𝐻 𝜕𝑡 −𝐻 −𝐻\np\nThe second term on the right-hand side represents advection, while the fourth term represents\nentrainment. If we take the entrainment term apart, we get\n𝑇̅ −𝑇 𝜕𝐻 𝑇̅ 𝑇 𝑇̅ −𝑇 𝜕𝐻\n−𝐻 −𝐻 −𝐻\n−( )( +𝒖 ∙𝛻𝐻+𝑤 ) = − 𝑤 + 𝑤 −( )( +𝒖 ∙𝛻𝐻)\n𝐻 𝜕𝑡 −𝐻 −𝐻 𝐻 −𝐻 𝐻 −𝐻 𝐻 𝜕𝑡 −𝐻\nAccording to the continuity equation, the vertical velocity, 𝑤 can be written as 𝐻(𝛻 ∙𝒖 ),\n−𝐻\nwhich suggests, the −\n𝑇̅\n𝑤 term can be written as −𝑇̅𝛻 ∙𝒖. This term represents the temperature\n−𝐻\nH\ntendency induced by the divergence or convergence of the horizontal velocity. Moreover, the\n𝑇−𝐻𝑤\n−𝐻\nH\nterm represents the net vertical heat advection through the bottom of the mixed layer.\nThe entrainment term in the classic mixed layer heat budget is in the form of the effective flux,\nwhereby the temperature difference (𝑇̅ −𝑇 ) multiplies vertical velocity (the term of 𝜕𝐻 +𝒖 ∙\n−𝐻 −𝐻\n𝜕𝑡\n𝜕𝐻\n𝛻𝐻+𝑤 has a unit of m/s). Given that the temporal variation of the mixed layer depth ( ) may\n−𝐻\n𝜕𝑡\nalso be affected by the vertical motion, the entrainment term likely contains more contribution from\nthe lateral advection. Therefore, the advection term in this equation may underestimate the contribution\nof the advection process.\nAnother simple way to understand this is when we applied the classic mixed layer heat budget\n𝜕𝐻\nequation to a box (the bottom of the box is not time-varying and is flat), the and 𝛻𝐻 terms\n𝜕𝑡\nvanished. By moving the net heat flux term to the left-hand side, we get:\n𝜕𝑀𝐿𝑇𝐴 𝑄 −𝑄 𝑇̅ −𝑇\n− 𝑛𝑒𝑡 ℎ = −𝒖̅ ∙𝛻𝑇̅ −( −𝐻 )𝑤 +𝑀𝑖𝑥𝑖𝑛𝑔\n𝜕𝑡 𝜌𝐶 𝐻 𝐻 −𝐻\np\nExpanding the second term on the right-hand side and rearranging, we have:\n𝜕𝑀𝐿𝑇𝐴 𝑄 −𝑄 𝑇̅ 𝑇\n− 𝑛𝑒𝑡 ℎ = −𝒖̅ ∙𝛻𝑇̅ − 𝑤 + −𝐻 𝑤 +𝑀𝑖𝑥𝑖𝑛𝑔\n𝜕𝑡 𝜌𝐶 𝐻 𝐻 −𝐻 𝐻 −𝐻\np\nwhere −𝒖̅ ∙𝛻𝑇̅ − 𝑇̅ 𝑤 is the lateral advection and 𝑇−𝐻𝑤 represents the vertical heat\n−𝐻 −𝐻\nH H\nadvection.\nA simple scale analysis could help to evaluate the relative importance of each term. In the Blob\nregion, the average horizontal velocity at the sea surface is 0.05 m/s, the maximum horizontal\ntemperature gradient is 9*10-6 ℃/m (1℃/110km* 0.001 km/m), the vertical velocity in the upper ocean\nis 10-4 m/s, the annual mean mixed layer depth is 50 m, SST = 12 ℃, 𝑇 = 12 ℃-0.3℃=11.7℃ (0.3℃\n−𝐻\nis the common criteria for defining the mixed layer depth), and suppose 𝑇̅ = 𝑆𝑆𝑇+𝑇−𝐻 = 11.85℃.\n2\nOn the right-hand side of the above equation, the first term (𝒖̅ ∙𝛻𝑇̅) is on the order of 4.5*10-\n7 ℃/s, the second term ( 𝑇̅ 𝑤 ) is on the order of 2.37*10-5 ℃/s, and the third term ( 𝑇−𝐻𝑤 ) is on the\n−𝐻 −𝐻\nH H\norder of 2.34*10-5 ℃/s.\n𝑇̅\nTherefore, both the lateral advection ( 𝑤 ) and the vertical advection at the bottom of the mixed\n−𝐻\nH\nlayer (\n𝑇−𝐻𝑤\n) plays an important role in the heat budget equation. However, if we write them in the\n−𝐻\nH\nform of the effective flux (−\n𝑇̅\n𝑤 +\n𝑇−𝐻𝑤\n, temperature difference multiplying vertical velocity),\n−𝐻 −𝐻\nH H\nthe two terms cancel out and generate a small term that is on the order of -3*10-7 ℃/s (-2.37*10-\n5 ℃/s+2.34*10-5 ℃/s), which is comparable to the term of 𝒖̅ ∙𝛻𝑇̅. This analysis suggests that if the\npredominance of a particular term is assessed based on its magnitude, both the vertical and lateral\nprocesses may be undervalued when the entrainment term is represented as the effective flux.\nThe mixed layer heat budget was not filtered to extract the most influential SSTA component in\nprevious study. Therefore, we assessed the heat budget within the frequency band that primarily\ninfluences the positive SSTA. Despite the classic mixed layer heat budget facing challenges in breaking\ndown the oceanic term, we consider its most basic form to be trustworthy. This indicates that the mixed\nlayer temperature anomaly (MLTA), the surface heat flux term, and the oceanic term—defined as the\ndifference between the two—are all reliable. Our analysis reveals that during the Blob period, the\noceanic term is positive, whereas the surface heat flux term is negative (the direction of the flux is\nfrom ocean to atmosphere). Consequently, we conclude that the oceanic advection plays a key role in\nthe formation of the Blob, while the surface heat flux is more important for its decay.\nIn this study, authors did not pay special attention to the differences of mechanisms between seasons.\nThank you for your comment. This is not our primary concern. Our analysis focuses on the\ninterannual timescale, and hence here seasonal variations are removed.\nWhat is the frontal dynamics? I think authors did not use too much frontal theory in this study.\nThank you for your comment. The process responsible for the northward movement of the front\nis based on classic wind-driven ocean circulation theory provided by the Sverdrup relation. The\ncurl\nof the wind stress causes downward Ekman pumping and stretching of the water column. This results\nin a northward movement of the water column to maintain potential vorticity conservation (D((f+ζ)/H)\n/Dt=0, where f and ζ represent the planetary and relative vorticity, and H denotes the thickness of the\nwater column). Based on the premise of geostrophic motion, the Rossby number << 1, or ζ<<f, the\npotential vorticity conservation represents the balance between planetary vorticity (f) and the layer\nthickness (H). More specifically, it is the shifts in the temperature front's position that are significant,\nrather than the front itself (for example, its intensity). Thus, it is not necessary to apply frontal theory\nin this context.\nFigure 2: I think there should be a more objective method to select northernmost and southernmost\nyears. Then, significance test should be added to confirm the robustness of your results. Currently, the\nsample numbers are small. Moreover, authors can add the climatological +0.47 line to show the\nchanges. Same issue for Fig. 3.\nThank you for your comment. The composite map provides a readily comprehensible overview of\nthe changes in ocean circulation and wind patterns. A significance test (p-value) was performed to\nshow the statistical significance of the correlation coefficient between the gyre boundary and SSTA\n(Fig. 2i), as well as between the gyre boundary and wind stress curl (Fig. 7a). Therefore, our\nconclusions are not influenced by the selection of the northernmost and southernmost years.\nSince the composite map is based on absolute values rather than anomalies, the years chosen for\nthe composite were selected based on specific criteria: first, we used data from entire years for\naveraging to avoid seasonal cycle aliasing; second, we noted a linear trend indicating a northward shift\nof the gyre boundary, so we aimed to select years that are evenly spread across each decade to prevent\nthe composite map from reflecting a long-term trend; third, we created anomalies for 2014 and 2019\nand compared them with the positive phase to ensure that the spatial patterns were similar.\nWe included the +0.47 contour line in each panel of Figure 2, and replaced the 0.47 contour line\nwith the zero contour line of the stream function derived from the wind stress curl.\nCurrently, I cannot understand well the texts in Line 135-137.\nWe appreciate you mention it. The description in lines 135-137 is problematic. Reviewer #1 also\nnoted this. We have now revised the text to “During the positive phase, the subtropical gyre located\neast of 175 °W expands northward, leading to a southwest-northeast oriented tilt of the gyre boundary\n(Fig. 2a). In contrast, during the negative phase, the gyre boundary shifts southward and assumes an\nalmost zonal orientation (Fig. 2d).”\nLine 215-218: Previous studies claimed the shallower mixed layer depth during the Northeast Pacific\nwarm blobs/marine heatwaves (Amaya et al. 2020 NC; Shi et al. 2022 GRL). Could authors give more\nexplanations on the deepening of thermocline here.\nThank you for your comment. In practical terms, if the subtropical gyre expands northward and\nencroaches upon the area occupied by the subarctic gyre, the existing characteristics of low sea level,\na shallower thermocline, and cooler sea surface temperatures, will be replaced by those of the\nsubtropical gyre, which include higher sea levels, a deeper thermocline, and warmer SST.\nFrom a theoretical standpoint, the principle of potential vorticity (PV) conservation explains why\nthe northward movement of the water column is typically linked to a deepening of the main\nthermocline. The equation representing PV conservation is D((f+ζ)/H) /Dt=0, where f is the Coriolis\nparameter and H is the water column thickness. Based on the premise of geostrophic motion, the\nRossby number << 1, or ζ<<f, the potential vorticity conservation represents the balance between\nplanetary vorticity (f) and the layer thickness (H). As the water column moves north, it enters a higher\nlatitude region where f increases. To maintain PV conservation, the water column must also be\nstretched, leading to the increase of H. Conversely, as the water column stretches first (H increases), it\nalso moves northward (f increases) subsequently to satisfy this relationship. In this case, the wind stress\ncurl causes Ekman pumping, which stretches the water column before it moves northward, contributing\nto the Sverdrup balance.\nThe main thermocline acts as the boundary between the upper ocean and the deep ocean. The\nstretching of the water column in the upper ocean causes the main thermocline to descend, which is\nthe dynamical response outlined here.\nThe neutral point appears a bit arbitrary. Authors tried to raise a new method to identify the onset of\nthe blob, but they should compare in details with the previous results. Why this method is more\nreasonable than others?\nThank you for your comment. In the revision, by carefully examining the 365-day low-passed\nSSTA, the neutral point for Blob 1.0 is identified as March 19, 2013, while the neutral point for Blob\n2.0 is determined to be May 6, 2018. Revisions have been made in accordance with these findings.\nOur analysis of the emergence of Blob 1.0 seeks to address the question posed by Amaya et al.\n(2020) regarding the differences in timing between Blob 1.0 and Blob 2.0. They noted, \"An important\ndistinction between the recent temperature anomalies (referred to as Blob 2.0) and those from\n2013/2014 (referred to as Blob 1.0) is that Blob 2.0 primarily intensified during the summer months.\nIn contrast, Blob 1.0 developed in the winter due to a persistent atmospheric ridge in the Northeast\nPacific, which weakened the climatological Aleutian Low and associated surface winds.\"\nWhile defining the onset as October is reasonable from an event perspective, it leads to entirely\ndifferent conclusions when analyzing mechanisms. A comparison of Fig. 2b and 2c reveals that the\nclassic mixed layer heat budget, after applying a 90-day low-pass filter, indicates that both ocean\nprocesses and surface net heat flux contribute to warming at the time marked by the black double-\nheaded arrow. In contrast, the results from a 365-day low-pass filter show that during the Blobs,\nocean processes consistently have a positive contribution, while surface net heat flux consistently has a\nnegative contribution (the direction of the flux is from ocean to atmosphere). This raises the question\nof which explanation to trust.\nThe primary fluctuations in SSTA in the Blob region occur at the interannual timescale, suggesting\nthat understanding the positive phase of the interannual component of SSTA can explain most of the\npositive anomalies. Therefore, we argue that focusing on the interannual timescale, particularly by\nanalyzing SSTA after a 365-day low-pass filter, is more suitable.\nMoreover, did Amaya et al. (2020) propose the onset of blob 2.0 in May 2018? I think they only\nfocused on the summer of 2019, rather than 2018.\nAmaya et al. (2020) proposed the Blob 2.0 was intensified in May 2019, and did not analyze (or\ndetermine) the onset of Blob 2.0. We have corrected the text.\nMinor comments:\nThe word “Iconic” is a bit hard to clearly understand. Could authors change into another word?\nThank you for your comment. We deleted this word.\nAbstract, Line 43: anomalous positive or negative wind stress curl drives the expansion of the\nsubtropical gyre. Please specify here.\nThank you for your comment. We revised this sentence to “The northward expansion (or\ncontraction) of the subtropical (subarctic) gyre is influenced by anomalous wind stress curl and can\nbe understood through wind-driven ocean circulation theory, specifically the Sverdrup relation.”\nBlob 1.0 and 2.0: It is not appropriate to appear in the abstract without any background here. Readers\nmay not know their similarities and differences.\nThank you for your comment. We deleted this sentence.\nAbstract, Line 52: It is better to specify the sea surface temperature, rather than temperature.\nThank you for your comment. This sentence is deleted in the revision.\nLine 61: Chen et al. (2021) classified two types of warm blobs and reported different mechanisms of\nthe warm blobs in different seasons.\nChen, Z., Shi, J., & Li, C. (2021). Two types of warm blobs in the Northeast Pacific and their potential\neffect on the El Niño. International Journal of Climatology, 41(4), 2810–2827.\nThank you for your comment. We cited this paper in the Line 58 in the revision.\nLine 64: For Blob 1.0, due to its multi-year persistence over 2013-2015, the growth months may be\ndifferent in these years.\nThank you for your comment. We removed \"growth\" since the two time points are centered on the\nonset of the Blobs.\nThere are some errors in the reference list, such as the missing page information, etc.\nThank you for your comment. We have now carefully checked and revised this list.\nLine 85-86: Why the advection is zeros here? Could you please draw a simple diagram to illustrate?\nMoreover, why do you give this ideal case? I feel confused to understand. How does advection outside\nthe computational domain explain the temperature tendency in the domain?\nThank you for your comment. We have now deleted this ideal case as it leads to misunderstanding.\nHere we explain the idea of advection that occurs outside the computational domain, a concept\nfirst introduced by Lee et al. (2004) and termed \"external advection.\" They incorporated the advection\nterm into the classic mixed layer heat budget, where the dot product of the temperature gradient and\nhorizontal velocity signifies heat redistribution within the computational domain (referred to as\ninternal advection). However, internal advection alone does not explain the overall movement of a\nwater parcel within the computational domain (control volume).\nThis issue arises from the distinction between the Eulerian and Lagrangian perspectives (Figure\nR2-1). From the Eulerian viewpoint, the computational domain used to calculate the net heat flux term\nremains fixed, while the water column beneath it is in motion. From the Lagrangian viewpoint, the\nwater column within the specified area (shown by the black parallelogram, A1) moves from one time\npoint (T1) to another (T2) due to advection, arriving at the position indicated by the red dashed\nparallelogram (A2).\nWe will assume that there is a transient net heat flux impacting the time point from T1 to T2. This\nnet heat flux leads to temperature changes between the area represented by the black parallelogram at\nT1A1 and the area depicted by the red dashed parallelogram at T2A2. In the budget calculation, the\ntemperature change is evaluated based on the temperatures in the black parallelogram at both T1 and\nT2, i.e., T1A1 and T2A1. The temperature difference utilized to calculate the temperature tendency\nshould accurately be represented as T2A2 - T1A1; however, the value that is computed is T2A1 - T1A1.\nClearly, the results obtained from the two methods differ. This discrepancy results from the overall\nmovement of the water parcel within the computational domain, and is termed external advection.\nAdditionally, internal advection refers to the movement of a single particle from one point in the\ncomputational domain to another.\nNext, we can evaluate how much this external advection influences the results. Typically, the\nclassic mixed layer heat budget is analyzed on a monthly basis. In the Blob region, the average\nhorizontal velocity at the sea surface is 0.05 m/s, which translates to a movement of 130 km over the\ncourse of a month (0.05 m/s * 0.001 km/m * 86400 s/day * 30 days/month). This suggests that nearly\n10% of the water column is displaced beyond the computational domain, and this effect is not\naccounted for in the classic mixed layer heat budget.\nFigure R2-1: Schematic view of the external advection. The black parallelogram (A1) indicates the\ncomputational domain and the red dashed parallelogram (A2) indicates the position where the water\ncolumn arrived at.\nLine 98-99: Please cite some related papers to confirm this.\nThank you", "reviewer_2": "I thank the authors for their great efforts in addressing the comments and concerns raised by the three reviewers. Although some of my comments have been well addressed, I may have more concerns to put forward in this round of review:\n1. The motivation or evidence of authors to focus on the interannual timescales is not clear, although interannual time scale/1-year filter is used in many places throughout the manuscript.\n2. In the revision, authors used a box model to try to give calculation for the heat budget of the blob marine heatwaves. However, based on the method description from Line 620, this method may not well distinguish horizontal and vertical advection processes. On the other hand, the relationship between the mixed layer heat budget and box model heat budget is not well addressed, which may affect the consistency and readability of this paper. In the following subsection b), authors give detailed discussion on the entrainment term, whose relationship with the box model also needs clear description.\n3. Related to comment 2, the response last time mainly addressed from the perspective of vertical processes, however, I think authors tried to explain the warm water transportation from the view of horizontal advection. But authors did not show any advantages of the box model from the view of horizontal advection.\n4. As other reviewers also pointed out the concern on the relationship between the shoaling mixed layer depth anomalies and deepened thermocline depth anomalies, I think authors should give the clear figures (spatial pattern and temporal evolution) to better clarify this issue.\n5. In terms of the circulation pattern, authors attributed to the TNH, rather than NPO. I think authors should note that the variance of TNH, based on the EOF analysis, is only around 5%, which may not well explain the frequent occurrence of marine heatwaves over the past decade. Moreover, the existence of TNH in the abstract as the main findings may not be good as they have been already reported by other studies.\n6. In the Fig. R1-3, I cannot see any obvious trend, although authors have processed the data.\nAdditional comments:\n1. Tracked-changes version: I cannot see any tracked changes. Please make sure that authors upload the correct file.\n2. Affiliation 7 did not belong to any author. Please check.\n3. Fig. 5e,f: Why did authors show the magenta lines here?\n4. Although authors added some statistical significance tests in the revision manuscript, there are other places that should have such test, for example, Figs. 2g, h.", "reviewer_3": "Most of the points I previously mentioned have been appropriately revised, so I suggest acceptance once the minor comments are addressed.\nMinor comments\nLine 91: In the phrase \"geostrophy influences ~\", geostrophy itself cannot induce vertical motion. However, vertical motion can occur when baroclinicity and the Ekman effect act together. This sentence may lead to misunderstanding, so I suggest revising it to explicitly include this explanation.\nFigure 2: Changing the colorbar for the Difference plots to a different color would make the distinction between (a)-(f) and (g)-(i) clearer."}, "author_response": "Thank you very much for your thoughtful review. We have now revised our manuscript based on your comments and those of the other two reviewers. Our detailed responses are provided below, with your original comments in black font and our replies in blue font."}, {"round": 3, "reviewer_comments": {"reviewer_2": "Thanks for authors to carefully address my concerns. I am mostly clear about my concerns. But I still have some minor comments before the manuscript being potentially accepted.\n1.Figure R2-1 can be moved to supplementary information, but now I cannot find any supplementary information.\n2.In the last paragraph of response to my second comment, it would be great if authors could add figures to show each term and their corresponding correlations.\n3.In the last paragraph of response to my third comment, I think authors could re-write it especially the last sentence and try to use the same marks (including consistent upper or lower cases with previous comment). Why is the first term on the right hand side derived from six boundaries of the box?\n4.I think the deepening of thermocline should be illustrated explicitly. Please show the figures.\nAs I was asked to take a look at the points raised by Reviewer #1 as well, I will also include those comments as below:\n1. I suggest authors to give specific line numbers for their modifications, which is not clear now.\n2. In terms of the response to #1 comment, authors did not prove the causality between the warming trend and displacement of the gyre boundary. Thus, “a major contributor …” needs to be rephrased.\n3. \"For the multi-season persistence, shifts in the Sverdrup balance are favored by the low frequency component of the atmospheric forcing, with the ocean acting as a low-pass filter to the atmospheric forcing\": By saying these, are ocean dynamics key to the multi-season persistence? I am not clear. Would it be possible to quantify the contribution to the SSTAs of the MHWs?\n4. \"ocean advection is the only factor contributing to the emergence of the two warm events\": How to reconcile with previous studies? I think this issue should be better elaborated in the manuscript."}, "author_response": "Thank you very much for your thoughtful review. We have now revised our manuscript based on your comments. Our detailed responses are provided below, with your original comments in black font and our replies in blue font."}, {"round": 4, "reviewer_comments": {"reviewer_2": "Thanks for addressing my concerns and comments previously. I am now glad to accept this paper."}, "author_response": ""}]

Decision letter