Increased stratification intensifies surface marine heatwaves north-east of Aotearoa New Zealand in New Zealand’s Earth System model

2.1. Model description

The New Zealand Earth System Model (NZESM; Williams et al. 2016) is a fully coupled Earth System model based on its parent model, UK Earth System model (UKESM) (Kuhlbrodt et al. 2018; Storkey et al. 2018; Sellar et al. 2019). Ocean physics is simulated through Nucleus for European Modelling of the Ocean (NEMO) (Madec et al. 2017), ocean biogeochemistry through Model of Ecosystem Dynamics, nutrient Utilisation, Sequestration and Acidification (MEDUSA) (Yool et al. 2013), sea-ice through The Los Alamos sea ice model CICE (Rae et al. 2015; Hunke et al. 2017) and atmospheric physics and chemistry, including aerosols, through Met Office Unified Model (UM) (Walters et al. 2019). The main difference between NZESM and UKESM is the dynamical oceanic downscaling that has been embedded in NZESM around New Zealand to improve the representation of the oceanic circulation (Behrens et al. 2020). The dynamical downscaling has been applied from 132.7°E to 143.7°W and 60.17 to 10.75°S, with a nominal resolution of 1/5° using Adaptive Grid Refinement in Fortran (AGRIF; Debreu et al. 2008). Outside the nested region, a global tripolar 1 × 1 model grid has been used to simulate ocean properties. Further NZESM model details, validation and impacts of the dynamical downscaling in NZESM can be found in Behrens et al. (2020).

2.2. Model simulations

For this study, we use four NZESM simulations, consisting of a historical simulation from 1950 to 2014 and three future projections with different greenhouse gas emissions shared socioeconomic scenarios (SSP1–2.6, SSP2–4.5, SSP3–7.0; Meinshausen et al. 2020) following the widely used Coupled Modelling Intercomparison Project 6 protocol. To reduce computational demand and allow a direct comparison with UKESM, initial conditions from related historical simulations of UKESM in 1950 were used to initialise the historical simulations of NZESM. The future projections cover the period 2015–2100, which represent low (SSP1–2.6), moderate (SSP2–4.5) and high (SSP3–7.0) future greenhouse gas emission scenarios. These simulations have been analysed for marine heatwaves around New Zealand (Behrens et al. 2022) and used to predict deep sea coral distributions (Anderson et al. 2020, 2022; Stephenson et al. 2023).

2.3. Diagnostics

In this study, we conduct a heat budget analysis, similarly to Behrens et al. (2019), over our target region to investigate why MHW intensities continue to increase in the future, despite the projected decrease in volume and associated heat transport into the region. We compute the oceanic heat transport over the study area (Fig. 1a and Table 1) with corner points (A–E), hereafter referred to as the target region. The enclosed region is dominated by two oceanic currents, the EAUC, which is partly sourced from the TF through the section (A–B), and the ECC, which carries water southward through section (D–E) (Fig. 1a). Flow across the other boundaries does not include any major oceanic currents but net heat transport is not negligible, as specified in the results section. There is net heat transport into the region across the western section and out of the region across the northern, eastern and southern section. Modelled meridional and zonal velocities in combination with potential temperature were used to quantify the oceanic heat transport in and out of our target region, using Eqn 1:

where H_T is heat advection (kg·m² s⁻³), ρ is the density of ocean water (kg m⁻³) (set to be constant at 1025 kg m⁻³), c_p is the specific heat (kg·m² s⁻²·kg⁻¹·K⁻¹), T is the temperature (K), v is the velocity (m s⁻¹), dl is the distance north or south depending on the section (m). The sum of the oceanic heat fluxes over the boundaries was then compared with oceanic heat content (OHC, kg·m² s⁻²) anomalies of the top 1000 m over the target region by integrating potential temperature over the enclosed volume:

The anomalies are the difference between the mean over 2050 and 2099 of each future pathway and the historical mean (1995–2014). Furthermore, the ocean–atmosphere surface heat flux over the region and heat transported to below 1000 m were diagnosed to close the budget and identify future drivers for MHWs in this region. This allowed us to identify the heat sources and sinks of the region. Volume transports across the sections were also computed to quantify changes to the oceanic circulation.

Fig. 1.

(a) Schematic of surface ocean currents in our study area (cyan box): East Australian Current (EAC), East Auckland Current (EAUC) and East Cape Current (ECC). Grey contour lines show ocean bathymetry (1000 and 3000 m). Table 1 provides the locations of the boundary for our study area indicated by the corner points A–E. Median marine heatwave intensity anomalies relative to 1995–2014 as projected from New Zealand Earth System Model by 2080–2099 for (b) low (SSP1–26), (c) moderate (SSP2–4.5), and (d) high (SSP3–7.0) greenhouse gas emission scenarios. The blue hatched regions indicate areas where the change in median MHW intensity is not statistically significant (P > 0.01) compared to the historical reference period. (c) and (d) adapted from Behrens et al. (2022).

Area averages over the target region were calculated to diagnose changes to the temperature, salinity profile and water column stratification using mixed layer depths and Brunt–Väisälä frequency. Mixed layer depths were diagnosed as the depth where the potential density is 0.01 kg m^–3 greater than the density at the surface. Lastly, model wind stress was assessed to quantify changes to wind-induced mixing over our target region.

The net heat and volume transport of the EAUC and ECC were determined by the maximum transport from integration of the heat and volume transport from the coast outwards for each time step. This approach accounts for the geographic shift of the current over time. As both currents are generally adjacent to the continental shelf, we assume the absence of a strong recirculation between the coast and the current.

Future changes to the above diagnostics were quantified over the period between 2050 and 2099 in comparison with the last 20 years of the historical simulation (1995–2014). Future MHW intensities (Fig. 1b–d and Fig. 2) were diagnosed using 5-day averaged SST data following the approach outlined in Behrens et al. (2022) by using a fixed baseline from 1982 to 2012. To demonstrate that MHW intensities are projected to become more severe with climate change and that this intensification is not solely driven by the long-term warming trend, a comparison was conducted where local linear trends were removed from the 5-daily averaged SST data for the individual analysed periods while keeping the fixed baseline period unchanged. Probability distributions were constructed for our target region and normalised so that the area under the curve equated to 1 for each probability distribution (Fig. 2). The probability distribution for MHW intensities with a fixed baseline and warming trend included (Fig. 2a) substantially changes shape in the future scenarios and implies a higher likelihood for more severe MHWs compared with the historical period. Even when the local linear temperature trends are removed, the probability of more severe MHWs is projected to increase in the future as indicated by the shift towards higher MHW intensities compared with the historical period (Fig. 2b). This MHW intensification provides the motivation for this study to further investigate the drivers for MHWs and how they change north-east of New Zealand.

Fig. 2.

Probability distribution for MHW intensities over our target regions from different climate change scenarios with NZESM: (a) fixed baseline 1982–2012; and (b) fixed baseline 1982–2012 and with local linear temperature trends removed.

3.1. Heat sources north-east of NZ

The heat from subtropical waters that is transported to the EAUC through the TF is the main source of heat into the target region, shown by the large positive heat transport between 34 and 33°S over the western section (Fig. 3a and Table 1). Investigating changes to the inflow of heat helps us understand the projected changes in oceanic heat and MHWs in the region. Along the southern section, the southward heat transport associated with the ECC west of 180°E is apparent. The northern and eastern sections in Fig. 3 have been scaled by a factor of three so their weaker variability is visible. They show some interannual variability, but the persistent heat transport is not associated with any permanent oceanic currents. The volume and associated heat transport of the TF, which feeds into the EAUC, is projected to decrease owing to the spin-up and extension of the South Pacific gyre (Sen Gupta et al. 2021). Over the historical simulation of NZESM between 1950 and 2014, the net heat transport of the EAUC shows a long-term decline by ~8% per decade and the net heat transport of the ECC a decline by ~15% per decade. The trends of the heat and volume transport for the EAUC and the ECC are shown in Fig. 4a, b respectively. Overall, the future projections suggest an overall weaker volume and heat transport relative to the historical period, but that the historical strong trends ease. However, only the high emission scenario (SSP3–7.0) shows a further declining trend in both quantities for the EAUC, whereas the other scenarios suggest a positive trend from 2015 onward. For the EEC, all scenarios suggest a positive trend in volume and heat transport beyond 2015.

Fig. 3.

(a–d) Historical simulations of the top 1000 m oceanic heat transport over the four boundaries of the target region; (e–h) interannual oceanic heat transport anomalies for SSP3–7.0 relative to 1995–2014 of the historical simulation. (i–l) Mean oceanic heat transport anomalies for SSP1–2.6, SSP2–4.5, SSP3–7.0 for 2050–2099 relative to the historical simulation (1995–2014). Positive (negative) values in the heat transport in the historical simulation represent heat transported into (out of) the region. Positive (negative) anomalies represent an increase (decrease) in heat transport. Note the scaling factors for the northern and eastern boundary.

Fig. 4.

The annual mean heat transport and volume transport in Sverdrups (1 Sv = 10⁶ m³ s^–1) of (a), (c) the EAUC, and (b), (d) the ECC; the solid line represents the trend of heat transport in the second half of the century for the future scenarios and the dashed line over the period 2015–2099.

The mean historic (1995–2014) net heat transport (top 1000 m) across the four sections, West, South, East and North, is 556.4, −241.4, −137.6, −163.9 × 10¹² W respectively (see Table 1), which indicates an oceanic heat convergence of 13.6 × 10¹² W. This heat convergence has flow-on effects on the heat content, heat fluxes to the atmosphere and fluxes to the deep ocean. The heat loss into the deep ocean, below 1000 m, calculated with NZESM, is six orders of magnitude smaller than the oceanic heat transport across the boundaries and atmospheric heat losses. Hence, the heat flux to the deep ocean is considered negligible for the heat budget shown in Fig. 5b. In the SSP1–2.6, SSP2–4.5 and SSP3–7.0 scenarios, 72.9, 45.2 and 42.1% respectively of the heat transport convergence leads through the warming ocean to increased heat fluxes to the atmosphere, while the remaining 27.1, 54.8 and 57.8% in SSP1–2.6, SSP2–4.5 and SSP3–7.0 maintains the ongoing warming of this region. Less than 0.003% of the heat transport convergence gets transported to the deep ocean below 1000 m. The cumulative heat fluxes of SSP2–4.5 and SSP3–7.0 do not suggest a dependency on the climate changes scenario. Therefore, the heat fluxes, atmospheric and top 1000 m ocean heat transports in and out of the region alone cannot explain the different rates at which MHWs are projected to intensify.

Fig. 5.

(a) Top 1000-m cumulative heat fluxes (Φ) for future projections based on oceanic heat advection (solid), atmospheric fluxes (dashed) and advection into the deep ocean below 1000 m (dot–dash) for SSP1–2.6 (green), SSP2–4.5 (blue) and SSP3–7.0 (red). (b) Net heat content changes derived from the difference between oceanic and atmospheric heat fluxes from (a). Positive heat flux anomalies indicate oceanic heat transport convergence and an increase of ocean heat content over the target region.

3.2. Heat transport changes

In this section, we examine how heat transport changes through the water column. Surface intensified warming has implications for temperature extremes, such as MHWs at the surface of the ocean in the euphotic zone. We focus on the primary inflow of heat across the western section (Fig. 1). Heat transports from the historical period show the surface-intensified flow of the EAUC next to the shelf of the North Island of New Zealand (Fig. 6a). The largest heat transport anomalies are observed within the top 250 m for all four sections.

Fig. 6.

Heat transport over all four sections surrounding the target region for mean historical (1995–2014) and future changes in SSP1–2.6, SSP2–4.5 and SSP3–7.0 (2050–2099) relative to the historical period: (a, e, i, m) Western; (b, f, j, n) Northern; (c, g, k, o) Eastern; and (d, h, l, p) Southern respectively. In the colour bar for the historical period, orange represents eastward and northward transport and purple represents westward and southward heat transport. The positive anomalies are increases in heat transport and negative anomalies decreases into and out of the study volume.

With the stronger greenhouse gas emission scenarios, the heat transport associated with the EAUC is projected to decline (Fig. 6e, i, m). Note in the SSP1–2.6 heat transports associated with the EAUC would be higher than the historical reference. The heat transport across the southern and eastern boundaries also weakens with increasing greenhouse gas emissions, while transport through the northern section increases. Integration of these anomalies reveals that declining EAUC heat transports are overcompensated by increased heat transport across the northern section with only mild declines across the southern section, leading to a heat transport convergence over this region in the top ~250 m. Please note that heat transports and projected changes are closely coupled to volume transports, with an r value between 0.98 and 0.99 and P values in the order of 10⁻⁷⁰ for each region and future scenario.

After investigating the spatial variations of the lateral heat fluxes into the target region, we turn our attention to the area-averaged temperature profile to investigate the cumulative impact of changes in these lateral heat fluxes (Fig. 7). We observe an ongoing warming of the top 250 m in all SSPs, reaching ~2°C in SSP1–2.6 and ~4°C in SSP3–7.0 by the end of the century. This depth range represents surface and subtropical mode waters in this region (Fernandez et al. 2018). Below this layer, between 250 and 700 m, a cooling signal emerges, which is stronger in high emission scenarios.

Fig. 7.

Area-averaged temperature anomaly profiles over the target region for (a) SSP1–2.6, (b) SSP2–4.5, and (c) SSP3–7.0 (2050–2099) relative to the historical period (1995–2014). Stratification changes over the target region quantified by the Brunt–Väisälä frequency for (d) SSP1–2.6, (e) SSP2–4.5, and (f) SSP3–7.0 relative to the historical period (1995–2014).

Changes in the temperature profile have impacts on the overall stratification in this region, which is shown in Fig. 7d–f. The stratification of the top 1000 m of the ocean is described by the Brunt–Väisälä frequency, which is a measure of vertical stability. The SSP1–2.6 scenario shows a weakening of the stratification over most of the century between 100 and 250 m, with positive anomalies above and below this depth. The stratification increases in the SSP2–4.5 scenario in the top 50 m in the first half of this century and down to 400 m in the second half of the century. The SSP3–7.0 scenario shows larger stratification changes than the lower emission scenarios; the stratification increases from 2030 in the top 100 m and extends down to 400 m in the second half of the century. Increased stratification in the upper layer of the ocean can trap heat near the surface, thereby intensifying the strength, duration and occurrence of MHWs (Oliver et al. 2021).

The historical simulation suggests annual mean mixed layer depths between 45 and 75 m over our target region for the historical simulation from 1995 to 2014, with deeper mixed layers to the south (Fig. 8a). The future projections show in Fig. 8b–d, these mixed layers becoming shallower by the end of the century by up to 15 m in the SSP2–4.5 scenario and up to 30 m in the SSP3–7.0 scenario. Negative anomalies mean the mixed layer depth is becoming shallower.

Fig. 8.

Mixed layer depths for (a) historical simulation (1995–2014), and anomalies in (b) SSP1–2.6, (c) SSP2–4.5, and (d) SSP3–7.0 over the period 2050–2099 relative to the historical period (1995–2014). Negative anomalies mean the mixed layer depth is becoming more shallow.

Finally, we investigate changes in wind stress because winds also influence oceanic stratification through turbulent mixing (Fig. 9). Stronger winds induce stronger mixing, which dissipates heat to deeper depths and prevents heat building up at the ocean surface. In the historical mean (Fig. 9a), we observe a band of low wind stress reaching from Australia to the north of New Zealand between 35 and 30°S. South of this region, wind stress increases. Future projections suggest a further decline in wind stress to the north and east of New Zealand, including over our target region. The pattern and magnitude of these anomalies is very similar between all three SSPs, which implies little sensitivity to greenhouse gas emissions. Ozone recovery is a possible cause for the changing wind pattern and is identical in each future scenario in NZESM. The projected wind decreases are consistent with the predicted future shoaling of the mixed layer (Fig. 8).

Fig. 9.

Wind stress magnitude for (a) historical simulation (1995–2014), and anomalies in (b) SSP1–2.6, (c) SSP2–4.5, and (d) SSP3–7.0 over the period 2050–2099 relative to the historical period (1995–2014). There is little difference between the three future scenarios.