Projections of Future Marine Heatwaves for the Oceans Around New Zealand Using New Zealand's Earth System Model

This paper investigates marine heatwave (MHW) characteristics in New Zealand's Earth System Model (NZESM) simulations for present-day conditions and how they are projected to change in the future in relation to anthropogenic greenhouse gas emissions. Three emission scenarios following the state of the art shared-socioeconomic pathways (SSPs, SSP1 2.6, SSP2 4.5, and SSP3 7.0) are each evaluated with a set of three ensemble members. These analyses are focused on the ocean around New Zealand, where NZESM captures boundary currents and mesoscale eddies, due to its high-resolution nested ocean grid. For present-day conditions, the model overestimates MHW intensity and underestimates the number of annual MHW days for subtropical waters, while some smaller positive biases are present in subantarctic waters compared to observations. Despite this, NZESM agrees with the observational pattern that more intense MHWs and more annual MHW days are found in subtropical waters compared to subantarctic waters. NZESM projects that MHW intensity will increase more strongly in subtropical waters compared to subantarctic waters, while the largest changes in annual MHW days are projected south of Australia and the Tasman Sea in the Subtropical Front (STF) frontal region, which suggests a southward shift of the STF under increased greenhouse gas emissions. Results using a high-emission scenario (SSP3 7.0) show an increase between 80 and 100% of median MHW intensities by the end of the century relative to the present-day for all analyzed coastal regions, and MHW conditions could become permanent year-round by the end of the century.


INTRODUCTION
New Zealand, with its two main islands, is surrounded by subtropical waters, which greatly impact its climate (Behrens et al., 2021). The prevalent westerly winds carry moist and warm air from the warm Tasman Sea toward New Zealand, where mountain ranges generate large rainfall upstream, and drier and warmer conditions downstream of the mountains, due to the Foehn effect (e.g., Wratt et al., 2000). Climate change is affecting New Zealand's oceans, where a southward expanding Subtropical Gyre and southward shifting Subtropical Front causes warmer sea-surface temperatures (SST) (Cai et al., 2005;Roemmich et al., 2007;Shears and Bowen, 2017;Sutton and Bowen, 2019;Behrens et al., 2021). In addition to the underlying continuous warming, several climate extremes such as marine heatwaves (MHWs) have been recorded around New Zealand over the past decades (Oliver et al., 2017;Salinger et al., 2019). The coastal regions around the Tasman Sea host a rich and diverse marine ecosystem, aquaculture facilities and commercial and recreational fishing grounds (Edgar and Barrett, 1999;Crawford, 2003;Gordon et al., 2010;Johnson et al., 2011;Townsend et al., 2018;Stenton-Dozey et al., 2021). All of these are impacted by extreme events, and hence new research analyzing and projecting changes in MHWs is beneficial to quantify future risks (Oliver and Holbrook, 2018;Chiswell and O'Callaghan, 2021). Blocking high pressure systems and changes in ocean heat content due to changes in the oceanic transports have been identified as the main drivers for MHWs Holbrook et al., 2019;Salinger et al., 2019Salinger et al., , 2020Li et al., 2020). With a warming ocean, caused by further increasing greenhouse gas (GHG) emissions, MHWs are expected to become more intense and longer compared to the present-day (Oliver et al., 2021). Trends in MHWs are not spatially uniform and regions with larger projected changes have been identified; usually where oceanic currents, such as the East Australian Current (EAC), are projected to change (Oliver et al., 2015;Gupta et al., 2021).
Updated projections of future MHWs from the new Coupled Model Intercomparison Project 6 (CMIP6) (Plecha and Soares, 2020;Qiu et al., 2021) against the older CMIP5  increase the confidence in projected future MHW conditions, due model improvements and refined future scenarios. These largely global assessments, however, do not allow for a focussed regional assessment of MHWs. Here, details in the regional oceanography, relating to oceanic currents and water mass properties play an important role in how MHWs arise and how they influence the ecosystem, in particular coastal regions. In this study we assess MHWs in the ocean around New Zealand in the New Zealand Earth System Model [NZESM, Williams et al., 2016] for present-day conditions and future conditions using the CMIP6 scenarios. The aim of this study is to generate a spatially detailed understanding of how MHWs could change around New Zealand in relation to GHG emissions. Recent studies have shown that UKESM (Kuhlbrodt et al., 2018), the basis of NZESM, has a high climate sensitivity to GHG emissions, with an equilibrium climate sensitivity (ECS) and transient climate response (TCR) of 5.3 and 2.8 • C, respectively (Meehl et al., 2020). Both values exceed the considered likely ranges of 1.5-4.5 • C for ECS and 1-2.5 • C for TCR (Nijsse et al., 2020), reflecting that NZESM will show a stronger surface warming response to GHG emissions globally.
The paper is structured as follows. Section Methods describes the model and methods. Section Results presents the main results. Here, Sections Regional MHW Intensity, Regional Annual MHW Days, MHW Duration, and MHW Seasonality provide a detailed overview of how MHW characteristics change for coastal regions around New Zealand and over the Tasman Sea. Section Discussion and Summary discusses and summarizes the findings of this paper.

METHODS
NZESM is a fully coupled earth system model based on UKESM, which uses NEMO (Madec et al., 2017) and MEDUSA (Yool et al., 2013) to simulate ocean physics and ocean biogeochemistry, the Unified Model (Walters et al., 2019) and Jules (Walters et al., 2019) for the atmosphere and land surface processes, and CICE (Rae et al., 2015;Hunke et al., 2017) for sea-ice. NZESM has been refined to capture the oceanic circulation for the New Zealand region better than its parent model, by embedding a 1/5 • resolution nested region in the ocean around New Zealand  and by incorporating a solar-cycle dependent photolysis rate (Dennison et al., 2019). The ocean grid cells around New Zealand are ∼17 km wide and allow oceanic mesoscale processes to be resolved without the use of Gent and McWilliams (1990) parameterization. More technical details of NZESM can be found in Behrens et al. (2020). Resolving these small-scale processes and improving transport of boundary currents in NZESM has led to reduced model biases in temperature and salinity around New Zealand. However, more realistic transports, in particular by the East Australian Current and improved circulation around New Zealand in combination with the elevated climate sensitivity of UKESM, have led to an exacerbated warm SST bias in the Southern Ocean .
Using NZESM, three historical ensemble simulations have been performed, started from different initial conditions and simulating the period from 1950 to 2014. Different start dumps (e.g., different years) of historical UKESM simulations were used to generate slightly different initial conditions for the historical simulations in NZESM. These three simulations have each been continued, from 2015 to 2100, following three climate change scenarios: SSP1 2.6, SSP2 4.5 and SSP3 7.0 (O'Neill et al., 2016), hereafter SSP 2.6, SSP 4.5 and SSP 7.0, which are part of the Tier 1 set of simulations within ScenarioMIP (O'Neill et al., 2016). These projections are intended to span a realistic range of possible future climates from aggressive reduction of GHG emissions (SSP 2.6) up to continued increases in GHG emissions (SSP 7.0). The number of scenarios and number of ensemble members in NZESM is constrained by available computational resources.
MHWs in NZESM and observationally-based NOAA OI-SST (https://psl.noaa.gov/data/gridded/data.noaa.oisst.v2.highres. html) have been diagnosed using key numerical guidelines outlined by Hobday et al. (2016). The widely-used Hobday et al. (2016) definition to detect MHWs is based on daily SSTs, a daily SST climatology and a 90th percentile threshold which needs to be exceeded for a minimum of five consecutive days to identify a MHW. The absence of daily SST data in NZESM prevented us to follow this approach. A pragmatic choice was taken to overcome this issue by using five daily averaged SST data from NZESM and NOAA OI-SST instead. These five-daily averaged SSTs were used following the Hobday et al. (2016) guidelines to construct a 5-daily climatology and to obtain the 90th percentile thresholds. A MHW was then identified when the actual SSTs exceeded the 90th percentile. No minimum duration limit was applied since the data represents mean SSTs over a 5-day period. The climatology was generated over a 30-year period from 1983 to 2012, following Oliver et al. (2018b). No additional smoothing has been applied to the baseline climatology, since the 5-daily averages reduced the level of high-frequency variability compared to daily SSTs (Supplementary Figures S1-S3). Overall our method detects the same features but is less restrictive than the Hobday et al. (2016) definition which affects the MHW characteristics. Using 5-daily averages instead of daily SST data reduces the median MHW intensities by about 0.4 to 0.8 • C but increases the median annual MHW days by 20-40 days.
In this paper the MHW intensity, annual MHW days, the length of MHWs and the timing during the year when MHWs occurs have been evaluated. The MHW intensity is the SST difference between the actual SST and the mean of the reference climatology. The annual MHW days is the number of days per year when the MHW criteria are met. The length of a MHW defines for how many consecutive days the MHW criteria is met. Since 5-daily mean SST data has been used in this study, the interval steps for annual MHW days and MHW length for an individual simulation of NZESM is a multiple of 5 days. For each simulation of NZESM these diagnostics were calculated separately.
In order to evaluate these MHW diagnostics and measure how they evolve in the future, 20-year periods have been defined. The first period spans the last 20 years of the historical period (1995-2014, hereafter considered as present-day conditions) and has been compared to 2040-2059 ("middle of the century") and 2080-2099 ("end of the century"). The three ensemble members in each SSP cover these 20-year periods with, in total, 60 years of 5-daily MHW statistics. Over these 20-year periods medians were calculated in NZESM and observations to describe the overall MHW ensemble characteristics and to explore the model range (see Supplementary Figures S5, S6 for results of individual NZESM simulations). To test if future projected MHW statistics change significantly in respect to the historical reference period a two sample and two-sided Kolmogorov-Smirnov test (Massey Jr, 1951) has been performed, where p-values of <0.01 have been used to identify significant change from the historic baseline. The significance test has been applied across the 60 year-long sample from all three ensemble members.
Making use of the fine oceanic model mesh of NZESM, a more detailed evaluation of MHW characteristics has been conducted for the Tasman Sea and three coastal regions around Tasmania and the North and South Islands of New Zealand. The Tasman Sea is defined as the oceanic region between 147 to 173 • E and 43 to 31 • S and has been identified together with the coastal water of Tasmania, where sea surface temperatures increase rapidly (Oliver et al., 2013(Oliver et al., , 2018a. The three coastal regions are defined by grid boxes within 50 km from the coast and water depths of <100 m. Using this definition the coastal region around Tasmania, North Island, South Island are characterized by 127, 189, and 219 grid boxes, respectively (Supplementary Figure S4).

Median MHW Intensity
Intense MHWs (>2.5 • C) in NZESM (Figure 1a) are visible in the EAC Extension, around Tasmania and western parts of the Tasman Sea (orange box). In addition, strong median MHW intensities (>2 • C) are also present to the North-East and East of the North Island of New Zealand and south of the Chatham Rise. MHW intensity in the subpolar waters, south of the STF, is in general lower compared to subtropical waters north of the STF.
The observed MHW intensity pattern from NOAA OI-SST (Figure 1d) shows similarities to the modeled pattern, while the overall modeled intensities are around 0.75-1.75 • C higher (Figure 1g) than NOAA OI-SST. The largest biases (>1.25 • C) are found in the EAC Extension, around Tasmania, north of the North Island and south of the Chatham Rise. In the observations the Tasman Sea, the coastal waters of Tasmania and region east of New Zealand are characterized by elevated MHW intensities (>1.5 • C). More intense MHWs are also found to the northeast and east of the North Island of New Zealand, with intensities exceeding 2 • C in some locations. Observed MHW intensities also show, as with the model, higher intensities in the subtropical waters compared to the subantarctic waters. Overall, the model shows for most regions a positive model bias in MHW intensity, caused by a larger warming trend during this historical period than seen in the observations. Note that the UKESM shows a predominantly negative bias for most of this region, with exception to North East Tasman Sea and to the North East Tasmania (Supplementary Figure S7).
Projections show a further intensification of MHW in the future (Figure 1), relative to the 1995-2014 period. In SSP 2.6, MHW intensity increases are largest (>0.5 • C) in the southern Tasman Sea, while intensities in the remaining subtropical waters increase between 0.25 and 0.5 • C for the period 2040-2059 (Figure 1b). For 2080 to 2099 (Figure 1c), the MHW intensities continue to increase further, with the Tasman Sea showing an increase between 0.5 and 0.75 • C, reflecting an increase of around 40%. This is also true for the region north and east of the North Island of New Zealand and around the Chatham Rise.
In SSP 4.0 (Figures 1e,f) the MHW intensity response is stronger than for SSP 2.6. Between 2040 and 2059 the ocean around New Zealand and over the Tasman Sea shows MHW intensities between 0.5 and 0.75 • C relative to the historic period. This intensification continues through to 2080-2099, where the Tasman Sea, the region around the North Island of New Zealand and the Chatham Rise, show median MHW intensities well above 1 • C, which corresponds to an increase of about 50%.
In SSP 7.0 the anomalies for 2040-2059 (Figure 1h) are very similar to SSP 4.5 during that time. But the intensity anomalies increase by around 1.5 • C for the later period 2080-2099 (Figure 1i) above anomalies seen in SSP 4.5. The ocean around Tasmania, around the North Island of New Zealand and around the Chatham Rise experience intensities above 1.75 • C and the remaining subtropical region experiences intensities above 1 • C, which reflects an increase of roughly 100% compared to the historical period. Overall, these diagnostics show that with increased GHG emissions MHWs become more intense, but not in a spatially uniform pattern. Subtropical waters, north of ∼45 • S, experience a stronger warming than subantarctic waters, and therefore more intense MHWs.

Median Annual MHW Days
The NZESM ensemble (Figure 2a) shows the largest number of annual MHW days (>50 days) in the southern part of the Tasman Sea and around Tasmania. Lower numbers (<40 days) of annual MHW days are seen north of the Tasman Front, east of New Zealand and in the sub-Antarctic waters. The observations show a similar pattern for the Tasman Sea but more MHW days east of New Zealand (Figure 2d). The model underestimates the number of annual MHW days (20-30 days) between 24 • S and the STF west of New Zealand, but up to 70 days east of New Zealand (Figure 2g). At the same time annual MHW days are overestimated (20-30 days) in subantarctic waters, south of the STF. Note, UKESM overestimates (>50 days) annual MHW days over the Tasman Sea region and toward the Subtropical Front (Supplementary Figure S8).
In SSP 2.6 for the period 2040 to 2059 (Figure 2b), the number of annual MHW days increases by at least 25 days in subtropical waters, with the southern Tasman Sea showing more than 75 MHW days per year, an increase of at least 75% to the historical period. For the period 2080-2099 the number of annual MHW days increases by about 50 days (Figure 2c), but there is a larger increase in the Tasman Sea and east of the Chatham Rise.
In SSP 4.5 the number of annual MHW days for 2040-2059 increases by around 25 days compared to SSP 2.6 (Figure 2e), with subtropical regions experiencing an increase of more than 75 annual MHW days in SSP 4.5, compared to the historical period. In particular, the southern Tasman Sea and the region south of it shows the largest increase in annual MHW days, of more than 125 days, compared to the historical period. This represents an approximate doubling compared to SSP 2.6. For 2080-2099 this increase in annual MHW days continues further and subtropical waters show more than 125 additional annual MHW days (Figure 2f). The southern Tasman Sea and the region east of the Chatham Rise experience more than 150 additional annual MHWs days at that point in time.
The number of annual MHW days in SSP 7.0 for 2040-2059 (Figure 2h) is very similar to SSP 4.5 but with an overall slight increase in annual MHW days. However, for 2080-2099 the number of annual MHW days increases substantially (Figure 2i) compared to SSP 4.0 and to the earlier time period. At that point subtropical waters experience between 150 and 200 additional annual MHW days (∼200% increase), and the regions of the southern Tasman Sea and east of Chatham Rise experience more than 225 additional annual MHW days, compared to the historical period.

Regional MHW Intensity
Around the North Island (1st row, Figure 3), NZESM (blue) overestimates MHW intensities on average by about 0.7 • C compared to NOAA OI-SST (black), for the historical period 1995-2014. NZESM does not capture the low MHW intensity (<0.25 • C), which causes the bimodal distribution in NOAA OI-SST. In the future scenarios, the historical distribution shifts to more intense MHWs, which leads to a broadening of the distribution. In SSP 7.0 the averaged MHW intensity reaches 4 • C, which reflects nearly a doubling compared to the historical period. The distributions for the coastal region around the South Island (2nd row) are very similar to those around the North Island. The model results show a better agreement with the NOAA OI-SST than around the North Island, and mean MHW intensities are only 0.6 • C higher, but the model does not show the bimodal observed distribution. The mean MHW intensities are lower compared to the coastal region around the North Island, suggesting a lower warming over the historical period. However, in SSP 7.0 for 2080-2099, mean MHW intensities are projected to reach 3.4 • C, compared to 1.9 • C for the historical period and reflecting an increase of 80% compared to historical intensities.
For the Tasman Sea (3rd row), the distributions are very similar to the coastal region of the North Island. The mean modeled MHW intensities are about 0.3 • C higher compared to NOAA OI-SST. The future distributions follow very much the evolution for the coastal region of the North Island and show mean MHW intensities of 2.9 • C in SSP 7.0 for 2080-2099, which is an increase of 80% to the historical period. The coastal region around Tasmania (bottom row) follows the Tasman Sea distribution with mean MHW intensity 0.7 • C higher than observations. The coastal region around Tasmania shows a similar warming as the Tasman Sea and South Island, with an increase of 75% by the end of century in SSP 7.0. Overall, the coastal region around the North Island experiences the largest absolute increase in future MHW intensity, compared to the other regions.

Regional Annual MHW Days
While the MHW intensity showed relatively little regional dependencies, the number of annual MHW days and probability distribution varies substantially between regions (Figure 4). NZESM underestimates the number of annual MHW days on average by about 17 days compared to NOAA OI-SST for the coastal region around the North Island (Figure 4a). Here, NOAA OI-SST shows on average 53 annual MHW days. The future simulations (Figures 4b-d) show increases in the number of annual MHW days, which reach 168 days on average in SSP 7.0 for 2080-2099-an increase of 360%. By that time the coastal region will experience at least 100 MHW days every year. NZESM shows a good match compared to the NOAA OI-SST for the coastal region around the South Island NZESM over the historical period (Figure 4e) with a mean model bias of about −5 days. The future projections show a broadening of the distribution toward a larger number of annual MHW days (Figures 4f-h). In SSP 7.0 for 2080-2099 this region experiences on average 174 annual MHW days-an increase of 460% to the historical period. For the Tasman Sea NZESM underestimates the number of MHW days by about 18 days on average (Figure 4i). The future projections (Figures 4j-l) show similar behavior to the coastal North Island, but with an increased probability of higher numbers (>200 days) of annual MHW days in SSP 7.0 for 2080-2099. In this scenario annual MHW days increase by 400% on average. The observations around Tasmania show a peak around 100 annual MHW days, which is not captured by NZESM, which shows a bias of −28 days on average. In SSP 7.0 for 2080-2099 this regions experiences 172 annual MHW days on average-an increase of 220% to the historical period.

MHW Duration
For the coastal region around the North Island NZESM shows a good agreement with NOAA OI-SST (Figure 5a), with the largest probability for MHW to last between 5 and 10 days. In the future projections a bimodal distribution emerges (Figures 5b-d). Short MHWs with 5-10 days duration are still most likely, but the secondary peak is showing MHWs to last between 100 and 250 days in SSP 7.0 for 2080-2099. The cause of this long-duration peak will be discussed in the next section. The duration of MHWs for the coastal region around the South Island (Figure 5e) aligns well between NZESM and NOAA OI-SST, although the likelihood for short MHWs (<10 days) is underestimated and-as for the coastal region for the North Island, the distribution develops a bimodal structure. In comparison to the North Island, the probabilities for the secondary peak are lower for the South Island and the distribution broader. The results for the Tasman Sea (Figures 5j-l) are very similar to those for the coastal region of the North Island, with modeled probabilities for short lasting MHWs (≤20 days) slightly overestimated. The coastal region around Tasmania (Figures 5m-p) largely follows the response of the Tasman Sea, but with a reduced probability for short lasting MHWs (≤30 days) and enhanced probability for longer MHWs (>200 days).
Overall, NZESM is able to simulate the duration of MHWs realistically and develops a bimodal distribution in the future with a proportion of MHWs projected to last longer than 100 days in SSP 7.0 by the end of the century. Relative increases are largest around New Zealand and over the Tasman Sea (>350%) in this scenario, with those around Tasmania nearer 200%.

MHW Seasonality
The modeled seasonality for the probability of MHWs (solid line) over all four regions agrees with the observed seasonal cycle (Figure 6)), with MHWs starting to emerge by December and vanishing by the middle of June, across all three coastal regions. Only the Tasman Sea shows a very low probability for MHWs during the winter for the historical period. The peak probability for MHWs is reached in middle of February, which marks the peak of summer in the southern hemisphere. That is also when median MHW intensity is the highest (dashed lines). The future projections show a distinct broadening of the probability peak, which reflects an extension of the season when MHWs occur, with highest intensities occur during February. This broadening in the MHW probability is more sensitive to the specific SSP scenario, while the response between all four regions is very similar. The results show that the extension of the MHW season is not uniformly split between the start and end of the season. In SSP 7.0 for 2080-2099 the MHW season starts in November about 1 month earlier, compared to the historical period (Figures 6d,h,l,p). However, probabilities at the end of season only gradually decline, with low values reached between July and September. That reflects a disproportional extension of the MHW season of 2-4 months at the end of the season. In SSP 7.0 by 2080-2099 all regions show MHWs during winter season, while probabilities are very low for the North Island and Tasmania. However, the other two regions suggest that some MHWs could persist for more than 365 consecutive days. Hints of that behavior were already seen in Figure 5, where these two regions showed elevated probabilities for MHWs to last more than 300 consecutive days. The median intensity for these MHW during winter seasons stays above 1.5 • C for the South Island, Tasman Sea and Tasmania while around 1 • C for the North Island. February median intensities vary between 4.8 and 5.5 • C between all regions.
The evolution of the seasonal profile explains the emergence of the bimodal distributions for MHW duration seen in Figure 5. In a warmer future, MHW lasting the full summer becomes virtually guaranteed, producing the long-duration peak. However, the other seasons are still not warm enough to produce long-persistent MHW events, instead producing the short-duration peak.

DISCUSSION AND SUMMARY
This study uses high-resolution (1/5 • ) output from the New Zealand Earth System Model to characterize present-day and future changes of MHW characteristics in the ocean around New Zealand, following a set of three socioeconomic pathways (SSP1 2.6, SSP2 4.5 and SSP3 7.0). NZESM and its parent model UKESM has a higher climate sensitivity (ECS of 5.3 • C and TCR of 2.8 • C) to greenhouse gases. The considered likely range is from 1.5 to 4.5 • C for ECS and 1 to 2.5 • C for TCR (Meehl et al., 2020;Nijsse et al., 2020;Tierney et al., 2020). The higher climate sensitivity of NZESM means that global mean surface temperature will increase more strongly with greenhouse gas emissions than considered realistic. Recent results suggest that in particular surface temperature changes in tropical regions are more directly impacted by the higher climate sensitivity than in the Southern Ocean (Huusko et al., 2021). Instead of a multimodel approach (Alexander et al., 2018;Holbrook et al., 2019;  Oliver et al., 2019;Hayashida et al., 2020;Plecha and Soares, 2020), an ensemble approach has been used in this regional study, with three ensemble members in each SSP, to explore the solution space of MHW projections in an eddy permitting system. Coarse models do underestimate variability in energetic regions, like western boundary currents and the Southern Ocean, where eddy activity and associated intrinsic variability (Penduff et al., 2011) is high and affect how they predict MHWs (Frölicher et al., 2018;Oliver et al., 2019;Hayashida et al., 2020).
Nevertheless, our overall findings align with previous work, which shows in general that MHWs will become more frequent and intense under future warming scenarios compared to the present day (Frölicher et al., 2018;Oliver et al., 2018bOliver et al., , 2019Hayashida et al., 2020;Sen Gupta et al., 2020). Our model results show that these changes vary regionally. Subtropical waters experience more frequent and more intense MHWs under present-day conditions compared to subantarctic waters south of the STF. Future projections show a larger warming in subtropical waters than in subantarctic waters, which causes MHW intensity to increase disproportionally between both regions. The projected pattern follows present-day observed temperature trends (Shears and Bowen, 2017;Sutton and Bowen, 2019) in this region. For annual MHW days the largest projected changes are observed south of Australia and of the Tasman Sea in the STF region, which points to a southward shift as a consequence of expanding subtropical gyres and changes in western boundary currents in the Southern Hemisphere (Oliver and Holbrook, 2014;Hu et al., 2015;Yang et al., 2016Yang et al., , 2020Gupta et al., 2021;Misra et al., 2021).
Due to the absence of daily SST data in NZESM, on which MHW detection is normally based on following the Hobday et al. (2016) definition, a pragmatic approach was taken by using 5-daily averaged SSTs instead. While both methods detect in general the same features, their MHW characteristic differ due the underlying data. Using 5-daily averaged SSTs data reduces the median MHW intensities between 0.4 and 0.8 • C but increases median annual MHW days between 20 and 40 days compared to daily SST data. This provides a guideline to place the results into context compared to research based on the Hobday et al. (2016) definition. Recent research has documented mean MHW intensities based on the Hobday et al. (2016) definition between 30 and 20 • S around 1.7 • C and mean annual MHW days around 10 days per year for 1995-2015 (Holbrook et al., 2022), which aligns with our comparison (see Supplementary Figure S1-S3). A further study using a CMIP6 ensemble showed a very similar pattern of projected change of annual MHW days and intensity around New Zealand (Qiu et al., 2021), providing confidence that our results are quantitatively similar despite being based on 5-daily averaged SSTs. Due to the fine-scale ocean grid of NZESM we attempted to characterize coastal (<50 km and <100 m deep) changes around Tasmania and around both Islands of New Zealand, because of the ecological importance of these regions. NZESM does resolve boundary currents and associated variability better than most global earth system models but still does not resolve and incorporate all key processes for these coastal regions (e.g., sub-mesoscale processes and tides). Despite that, it does provide a first closer look at how MHWs will impact these environments in the future, but further downscaling activities would be required to project these results onto local near-shore scales.
The regional diagnostics reveal that future MHW intensity increases most strongly in coastal waters around the North Island, and comparatively less around Tasmania, Tasman Sea and the South Island. However, the Tasman Sea and the coastal region around the South Island are regions where under high CO 2 emission scenarios MHWs can become year-round features relative to a fixed baseline climatology from 1983 to 2012. A moving baseline climatology would be useful when studying the ecosystem response to MHWs and assuming the ecosystem adapts quickly to a new baseline (Oliver et al., 2021). Since the focus here is not to investigate such ecosystem changes, it motivated the use of a fixed baseline and so a moving baseline climatology is out of scope of this study. Furthermore, our results show that the MHW season, which is centered around summer, will disproportionally extend into the autumn season in the future.
While NZESM has been purpose-built for earth system studies around New Zealand and the Southern Ocean, with a few previously identified model biases having been reduced , some biases still persist. The model overestimates MHW intensities and underestimates annual MHW days in the subtropical region under present day conditions. The MHW bias pattern between NZESM and UKESM differ, where UKESM shows larger biases around Tasmania and the North Island of New Zealand than NZESM (Supplementary Figures S7, S8). While some future changes in NZESM have been evaluated relative to the historical period, which would eliminate a timeinvariant bias, it remains unknown how the underlying bias in NZESM evolves with time and impacts future projections. Here, the comparison with other models can help to identify robust change, while acknowledging there is no model without a bias. However, this exercise is considered out of scope for this particular study. We note the Tasman Sea and region east of New Zealand has experienced intense MHWs in the summers of 2015/16 (Oliver et al., 2017), 2016/17 Salinger et al., 2019) and 2018/19 (Salinger et al., 2020). Due to our choice of the baseline period these events were not included in our observational records. However, simulated climate variability in NZESM does not align in time with observed climate variability so potentially some of these intermittent MHW events have been included in the model results, affecting the described model bias. A larger model ensemble (>3 members) would have been useful here but was not feasible due to computational constraints. The small ensemble size prevented us from looking at extreme MHW events (e.g., maximum intensities) and therefore predominantly median statistics have been applied to diagnose conservative and robust changes.
In summary, the model projections show that distributions of MHW intensity, MHW duration and annual MHW days will broaden to more intense, longer and more frequent MHWs in the future, with the amount of change dictated by the emission scenario. While median MHW intensity increases of about 2 • C in SSP3 7.0 for 2080-2099 appear relatively modest the likelihood of more severe MHWs (>6 • C) is increased, due to changes in the underlying probability distribution. Here a set of ensemble members in one model framework is helpful to explore the solution space in a non-linear, eddying system. It is a future challenge for models when approaching eddy resolving resolution to disentangle forced changes from intrinsic variability.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
EB performed the analyses and compiled the paper draft.
JW performed the NZESM simulations. GR, SR, JW, OM, and DS provided feedback on diagnostics and contributed to the paper draft. All authors contributed to the article and approved the submitted version.

FUNDING
NOAA High Resolution SST data provided by the NOAA/OAR/ESRL PSL, Boulder, Colorado, USA, from their website at https://psl.noaa.gov/. This project obtained funding from the Ministry of Business, Innovation and Employment through the Deep South National Science Challenge (C01X1902).

ACKNOWLEDGMENTS
We acknowledge the support of NeSI High Performance Computing Facility and its team. In addition, we acknowledge those involved in the development of UKESM and the MetOffice team. Furthermore, we acknowledge the time, effort, and constructive comments of all 4 reviewers to help us improve this manuscript. Finally, EB would like to acknowledge my partner, son, and daughter.