An Optimal Precursor of Northeast Pacific Marine Heatwaves and Central Pacific El Niño Events

The intensity of Northeast Pacific marine heatwaves (MHWs) has been related to local stochastic atmospheric forcing with limited predictability, but their evolution and persistence may be controlled by large‐scale climate influences. A Linear Inverse Model containing both sea surface temperature and sea surface height (SSH) anomalies is used to identify the “optimal” conditions for observed Northeast Pacific MHW events that developed two‐to‐four seasons later. These optimal initial conditions include SSH anomalies that are responsible for most of the MHW growth, suggesting the relevance of subsurface ocean dynamics. Moreover, Northeast Pacific MHW growth occurs as part of a basin‐scale dynamical mode that links the North Pacific to central equatorial Pacific El Niño events, whose subsequent development may lengthen MHW duration.

be more sensitive to equatorial SST anomalies (SSTAs) in the central-western Pacific (Capotondi, Sardeshmukh, et al., 2019), highlighting the importance of accounting for this diversity of ENSO expressions when examining the ENSO-MHW connection.
Others (Amaya et al., 2021;Bond et al., 2015) have emphasized the role of persistent local atmospheric forcing, together with enhanced upper-ocean stratification and reduced mixed-layer depth as causes of the extreme intensity and persistence of Northeast Pacific MHWs. Recent results (Xu et al., 2021;XNCD hereafter) show that the extreme warm conditions that have recently occurred in the Northeast Pacific are part of a continuum of events with a broad range of intensities/durations (Scannel et al., 2016;XNCD), and with similar temporal evolutions. The duration of these events appears to be controlled by the timing of ENSO development in the equatorial Pacific, while the intensity is primarily influenced by local atmospheric forcing. Similar evolution of these events suggests that large-scale climate conditions may help drive the development of Northeast Pacific MHWs. In this study, we thus build on XNCD to: (a) Ascertain that such large-scale conditions do exist and can provide some degree of MHW predictability; (b) Elucidate their dynamics, and the relative role of North and tropical Pacific influences; and (c) Examine the role of local and remote oceanic processes in the Northeast Pacific MHW onset.
The paper is organized as follows: In Section 2 we describe the data and methodology used, while in Section 3 we examine the composite evolution of northeast Pacific MHWs, and determine the large-scale precursors of MHWs. Conclusions are drawn in Section 4.

The Linear Inverse Model
The time evolution of a climate state x may often be approximated by the stochastically forced linear dynamical system (Penland & Sardeshmukh, 1995): where L is a dynamical operator encapsulating the linear dynamics of the system, including an implicit linear parameterization of "fast" (decorrelation timescales much shorter than the dominant dynamical timescales) nonlinear processes, and ξ represents a stochastic white noise forcing that may be spatially coherent. For a given initial state x(t), the most likely state at a later time t + τ is: where we define the propagator operator G(τ) = exp(Lτ). Determining Equation 1 from observed covariances results in a Linear Inverse Model (LIM; Penland & Sardeshmukh, 1995). LIMs are typically low-order models that are constructed from the leading Principal Component time series of a reduced empirical orthogonal function (EOF) space, but with seasonal forecast skill comparable to state-of-the-art global climate models . Further details concerning the LIM and its construction may be found in Penland and Sardeshmukh (1995), Newman et al. (2011), Capotondi andSardeshmukh (2017), and XNCD.
In general, L is stable, so perturbations do not grow exponentially without bound. However, asymmetries in L (e.g., due to geographical variation of the climatological time-mean state and/or the interactions between variables) mean that its eigenvectors are generally not orthogonal. Under certain circumstances, this allows anomalies to undergo amplification over finite time intervals when different eigenvectors with similar spatial patterns, but different temporal scales, evolve from destructive to constructive interference (e.g., Farrell, 1988;Penland & Sardeshmukh, 1995). The eigenmode combinations that yield such transient growth can be determined through an analysis of the propagator matrix for the chosen time interval. This amplification can be measured over the whole domain (based on the so-called L2 norm; Text S1 and S2 in Supporting Information S1) or for specific patterns of interest. Importantly for our problem, we can use the LIM methodology to objectively identify initial conditions that optimally grow into a MHW pattern, referred to here as the "MHW-norm" (Text S1 in Supporting Information S1). This approach is applied for different "optimization times" τ, that is, different growth periods, thus highlighting MHWs of different durations. SSH provides a measure of upper ocean heat content and system's memory. Thus, our approach differs from that of XNCD, which only included SSTAs in their LIM's state vector. Including SSH in the LIM improves its ability to capture the evolution and diversity of ENSO events (Capotondi & Sardeshmukh, 2015;Newman et al., 2011), and might be expected to affect the development and timing of MHWs. Other specific details of the LIM used for this paper are in the supplemental material (Text S1 in Supporting Information S1).
Monthly SST and SSH data were obtained from the European Centre for Medium-Range Weather Forecasting (ECMWF) Ocean Reanalysis System 4 (ORAS4; Balmaseda et al., 2013), available to us (at the time of our analysis) from January 1958 to December 2015. To examine the large-scale atmospheric fields associated with the development and evolution of the optimal precursors' patterns we used sea level pressure (SLP) from the National Center for Environmental Prediction (NCEP)/National Center for Atmospheric Research (NCAR) reanalysis (Kalnay et al., 1996).

Constructing the Northeast Pacific MHW-Norm
Following XNCD's approach, we first construct a Northeast Pacific index ("NEPac"; Figure 1a) from 3-month running mean SSTAs averaged within the region 35°-55°N, 160°-130°W (box in Figure 1b, lag 0). As in previous studies of persistent climate anomalies (e.g., Breeden et al., 2020; XNCD), we define a MHW event as occurring when the NEPac index first exceeds a given amplitude and then remains above it for a given time period. Using XNCD's "intensity" and "duration" thresholds, one standard deviation and 5 months, respectively, we find seven separate events. Note that XNCD found a few additional events for the same criteria due to the use of a longer time period and a different SST data set. The "onset" time of each event, when the index first exceeds the amplitude threshold, is marked by a light-blue vertical line in Figure 1a. The composite evolution of the SSTAs for these events, from 12 months prior to onset time (lag 0) to 12 months afterward, shows weak positive anomalies in the Gulf of Alaska and in the subtropical North Pacific at lag −12 months, while negative anomalies are seen in the eastern-central equatorial Pacific (Figure 1b). The initial anomalies in the Gulf of Alaska grow in intensity over the subsequent months, as an El Niño develops in the central equatorial Pacific (CP El Niño). A similar evolution is obtained when the 2013-2016 event is excluded from the composite (not shown).
The composite MHW evolution also shows a pronounced SSH signature (Figure 1c), which was not examined in XNCD. Anomalies in the tropical Pacific underscore the El Niño development with the initial negative equatorial anomalies (shallower thermocline) slowly being replaced by positive anomalies, consistent with off-equatorial downwelling Rossby waves (positive SSH anomalies in the off-equatorial tropics) propagating to the western Pacific, continuing along the boundary to the equator as coastal Kelvin waves, and leading, 12 months after MHW onset, to the east-west SSH dipole typical of mature ENSO phases (e.g., Capotondi, Deser, et al., 2020;Meinen & McPhaden, 2000). Note that XNCD found that more extreme threshold criteria led to a higher amplitude composite and also a reduced sample size but did not qualitatively impact their analysis.
Finally, to construct the MHW norm, we create a "MHW onset composite" by retaining only anomalies north of 25°N in the lag 0 SST composite (Figure 1d). This pattern is then projected on the reduced EOF space of the LIM, resulting in the MHW-norm pattern (r SST in Text S1 of Supporting Information S1) shown in Figure 1e. SSH anomalies (SSHAs) associated with the initial MHW composite have a very small signature north of 25°N when represented in the reduced EOF space, and therefore have not been considered in the definition of the MHW-norm.

Optimal Structures Under the MHW-Norm
The two leading optimal structures maximizing SST and SSH growth over the full domain, according to the L2 norm, correspond to canonical and Central Pacific (CP) ENSO events (Text S2, Figures S1, and S2 in Supporting Information S1). We then isolate the specific initial condition x MHW (0) maximizing Northeast Pacific MHW growth using the pattern in Figure 1e as our target (the "MHW-norm," Text S1 in Supporting Information S1).
To quantify MHW anomaly growth, rather than overall anomaly growth, we first define MHW amplitude as the projection of an anomaly onto the MHW-norm. MHW growth is then defined as the difference between the MHW amplitudes for an evolving anomaly at its optimal (initial) and maximum (evolved) states. This varies for each optimization time interval τ, as shown in Figure 2a, with a rapid increase at short optimization times, peak potential growth at around 7 months, and then a slow decrease at longer optimization times. Longer τ values allow more time for amplification from a unit initial condition, but growth cannot occur indefinitely. The optimization time corresponding to the peak growth thus represents the interval over which the greatest possible growth can occur, with longer time intervals experiencing less net growth over that interval.
The SST and SSH anomalies for the MHW τ = 9 months optimal (that is, initial conditions optimizing MHW growth over a 9 months interval) are shown in Figures 2b and 2c, respectively, along with their subsequent evolution beneath. The initial patterns and evolutions obtained at optimization times between 6 and 12 months (where growth is largest) are very similar (e.g., Figure S3 in Supporting Information S1). Initially (t = 0), weak SSTAs are seen along the Gulf of Alaska rim, in the northern subtropics and along the equator, while larger SSH signals are seen, notably in the off-equatorial northern tropics (likely associated with Rossby wave activity) and along the Kuroshio Extension. Weak positive SSH anomalies are also present in the Northeast Pacific. As time progresses, large SSTAs develop in the center of the Gulf of Alaska, concurrent with the development of the Pacific Meridional Mode (PMM; Chiang & Vimont, 2004) in the subtropics (positive SSTAs extending southwestward from the coast of California toward the equatorial Pacific) and CP El Niño conditions along the equator. SSHAs also evolve into conditions typical of CP events (Capotondi, Wittenberg, et al., 2020): larger heat content (positive SSH) in the central-eastern equatorial Pacific and reduced heat content (negative SSH) in the far western Pacific.
Regression of SLP on the combined SST/SSH optimal time series indicates initial conditions characterized by the dipole pattern typical of the North Pacific Oscillation (NPO; Linkin & Nigam, 2008;Rogers, 1981), as previously noted for individual observed MHW events (Amaya et al., 2020;Bond et al., 2015), which slowly evolves into a deepened Aleutian Low, as El Niño develops along the equator (Di Lorenzo & Mantua, 2016). Key characteristics of the MHW-Optimal are compared in Figure 3 for different optimization time intervals (τ = 6, 9, or 12 months). In each case, the MHW amplitude grows for a period of approximately τ months, by which time it reaches its maximum value (Figures 3a-3c). Exploiting the linearity of the LIM, we examine the relative contribution of the initial SST and SSH signals to the MHW growth by zeroing either the SST or SSH components of the initial MHW optimal. While the absence of the initial SSTAs only slightly reduces MHW growth (red lines in Figures 3a-3c), removal of the initial SSH signal notably hinders growth (blue lines). Inspection of the initial evolution of the SSTAs in the two cases (SSTA = 0 or SSHA = 0 at t = 0) shows that the presence of initial SSHAs allows the generation of SSTAs already at t = 1, even if they were set to zero at t = 0 ( Figure S4 in Supporting Information S1). The time evolution of SST and SSH anomalies averaged within an eastern Pacific longitude band (dotted meridional lines in Figure 1b) shows anomaly growth both in the northeast Pacific MHW region and, slightly lagged, in the tropics (middle row of Figure 3). For longer optimization times, MHW anomalies have reduced amplitude but increased duration (cf., Figures 3d and 3f), with year-long events having a double maximum similar to the 2014-2016 event. Enhanced persistence and an eastward shift of MHW anomalies is associated with longer-lasting warm conditions in the equatorial Pacific, as found by XNCD.
The potential growth captured by the MHW optimals is observed to occur. Comparison of the projection of the SST/SSH fields on the optimal patterns at time t and the observed NEPac index at time t + τ (Figure 3, bottom panels) shows fairly linear relationships, with statistically significant correlation coefficients 0.61, 0.54, Units are arbitrary, but normalization is consistent among panels. (g-i) The projection coefficients of observed anomalies on the corresponding optimal patterns, versus the observed NEPac index (g) 6, (h) 9, or (i) 12 months later. The dashed lines indicate the linear fit using a least squares approach. Correlations (insets) are significant at the 99% level based on an effective number of degrees of freedom (Bretherton et al., 1999). Black dots indicate the NEPac index at the onset times of the seven identified MHWs (y-axis) versus the projection coefficients τ months earlier (x-axis).
10.1029/2021GL097350 7 of 10 and 0.50 for optimization times of 6, 9, and 12 months, respectively. Scatter is consistent with the presence of local stochastic atmospheric forcing of the MHW events, which may impact individual event amplitude (Bond et al., 2015, XNCD).

Large-Scale Dynamics Linking MHW and CP El Niño Events
The connection between Northeast Pacific and equatorial CP warming raises the question of how the North Pacific and the tropics are linked. To further investigate this issue, we examine the eigenmodes of the LIM dynamical operator L to identify which most contribute to MHW optimals and their growth (Text S3, Figures  S5, and S6 in Supporting Information S1). The eigenmodes of L ( Figure S5 in Supporting Information S1) are all stable and can be either stationary (real) or propagating (complex conjugate pairs). Since we are interested in events lasting at least two seasons, we focus here on eigenmodes with e-folding times greater than 5-6 months, although other more strongly damped eigenmodes (shorter e-folding times) may contribute to initial anomaly growth and evolution. We identify three such eigenmodes ( Figure S6 in Supporting Information S1) that correspond to different flavors of ENSO: an "NP-CP" eigenmode with SSTAs spanning both the North Pacific and the central equatorial Pacific, and two eigenmodes whose anomalies are primarily confined within the tropics. Of these, the NP-CP eigenmode appears to play the key role in MHW growth, with its most energetic phase ( Figure  S6b in Supporting Information S1) dominating the mature phase of the MHW optimal (e.g., Figure 2b, t = 9) for both SST and SSH. Since for any nonnormal linear dynamical system, the optimal initial condition for asymptotic growth of an eigenmode is not itself, but rather its adjoint (i.e., the corresponding eigenvector of L † , see Text S3 in Supporting Information S1 for details ;Farrell, 1988), we compared evolution starting from an MHW optimal to that starting from an initial state corresponding to the adjoint of the NP-CP eigenmode (Figure 4). The pattern of the initial adjoint, with appropriate normalization and phase rotation, is quite similar to the MHW optimal (top row of Figure 4), and the resulting evolution of both is also strikingly similar. That is, the initial condition leading to optimal development of Northeast Pacific MHWs is essentially the initial condition leading to optimal development of the NP-CP eigenmode. Similar comparisons are found at different optimization times ( Figure S7 in Supporting Information S1), although the correspondence between adjoint and optimal is reduced at shorter τ, since the influence of other faster decaying eigenmodes is greater then.

Conclusions
A LIM is used to objectively identify large-scale precursors conducive to the dynamical development of Northeast Pacific MHWs on seasonal time scales. Notably, we find that Northeast Pacific MHWs and subsequent CP El Niño events may be dynamically linked through the optimal evolution of a single dynamical mode, the "NP-CP" eigenmode, whose most energetic phase includes anomalies in the Northeast Pacific, the northern subtropics and the central tropical Pacific. The spatial structure of the mode is consistent with the two-way teleconnection between the North Pacific and the tropics identified in previous studies as a source of extended climate memory and decadal variability Zhao & Di Lorenzo, 2020). The evolution of the mode is also consistent with the development of those El Niño events that impact US West Coast SSTs (Capotondi, Sardeshmukh, et al., 2019). This mode was also identified as a component of the Pacific Decadal Oscillation (Newman et al., 2016). Moreover, the pronounced similarity of the optimal precursor of Northeast Pacific MHWs to the adjoint of the NP-CP eigenmode is particularly striking. Since the "stochastic optimal" noise forcing of an eigenmode is related to its adjoint (Moore & Kleeman, 1999), this similarity suggests that some MHW and CP El Niño events may be first triggered by the same initial noise event, with subsequent deterministic development from MHW through CP El Niño due to dynamical processes that are not yet fully understood.
The LIM used in this study does not explicitly account for seasonal variation of the deterministic dynamics, which we know are important for both ENSO evolution, and for the processes underpinning MHW growth. The onset times (September-to-February) of six of the seven MHWs identified in our record (Figure 1a) are consistent with the seasonal development of the Pacific Meridional Mode, and the subsequent evolution of CP El Niño events, as discussed in a recent study that used a seasonally varying LIM (Vimont et al., 2022), suggesting that some seasonal aspects of MHW and CP ENSO evolution may be implicitly captured by our LIM. A seasonally varying LIM  will be considered in future studies.