We use temperature, currents, and MLD from the GLORYS2-V4 Global Ocean Ensemble Physics Reanalysis (hereafter GLORYS; https://data.marine.copernicus.eu/product/GLOBAL_REANALYSIS_PHY_001_031/description) (Garric and Parent, 2017) for the period January 1993 to August 2019. This reanalysis is based on the NEMO ocean general circulation model version 3.1 (https://www.nemo-ocean.eu/). The ocean model assimilates satellite SST from AVHRR and AMSR-E (1/4° resolution), along-track sea-level anomalies derived from satellite altimetry, and temperature and salinity profiles from ARGO floats since 2002. The ocean model has a horizontal resolution of 1/4° curvilinear orthogonal grid and 75 vertical levels. NEMO’s surface boundary conditions are taken from ERA-Interim global atmospheric reanalysis (https://www.ecmwf.int/en/forecasts/datasets/reanalysis-datasets/era-interim) and includes zonal and meridional surface wind, latent and sensible heat fluxes, and net shortwave and longwave radiation. The ERA-Interim database, also used in the present study, has a temporal and spatial resolution of 6 h and approximately 80 km (3/4°), respectively, and 60 vertical levels from the surface to 0.1 hPa (Berrisford et al., 2011; Dee et al., 2011). Only sea surface data over the period from 1993 to August 2019 were used in this study. To achieve the same horizontal resolution as GLORYS (i.e., 1/4°), the ERA variables were interpolated using an Akima spline.

MHWs were defined using the mixed layer mean temperature (MLT) following the methodology proposed by Hobday et al. (2016). When the daily MLT in a specific grid point exceeds a climatological threshold value for a period of at least 5 days, it is considered to be impacted by a MHW. Daily threshold values correspond to the 90th percentile of the MLT for each grid point. For the calculation of the 90th percentile for a specific day of the year, we used an 11-day moving window centered on the specific day and the 26-year-long time series (see Figure S1 in Supporting Information). Finally, the threshold values were smoothed using a 31-day moving mean (Hobday et al., 2016).

To characterize the MHW events, different metrics were used: duration, spatial extent, and mean and maximum intensities (Hobday et al., 2016). We also quantified the frequency of MHWs inside the study region based on the number of events per year.

The equation for the tendency of the MLT is used to evaluate the relative importance of the main MHW drivers. The different terms involved in this equation are associated with different ocean and/or atmospheric processes and conditions, and they can contribute to increasing or decreasing the MLT. Following Oliver et al. (2021), we write the heat budget in the mixed layer as

where T ¯ is the vertical mean temperature in the mixed layer, t is time, u=(u,v) is the horizontal velocity (obtained from GLORYS), u_{- h} is u at z=–h (the base of the MLD), and w is the vertical component of the velocity estimated from the wind stress curl (Ekman pumping, EP). Close to the coast, we add a term directly related to the cross-shore Ekman transport (M), assuming that the cross-shore Ekman transport is locally compensated by a vertical transport inside a distance given by the first baroclinic Rossby radius (Chelton et al., 1998). EP and M were calculated from the ERA-interim surface wind field. The vectorial operator ∇ is in the horizontal plane, Q are the ASHFs (related to shortwave and longwave radiations, and sensible and latent heats), ρ is the density of seawater, and C_p is the specific heat of seawater at constant pressure. The top bar indicates vertically averaged quantities in the mixing layer, a subindex h indicates that the quantity is evaluated at a depth z=–h(x,y,t), i.e., at the base of the mixed layer.

The first term on the right-hand side of (1) is the horizontal heat (temperature) advection vertically averaged in the MLD. The second term comprises the total change of the MLD (Dh/Dt) and the entrainment of temperature into the mixed layer (herein, entrainment). The third term comprises the sum of the four terms that conform the net ASHF: the two radiative terms –the short wave (Q_SW ) and longwave (Q_LW ) fluxes– and the two turbulent heat flux terms –sensible (Q_SENS ) and latent (Q_LAT ) heat. Herein, we call this third term air-sea heat fluxes (ASHFs). The last term is a residual associated with turbulent horizontal and vertical mixing, which is assumed to play a minor role compared to the other terms (e.g., Holbrook et al., 2019; Marin et al., 2022). We identify the first three terms of the right-hand side of (1) as the possible drivers of the MHW events. Time series of HA, ASHF, and entrainment are calculated by integrating only periods and zones affected by MHWs. The resulting time series for each term represents the partial contribution to the total temperature change over the region impacted by MHWs per day (in °Ckm² day^-1).

The relative importance of the drivers associated with each MHW is calculated by comparing the daily magnitudes of a specific driver with the daily time series of the MLT tendency (i.e., ∂ T ¯ / ∂ t from GLORYS), both integrated over the MHW-impacted area. Specifically, the comparison considers the 75th percentile of the driver and the 50th percentile of the MLT tendency, a procedure explained as follows: percentiles are calculated based on the respective daily time series during each MHW, e.g., percentiles for a MHW lasting 40 days are calculated from the respective distributions of these variables within those 40 days. The choice of these percentiles for both the MLT tendency and the drivers are discussed in the supplementary material ( Figures S4 and S5 in Supporting Information). If the 75th percentile of a certain driver is higher than the 50th percentile of the MLT tendency for a given MHW, we consider that the driver plays a significant role in the generation of this MHW event. When more than one driver plays a significant role in the generation of the same MHW, we compare the magnitude of these drivers to evaluate the dominant driver. Thus, the driver whose 50th percentile exceeds the 75th percentile of the other relevant driver is considered dominant. This means that the warming generated by the dominant driver exceeds the warming generated by the other driver. When neither of the drivers is dominant and all play a significant role, we refer to the MHW as a combined type.

The MHW frequency ranged from ~1 to 3 events per year ( Figure 1A ). The duration of MHW ( Figure 1B ; calculated using definition1) fluctuate in a range between 16 (25th percentile) and 20 days (75th percentile). The duration and frequency of the MHWs was inversely correlated, i.e., longer (shorter)-lasting events are less (more) frequent. Among the different metrics, intensity has been recognized as the most relevant for assessing MHW severity (e.g., Hobday et al., 2018; Oliver et al., 2021) and is frequently used for the overall MHW classification. In our study region, maximum intensity typically ranged from about 1.5 to 1.7°C ( Figure 1D ; the lower and upper limits correspond to the 25th and 75th percentile, respectively), whereas mean intensity was 1.2 ± 0.2°C ( Figure 1C ). The intensity increased consistently toward the coast ( Figures 1C, D ), particularly when plotting maximum annual intensities ( Figure 1D ). Unlike intensity, frequency and duration do not have a clear spatial pattern.

Between January 1993 and December 2018, two extreme, eight severe, 31 strong, and 18 moderate MHWs occurred off Chile ( Figure 2 ). These MHWs showed a wide range of duration, spanning from 16 to 20 days (9 to 44 days) according to definition 1 (2) (see section 2.3). The spatial extent of the MHWs (i.e., the sum of all pixels impacted by each MHW) ranged from 3×10³–4×10⁴ km² ( Figure 2 ). In general, it is expected that the higher the category of a certain MHW, the greater its intensity, duration, and spatial extent. The category of MHW correlates well with the duration (r=0.73), maximum intensity reached during that category (r=0.77), and spatial extent (r=0.85) of the events. Although this is true for many MHWs, some events show considerably greater duration or spatial extent than others of the same category ( Figure 2 ). For instance, the strong 249-day-long MHW labeled 2015 is much longer than all the strong and severe MHWs, and the area of the severe MHW labeled 2016 covers 18700km², an area greater than that of all the severe MHWs.

Based on a heat budget analysis in the surface mixed layer, we examined the main processes involved in the formation and evolution of the MHWs. To estimate the relative importance of the different drivers represented in (1) and their seasonal contributions, we calculate their respective seasonal cycles as daily averages for the 26-year period ( Figure 3 , black curves). Note that the seasonal cycles for all drivers were calculated using only periods and regions impacted by MHWs. During most of austral spring and summer (September-February), radiative forcing is commonly the dominant term ( Figure 3A ), whereas HA is the dominant term in austral fall and winter (March-August) ( Figure 3B ). Entrainment is usually weakly negative during MHW events, reducing the mixed-layer temperature. In general, the contribution of entrainment is very small (<8% of the magnitude of ASHFs and HA separately, Figure 3C ) and practically negligible in the heat balance. Solar radiation is reduced in fall and winter, increasing the relative importance of HA ( Figure 3B ). As the magnitude of HA depends on both the magnitude of the current normal to the temperature gradient and on the magnitude of this gradient and given the large magnitude of the offshore temperature gradient, relatively small positive (onshore) current anomalies are important for increasing the HA term during some MHWs. Nevertheless, large southward current anomalies along the coast may also increase HA. These particular cases are illustrated below for the long MHW related to the El Niño 1997-1998 off central Chile. We show that the magnitude of MHWs was greater from late austral spring to the end of summer ( Figure 3D ) and lowest from September to November. This magnitude is associated also with the seasonality of the frequency of the MHW observed in the study region, which are larger in austral summer and fall and minimum in spring. Note that the secondary maximum of MHW magnitude observed in April can be related to the increasing of HA.

Figures 3A–C also show the contribution of the different terms of (1), distinguishing their positive (red curves) and negative (blue curves) contributions to the heat balance. The positive (negative) contribution of the driver occurs when only values favorable to the warming (cooling) of the driver are considered. This allows better analysis of the nature of the MHW in the region. For instance, a MHW event could be dominated by two different mechanisms in different subareas within the study region, and the same driver can be adding and subtracting heat from the surface layer within different subareas of the same MHW, which would result in a rather small total effect of the driver after averaging over the whole region affected by the MHW. In such cases, the double nature of the MHW would be missed. In summary, the climatological heat balance within the MHWs showed that ASHFs practically always contribute to increasing surface layer temperature (i.e., the red and black lines are similar in Figure 3A ), whereas HA shows both positive and negative contributions –represented in Figure 3C by red and blue lines, respectively– with a predominance of positive values (black line). Entrainment plays only a secondary role and commonly contributes to reducing MLT. Then, the main drivers of the MHWs are HA and ASHFs.

Daily anomalies of ASHFs, HA, and entrainment vary widely within the MHWs, with fluctuations that are about one order of magnitude larger than the amplitude of their seasonal cycle ( Figure 4 and Figure S3 in Supporting Information). The red line shows the standard deviation, which was calculated using a moving window of 365 days (using this moving window, we avoid overweighting the influence of ENSO). The values over each red line correspond to periods when the drivers are more important; these periods are used to illustrate the specific processes associated with the different drivers. Note that most of the peaks in the ASHF and HA anomalies ( Figures 4A, B ) are related to El Niño. This phenomenon is most evident in the extreme events of 1997-1998 and 2015-2016, but the impact of the “Coastal El Niño” is also clearly visible during the austral summer of 2017.In general, the spatial pattern of the HA anomalies ( Figure 5A ) resembles the observed MHW patterns of both frequency and duration ( Figures 1A, B ). These patterns result from the combination of ASHF and HA. Despite the ASHF is the main driver of most MHWs, the HA has more spatial structure in the oceanic region than the relatively homogenous ASHF pattern (see Figure 7A below). Note that there are periods not related to ENSO when HA and ASHF anomalies are larger than their respective threshold red curves.

During El Niño, the magnitude of MHWs (mean value ~6.9×10⁴°C km²) is greater than during La Niña (mean value ~5.7x10^4° km²) and neutral periods (mean value ~4.3x10⁴ km²) ( Figure 6A ). Interestingly, the magnitude of MHWs during La Niña is greater than during neutral conditions ( Figures 6C, D ). During El Niño, the most intense MLT anomalies are observed near the coast ( Figure 6B ), whereas during La Niña, larger magnitudes are present in the offshore region ( Figure 6C ). In summary, ENSO cycles increase the probability of more intense MHWs in our study region.

The warming pattern (in °Cd^-1) related to the ASHFs shows positive values across the region. Larger values are found near the coast ( Figure 7A ), where the mixed layer is shallower ( Figure 7B ).Among the different terms involved in the air-sea heat exchange ( Figure 8 ), the term related to latent heat flux ( Figure 8C ) is the most relevant in creating the spatial pattern observed in Figure 7A . In contrast to the pattern observed from the latent heat flux anomalies ( Figure 8C ), patterns of short and longwave radiation anomalies show mainly meridional variability and smaller magnitude ( Figures 8A, B ). The sensible heat flux anomalies are spatially more homogeneous and slightly positive ( Figure 8D ).

The spatial pattern of the HA ( Figure 5A , in °C/d) shows that all the region tends to be warmed up by this driver, but the largest values are observed along a narrow coastal band. This warming is more intense south of 32°S, affecting a larger offshore region. The warming pattern is related to both eastward (onshore) current anomalies and a large zonal temperature gradient in all the region, being more intense south of 32°S ( Figure 5B, C ). It is worth noting that the HA is, in general, smaller than the ASHF in our study region and entertainment anomalies are commonly much smaller (<8%). As stated before, a MHW could be generated and maintained by different drivers acting at different times and places. Most long MHWs are driven by a combination of ASHF and HA.

To better illustrate our results, we analyze two extreme MHWs: one during 1997-1998, related to the large El Niño event, and the other during 2016-2017, related to the “coastal El Niño”. The 1997-1998 El Niño was the strongest on record (e.g., McPhaden, 1999) and generated the most extreme and longest MHW observed in our study region ( Figure 2 ). During this MHW, both HA and ASHFs were relevant. Previous to the 1997-1998 El Niño, during January-February, a MHW impacted the region ( Figure 9A ). Its growth was initially driven by an increase in ASHFs ( Figure 9C ), whereas the contribution of HA remained relatively low. After that, the magnitude of the MHW dropped nearly to 0 but remained positive. During April-June 1997, the magnitude of the MHW began to grow again ( Figure 9A ) related to the onset of the strong El Niño in the region ( Ulloa et al., 2001 ). At that time, the HA increased considerably ( Figure 9B ), whereas the magnitude of ASHFs was close to zero. In austral winter, the ASHFs started increasing and peaked in November 1997 ( Figure 9C ), contributing to maintaining the long MHW observed in our region. Then, in December 1998, near the peak of the ENSO event, both ASHF and HA increased ( Figures 9B, C ). Thus, the extreme MHW of 1997-1998 observed off central Chile was driven by a combination of ASHF and HA. The pattern of the ASHFs term shows that during January-February 1997, ASHFs contributed to the MHW in most of the offshore region ( Figure 9E ). Then, during the onset of the 1997-1998 El Niño (April-June 1997), HA was the dominant driver with intensity showing a variable spatial pattern in the offshore region ( Figure 9F ). In contrast, during December 1997, the HA dominated very close to the coast increasing the MHW intensity ( Figure 9G ). This last pattern is consistent with large coastal trapped waves –forced by equatorial Kelvin waves– arriving in central Chile (e.g., McPhaden and Yu, 1999). In addition, the ASHF also made an important contribution in December 1997, similar to the previous period (August-November 1997), with a similar warming pattern ( Figure 9H ). Prior to the 1997-1998 El Niño, in January-February 1997, a large MHW impacted the region ( Figure 9A ). This MHW was mainly driven by ASHF, with a minor contribution from HA ( Figures 9B, C ). During this period, satellite SST data showed positive anomalies over a wide region of the South Pacific centered around 30°-35°S, reaching the coast of central Chile (c.f., Shaffer et al., 1999).

During the coastal El Niño of 2016-2017, a large MHW occurred off central Chile ( Figure 10 ). In November and December 2016, ASHF dominated the surface heat budget, driving this MHW ( Figures 10A, C first yellow-shaded period). The spatial pattern of ASHF-related warming in this period is intense and relatively homogeneous, affecting the whole zone ( Figure 10D ). Later, both ASHF and HA become relevant in supporting this MHW ( Figures 10A–C ; gray-shaded period). This occurs simultaneously with the maximum expressions of coastal El Niño. During the last stage of this MHW, the HA was stronger and dominated the heat balance ( Figures 10A, B April-May 2017; second yellow-shaded period). Finally, HA and MHW decay together during mid-2017. The heating pattern generated by HA during this last period shows large mesoscale variability with strong warming nuclei, mostly near the coast ( Figure 10G ).

Based on the main drivers, we can distinguish three types of MHW: 1) mainly forced by ASHF, 2) mainly forced by HA, and 3) combined-type, in which the relevance of both ASHF and HA is similar. Figure 11 presents a summary of the classification of MHWs. ASHF-type MHWs were the most frequent –more than 65 out of 147 events (44%) were observed in our 26-year-long study period ( Figure 11A ). These MHWs had a mean intensity of ~1.2°C and a maximum intensity of ~2.5°C ( Figure 11C ). For those MHWs that are dominated by a particular driver, HA exhibits a higher mean partial contribution to the total magnitude than ASHF (in average, 814 vs 605°C km² day^-1, respectively). During ASHF-type MHWs, the mixed layer is considerably reduced, reaching a typical depth of only ~18 m ( Figure 11D ). Thus, ASHFs can warm up more quickly in a region with a reduced mixed layer, making this driver more efficient to increase the MLT. This type of MHWs occurs when the wind is reduced, thus decreasing the mixing in the surface layer and, in turn, the MLD. ASHF-type MHW are characterized by a relatively short duration, ~24 days ( Figure 11B ), whereas their average spatial extent covers ~2.8×10⁴ km² ( Figure 11D ). The largest MHW of this type impacted ~34% of the whole study region on average.

The HA-type MHWs are considerably less frequent ( Figure 11A ): only 33 out of 147 (23%) events of this MHW occurred during the study period. As stated above, the partial contribution of HA to the total magnitude is on average ~814°C km² day^-1 during HA-type MHWs, which is larger than the partial contribution of air-sea heat fluxes during the ASHF-type MHWs. However, in this case, the mixed layer is thicker (~47 m depth) ( Figure 10D ). These MHWs have smaller maximum (~2.0°C) and slightly smaller mean intensities (~1.1°C on average) compared with the other MHW-types ( Figure 11C ). The mean duration of AH-type MHWs is ~57 days, about twice the mean duration of the ASHF-type MHWs. Moreover, their average spatial extent involves an area of ~2.4×10⁴ km² ( Figure 11D ), including a particular MHW that reached an average extension of 27% of the entire region analyzed.

During the study period, 36 combined-type MHWs were observed, representing 25% of the total 147 events. Hence, this group is more (less) frequent than HA-type (ASHF-type) MHWs ( Figure 11A ). The combined-type MHWs showed a mean intensity of ~1.2°C and maximum intensities (~2.6°C) similar to the ASHF-type MHW ( Figure 11C ). The mean partial contribution of the combined heat fluxes (ASHF + HA) to the total magnitude was 650°C km² day^-1 ( Figure 11A ). These MHWs had the longest duration (~65 days) ( Figure 11B ). As shown in the study cases above, the different drivers may act simultaneously or during different periods, thus extending the MHW duration. Consequently, combined-type MHWs had the largest spatial extents (~4.5×10⁴ km² ( Figure 11D ), being considerably superior to the other two types of MHWs. The largest MHW reached an average spatial extent of about 39% of the entire region and corresponded to the event observed during the 1997-1998 El Niño, analyzed in the previous section. The MLD during this type of MHW was on average ~36 m ( Figure 11D ), an intermediate value between those for the other two types. Finally, MHWs that do not exhibit a main contribution by the drivers analyzed presented the lowest frequency: only 12 events out of 147, i.e., 8%. They were generally weaker than the MHWs analyzed above ( Figures 9 and 10 ), showed the shortest duration (~8 days on average) ( Figure 11B ), and the smallest spatial extensions (~0.4×10⁴ km²) ( Figure 11D ).