Drivers of the extreme North Atlantic marine heatwave during 2023

Observational and reanalysis data

The atmospheric and air–sea flux fields analysed in this study are taken from the European Centre for Medium-Range Weather Forecasts (ECMWF) Reanalysis Version 5 (ERA5)³⁸, using daily mean fields derived from hourly data and monthly averaged fields from 1981 to 2023 over a 1/4° horizontal resolution grid. Variables analysed from ERA5 include surface wind speed (m s⁻¹) and surface wind vectors, surface air–sea heat fluxes (W m⁻²), mean sea-level pressure (hPa) and cloud cover fraction. SST (°C) is also taken from ERA5. The net air–sea heat flux is the sum of incoming shortwave radiation less outgoing latent, sensible and longwave radiative heat fluxes.

The surface net shortwave (solar) radiation corresponds to the amount of solar radiation that reaches the surface of the Earth (both direct and diffuse) minus the amount reflected by the Earth’s surface (which is governed by the albedo). In the ERA5 reanalysis, only climatological aerosol effects and aerosol changes due to large volcanic eruptions are incorporated into the estimate of the net surface shortwave (solar) radiation. In particular, the incoming solar radiation is partly reflected back to space because of both clouds as well as particles in the atmosphere (aerosols), but without interannual variability in aerosol effects aside from volcanic eruptions. The remainder is incident on the Earth’s surface, for which some of it is reflected owing to surface albedo effects. ERA5 assimilates cloud liquid water from a variety of satellites, as detailed in the ERA5 online documentation, formulating cloud fraction and cloud albedo as surface and single-level parameters of the reanalysis. To test the sensitivity of our results to aerosol effects from reduced shipping emissions, we further analysed the observed and modelled MLT budgets by adding 1 W m⁻² to the solar radiation terms for all months in 2023, and also to just the summer months May–August when solar radiation effects are at their maximum. The added heat flux of +1 W m⁻² is approximately double the basin-average estimate for the impact of reduced shipping emissions over the North Atlantic of 0.56 W m⁻² in ref. ³¹, which is at the high end of available estimates. In all cases our results remain robust to this addition of 1 W m⁻² solar radiation, with negligible changes to the results shown (Extended Data Fig. 11).

Subsurface ocean temperature and MLDs are obtained from an observation-based dataset; the Institute of Atmospheric Physics (IAP) hydrographic data^39,40,41 version 3 available from 1981 to 2023 with 1° horizontal resolution (hereafter referred to as ‘IAP data’). The main latitude band of interest spans the North Atlantic from the equator to 60° N, which mostly avoids regions of low data coverage, such as the high-latitude regions, particularly ice-covered areas. Sparse coverage by Argo floats is, however, a problem in marginal seas, although the gridded IAP climatology merges both Argo and other historical hydrographic observations, including XBTs. Nonetheless, reliable reconstruction of interior ocean temperature and MLD is only available at monthly mean resolution. The Gibbs SeaWater (GSW) Oceanographic Toolbox⁴² from the TEOS-10 software⁴³ is applied to convert the IAP in situ temperature to conservative temperature⁴⁴. Static stability of the water column on timescales of months is achieved by applying vertical stabilization software⁴⁵. The surface-referenced potential temperature, θ, and potential density, ρ, are then calculated using functions from the GSW Oceanographic Toolbox.

Ocean MLD is calculated as the depth at which the monthly averaged potential density in the water column first exceeds the density at 10-m depth by 0.125 kg m⁻³, following previous work⁴⁶. The MLD anomalies and the related MLT budget analyses were repeated using a density threshold of 0.03 kg m⁻³ and the results presented here are robust to this choice. MLD anomalies were also re-evaluated using a reference depth of 5 m instead of 10 m and robust patterns were obtained. The separate contributions of T and S to the surface density trends were also calculated (Extended Data Fig. 8). This was done by combining the monthly climatological salinity field (averaged during 1981–2010) with the time-varying temperature fields to derive the temperature contribution. When estimating the salinity effects on density trends, the opposite was done, namely, using the climatological mean temperature fields alongside the time-varying salinity fields.

The robustness of the estimated 2023 anomalies in MLD derived from IAPv3 data was evaluated by examining the corresponding MLD anomalies in the IAPv4, ORAS5 and ensemble mean of EN4.2.2 datasets (EN4-ESM). Overall consistent patterns of record shallow MLD anomalies were obtained during 2023 in all four datasets (Extended Data Fig. 10). However, there are variations in the long-term multidecadal trends in MLD and upper ocean temperature gradients estimated by the different products (Extended Data Fig. 10), consistent with previous findings for the North Pacific³⁶. The differences in long-term trends are in part the result of data sparsity before the Argo period, although variations can also be seen in the more recent record. Note also that the EN reanalysis is less reliable for long-term trend estimates⁴⁷, as missing values relax to climatology in the absence of any observations. Nonetheless, in all ocean reanalyses considered, 2023 sees a sudden anomalous shoaling of the mixed layer compared with previous years.

Ocean model simulations

Two ocean model simulations forced by different observation-based atmospheric reanalysis fields, namely, ERA5 and JRA55-do, were integrated to evaluate the driving mechanisms of the North Atlantic peak warming during 2023. The two model simulations both capture the temporal evolution of the onset of surface layer warming over the North Atlantic during 2023, particularly the ERA5-forced simulation (Fig. 1a). The model configurations are taken from the Australian Community Climate and Earth System Simulator Ocean Model Version 2 (ACCESS-OM2) with the configurations run at a nominal 0.25° horizontal resolution with 50 vertical levels (surface grid cell thickness of 2.3 m). ACCESS-OM2 is a global ocean–sea-ice model driven by a prescribed atmosphere; the version used here has been described previously⁴⁸, with improvements described by follow-up work⁴⁹. The ocean model component is the Modular Ocean Model (MOM) version 5.1 (ref. ⁵⁰), coupled to the Los Alamos sea ice model version 5.1.2 (ref. ⁵¹). Vertical mixing in the boundary layer is parameterized using the K profile parameterization (KPP⁵²). Further model details can be found in previously published work^48,49.

The forcing of the two ocean model simulations is as follows. In the main run analysed here, referred to as the ERA5 simulation, the model is forced by ERA5 atmospheric reanalysis fields³⁸. In the second simulation, referred to as the JRA55 run, the Japanese 55-year Reanalysis (JRA55-do) dataset for forcing ocean–sea-ice models is used. In the JRA55 run, after spinning up over six cycles forced by the 1958–2018 JRA55-do v1.4 fields (following the OMIP-5 spinup protocol⁵³), the sixth cycle is extended using JRA55-do v1.5.0 for 2019 and the delayed-mode JRA55-do v1.5.0.1 to the end of 2023 (ref. ⁵⁴). The ERA5 run is initialized at the end of 1957 and then run over 1958–2023 using ERA5 atmospheric forcing, with runoff taken from JRA55-do (as for the JRA55 run). Both model runs were initialized using World Ocean Atlas 2013 v2 temperature and salinity fields^55,56.

The model data analysed here include SST and surface heat fluxes from both the ERA5 and JRA55 simulations and, for the ERA5 run, all variables required to reconstruct a surface MLT budget (detailed below) both during 2023 and also during all years of the model simulation from 1981 onwards, to construct the baseline 1981–2010 climatological mean MLT budget from which we evaluate the 2023 anomaly fields. The North Atlantic SST evolution of the two model runs is, overall, consistent with observations (Fig. 1a), as is the pattern of MLD shoaling across the North Atlantic during 2023 (Extended Data Fig. 6). The magnitude of air–sea surface heat flux anomalies during 2023 is also consistent with observations, as revealed in the breakdown of the anomalous net surface heat flux contributions to the MLT budget during May–August (compare Figs. 4a and 5a). The geographic pattern of net surface-heat-flux-driven MLT warming also compares well between both model runs and observations.

Surface layer temperature budget

A budget for the upper ocean MLT can be written as

in which θ_h is the MLT, t is time, Q_net is the net surface air–sea heat flux (W m⁻²), Q_h is the heat flux corresponding to shortwave radiation penetration across the base of the mixed layer (also in W m⁻²), ρ₀ is seawater density (taken to be 1,027 kg m⁻³), c_p = 3,992 J K⁻¹ kg⁻¹ is the heat capacity of seawater, h is the MLD (m), advection represents net heating owing to three-dimensional ocean advection, entrainment represents the heat flux associated with mixed layer shoaling or deepening and the mixing term includes the net heat flux owing to parameterized subgrid-scale processes such as vertical mixing and eddy-induced mixing. The entrainment term can be written as:

representing entrainment associated with mixed layer shoaling or deepening, with θ_ent the temperature of detrained/entrained water and θ_h the MLT.

The shortwave penetration across the base of the mixed layer (Q_h) in the model is computed online and output as a heat budget diagnostic by MOM5, in which shortwave radiation is distributed into the ocean interior as a function of depth⁵⁷, modulated by the seasonal climatological chlorophyll distribution⁵⁸. The penetration of shortwave radiation across the base of the mixed layer for the observationally based estimate is taken as Q_h = Q_SW. F(z), in which Q_SW is the shortwave radiation across the air–sea interface from the ERA5 reanalysis and F(z) is the exponential decay function⁵⁹,

where z is depth (positive downwards), R = 0.58, h₁ = 0.35 m and h₂ = 23 m.

The two ocean model simulations give, overall, robust results in terms of the magnitude of the anomalous surface heat flux contributions to the MLT tendency terms (figure not shown), including the breakdown into shortwave, longwave, sensible and latent heat flux terms. The focus of our more detailed MLT budget analysis is the ERA5 run, which was integrated with all required terms saved to quantify monthly variations in the MLT balance, including three-dimensional ocean advection and mixing, as described below. The temperature budget of the ERA5 model can also be compared with estimates using the observationally based ERA5 reanalysis combined with IAP estimates of surface ocean temperature and mixed layer fields. However, insufficient observational data are available to constrain the other surface MLT budget terms, such as ocean advection and mixing, as well as vertical entrainment. We can nonetheless compare the surface air–sea heat flux contributions across the model and observed reanalysis fields.

All terms in the model MLT budget, including all components of ocean advection, mixing and entrainment, can be diagnosed by integrating the fully closed cell-by-cell model heat budget^50,60,61 over the surface mixed layer (defined here using the same definition as the observed MLD, described above) and dividing by the surface MLD. Here this integration is performed using monthly averaged heat budget diagnostics and the monthly averaged MLD. The entrainment term, the only term not explicitly included in the three-dimensional model heat budget, is computed by residual as the difference between the tendency of the MLT (computed by taking the difference of snapshots of the temperature and the MLD at the beginning and end of each month) and the total mixed layer heat content tendency (that is, the sum of all point-by-point heat budget processes, divided by ρ₀c_ph). Such a calculation produces a closed budget but neglects correlations between submonthly variations in the MLD and cell-by-cell model heat budget terms. However, a comparison performed in 2023 between these monthly averaged budget terms and similar terms obtained using daily averages (only available for 2023) showed only small differences (order less than 5%; figure not shown).

The MLT budget analysis for observations (Fig. 5) is derived by comparing the estimated mixed layer heat storage rate to the heating and cooling driven by surface heat fluxes alone (less shortwave penetration at the base of the mixed layer), with the residual implicitly including all ocean advection, entrainment and mixing terms, as well as any errors owing to incomplete data coverage in space and time, uncertainty in estimating MLDs and temperatures and uncertainty in the ERA5 net air–sea heat fluxes (both in 2023 and in the baseline period 1981–2010). As for the model simulations, the observational budget also neglects submonthly correlations between the surface heat flux and MLD variations, although analysis of the model budget at daily and monthly timescales suggests that this term is small. The MLT tendency, ∂θ_h/∂t, is calculated from monthly observations as a centred second-order difference of monthly mean MLT. With only monthly MLD values available from observations, this results in substantial smoothing of the MLT warming signal compared with the daily resolved model MLT budget.

The pattern of residual terms seen in the observed calculation reveals features that include unresolved temperature changes due to ocean circulation effects both laterally, such as in the Gulf Stream, and also vertically, because of entrainment and vertical mixing at the base of the mixed layer. By contrast, the MLT budget of the model explicitly resolves the contributions from ocean circulation and mixing, and so all terms can be shown.

The shortwave radiation anomalies in the MLT budget are further diagnosed by recalculating these anomalies during 2023 but separately considering the effects on MLT of 2023 anomalies in MLD, surface incoming shortwave radiation and radiation through the base of the mixed layer. This calculation simply recomputes the 2023 MLT budget holding all other variables at their climatological mean, then separately including the 2023 anomalies in MLD (MLD′), surface shortwave heat flux (Q_SW′) and anomalies owing to the radiative flux across the base of the mixed layer (Q_SW,H′). The resulting values are shown from the ERA5-forced model run and observations in Figs. 4a and 5a and Extended Data Fig. 9, respectively. This decomposition indicates the sign and magnitude of the shortwave radiation component of the MLT budget anomalies assuming that only one of these three factors (MLD′, Q_SW′ and Q_SW,H′) varied with its 2023 values. For example, the MLD′ values indicate the size of the 2023 shortwave radiation term assuming that the incoming shortwave radiation, and the radiation through the base of the mixed layer, followed the climatological mean. The Q_SW′ term conversely takes MLD and Q_SW,H to follow the climatological mean, with Q_SW varying with its 2023 values. Last, Q_SW,H′ assumes that Q_SW and MLD follow the climatological mean, with only the radiative flux through the base of the mixed layer evolving with its 2023 values. This approach is a simple way to evaluate what factors were the most important in generating the summertime warming in the North Atlantic during 2023, although the contributions are not additive to unity by definition.