Forecasting extreme marine heat events in key aquaculture regions around New Zealand

Keywords: ACCESS-S, aquaculture, climate change, marine heatwave, model skill, ocean warming, seasonal forecasting, Tasman Sea.

The Western Pacific and Tasman Sea regions have been identified as climate hotspots (Law et al. 2018; Oliver et al. 2017). The rate of warming around New Zealand (NZ) and in the Tasman Sea is as much as 0.4°C per decade (Sutton and Bowen 2019), with subsurface waters warming at nearly four times the global average rate (Oliver et al. 2017). With a warming climate, extreme marine heat events are predicted to increase in both frequency and severity (Oliver et al. 2014; Rickard et al. 2016). This will have increasingly serious implications for regional marine ecosystems and the industries that rely on them (Chiswell and Sutton 2020).

Due to the latitudinal extent of Aotearoa (NZ) and its mostly narrow continental shelf regions (Stevens et al. 2019), marine temperatures in coastal waters are strongly influenced by a combination of major large scale oceanic fronts and their dynamics (table 1 of Chiswell et al. 2015), and processes such as wind stress and fluctuations in the Tasman Front in the Tasman Sea and the wider Western Pacific (e.g. Behrens et al. 2020). Although the western seaboards of both the North Island and South Island face the complex and weakly dynamic Tasman Sea, the rest of the coast is heavily influenced by strongly dynamic waters. The Tasman Front and East Auckland Current strongly influence the east coast North Island waters (Fig. 1). To the south, the Subtropical Front controls the east coast South Island conditions. The central zone around Cook Strait acts as a confluence of all of these systems and is highly variable (Stevens et al. 2019).

Fig. 1.  (a) Map of the NZ region showing the main oceanographic features (Subantarctic Water (SAW), Subtropical Water (STW), Tasman Front (TF), East Auckland Current (EAuC), Southland Current (SC), Fiordland Current (FC), Westland Current (WC)) and the selected three regions of (b) Hauraki Gulf (Cape Reinga (CR)), (c) Western Cook Strait (where the strait itself is identified as CS and the d’Urville Current (dC)) and (d) Foveaux Strait (FS) (between the South Island and Stewart Island Rakiura (SIR)). Shading in subplots (b)–(d) indicate the mean correlations across each month in the year between the regionally averaged sea surface temperature anomaly (SSTA) index for each region and SSTA values within the region.

Several well-documented marine heatwaves (MHWs) have occurred in recent summers in the Tasman Sea and NZ coastal waters (Salinger et al. 2019; Salinger et al. 2020). Coastal MHWs can compound other stressors, leading to detrimental impacts for local ecosystems and industries. One example of this in the NZ context is the combined effects of MHWs and poor water clarity on giant kelp in coastal regions (Tait et al. 2021). Hobday et al. (2016a) define a MHW as five consecutive days above the daily 90th percentile, meaning events can occur year round. Based on this definition, three MHWs have occurred in the Tasman Sea in recent years in the summers of 2015/2016, 2017/2018 and 2018/2019 (Oliver et al. 2017; Perkins-Kirkpatrick et al. 2019; Holbrook et al. 2020). In particular, the 2015/2016 MHW in the Western Tasman lasted 250 days, with sea surface temperatures (SSTs) reaching nearly 3°C above climatology (Oliver et al. 2017). These events were associated with low cloud cover and low wind conditions, so that the heat remained confined to the upper ocean, which then gave rise to a MHW signal in the SSTs (Salinger et al. 2019).

Behrens et al. (2019) examined heat content in the upper 2000 m of the Tasman and identified different behaviours for the upper 250, 250–750 and >750 m. In particular, they found that the upper 250 m exhibited a cool phase during the 1990s but then entered a warm phase from 1990 to 2002 and from 2012 to the present. Elzahaby et al. 2021 found that MHWs in the Western Boundary Current (WBC) jet of the East Australia Current are predominantly driven by air–sea heatflux, whereas in the WBC extension, they are advection-driven. They also note that the deepest and longest MHWs are advection-driven and are more prevalent in autumn and winter, whereas air–sea heatflux-driven MHWs are shallower and occur predominantly in summer. This highlights the variety of ‘flavours’ of MHWs that can occur in the Tasman Sea: both shallow, atmospherically driven events, and deeper events driven by ocean currents, the nature of which then impact MHW predictability. Both flavours of MHWs have impacted both marine and terrestrial ecosystems and related human activity, such as fishing and aquaculture.

The 2015/2016 Tasman MHW event had significant ecological impacts on giant kelp (Macrocystis pyrifera; Tait et al. 2021) and a range of wild and cultured species, including oysters, paua (abalone), salmon and kingfish (Oliver et al. 2017) within the Tasman Sea and along the Australian and NZ coasts. The 2017/2018 event on the southeastern seaboard saw substantial mortality of kelp (e.g. Duvillaea spp., Thomsen et al. 2019; Macrocystis pyrifera, Salinger et al. 2019) and salmon (Salinger et al. 2019). Chiswell and O’Callaghan (2021) also examined the impacts of MHWs on upwelling and primary production along the South Island west coast, which can in turn affect productivity. These types of extreme ocean heat events have significant implications for aquaculture and fisheries, which are major industries in NZ, contributing ~NZ$1.1 billion to the economy (Stats 2017). King salmon (Oncorhynchus tschawytscha, also known as Chinook salmon) are the largest type of Pacific salmon and the only salmon species farmed in NZ. Optimal water temperatures for Pacific salmon range from 10 to 17°C, with elevated temperatures being related to enhanced mortality due to disease (Brosnahan et al. 2019). Shellfish are also vulnerable to thermal stress, including green-lipped mussels (Perna canaliculus) and oysters (Ostrea chilensis) (Sorte et al. 2019) species, which are farmed in NZ. As with elsewhere around the globe, these industries have been impacted by MHWs in recent years, leading to increased disease, fish kills, losses in productivity and, ultimately, reduced profit. The impacts are not solely caused by direct exposure to elevated temperatures; increased temperatures can drive enhanced stratification, which, in turn, leads to reduced oxygen or nutrient availability (e.g. Chiswell and Sutton 2020).

Given sufficient warning, these industries can mitigate some of the effects of marine heat events. Forecasts on multi-week to seasonal timescales have been shown to benefit the proactive management of marine industries (Hobday et al. 2016b; Tommasi et al. 2017). These have included commercial fisheries (e.g. southern bluefin tuna; Hobday et al. 2011; Eveson et al. 2015), recreational fisheries (Brodie et al. 2017), aquaculture (salmon, prawns; Spillman and Hobday 2014; Spillman et al. 2015) and coral reef management (Spillman 2011; Smith and Spillman 2019). These forecasts are important for providing advance warning of suboptimal and extreme ocean temperature conditions, and thus, they allow time for industries to implement proactive response plans to minimise and mitigate the impacts of such events.

The Australian Bureau of Meteorology produces operational global seasonal forecasts out to 6 months from the coupled ocean–atmosphere ensemble prediction system, the Australian Community Climate and Earth-System Simulator-Seasonal v1.0 (ACCESS-S1; Hudson et al. 2017). Seasonal forecast skill for SST and heat content to 300 m around NZ has been assessed and was found to be higher for inshore regions (depth < 300 m) than offshore for 2 weeks to 6 months into the future (de Burgh-Day et al. 2019). Model SST skill is higher to the west of NZ where currents are relatively slow moving and well separated, and to the north where skill is derived from well-simulated teleconnections with the Eastern Pacific and the El Niño Southern Oscillation (ENSO) (de Burgh-Day et al. 2019). Conversely, to the east of NZ, SST skill is lower due to complex dynamics, including the convergence of a number of currents and the meeting of the Subtropical and Subantarctic fronts (see Fig. 1 for major currents in the region). Skill is also reduced in the south where ACCESS-S1 is overly reactive to La Niña conditions and is impacted by a warm bias in the Southern Ocean (de Burgh-Day et al. 2019).

In this study, we assessed model forecast skill in three key coastal regions which either contain, or directly influence, regions of significant aquaculture activity (NZGAS 2019): Hauraki Gulf, Western Cook Strait and Foveaux Strait (Fig. 1). The regions defined for this study were intended to give an indication of the skill of ACCESS-S1 for aquaculture industries around NZ. The regions span a diversity of ocean structures and dynamics, but it was not possible to investigate the impact on skill on a smaller scale due to the resolution limitations of the ~25 km grid of the ocean component of the model. The extent of each region was therefore intended to be as targeted as possible without being affected by the limitations of model resolution. The Hauraki Gulf sits on the northern edge of the North Island and supports mussel and oyster farming industries as well as proposed king fish farming. Here the Western Cook Strait region was assumed to directly influence conditions in the Marlborough Sounds, a network of ancient sunken river valleys, where extensive salmon and mussel farming occurs due to nutrient rich and quiescent waters (Zeldis et al. 2013; Chiswell et al. 2017). Foveaux Strait hosts a dredge oyster fishery industry and is the primary upstream supply of clean water for nearby embayments that support aquaculture industries (e.g. salmon aquaculture in Big Glory Bay, Stewart Island Rakiura). It is also a potential site for further development of aquaculture industries in the future (Camara and Symonds 2014) as warming trends continue to shift habitat suitability for various species. First, we determined the ensemble mean forecast skill in these three key regions as well as the skill for extreme heat events (greater than the 90th percentile). Secondly, we considered the utility of real-time ACCESS-S1 forecasts in providing advance warning of extreme marine heat events around NZ in the context of the Tasman Sea marine heatwave event in the 2018/2019 summer.

Daily observed satellite SST values for 1990–2012 and 2017–2018 from the Reynolds Optimum Interpolation Sea Surface Temperature V2 (OISST-V2; for 1990–2012) and V2.1 (OISST-V2.1; for 2017–2018) daily Advanced Very High-Resolution Radiometer (AVHRR) 0.25° analysis were used for forecast skill assessment (Reynolds and Smith 1994; Reynolds et al. 2002). Observed daily SST values were first averaged to create monthly means and then interpolated onto the ACCESS-S1 ocean model grid. SST anomalies (SSTAs) were calculated by removing the appropriate monthly climatology, which is the long term monthly mean over the period 1990–2012. Persistence forecasts were also generated by persisting observed SSTA from the month prior to the model forecast start date out to 6 months into the future (see Spillman and Alves 2009).

SSTA indices were calculated by regionally averaging observed SSTA values within each of the three regions of interest: Hauraki Gulf, Western Cook Strait and Foveaux Strait (Fig. 1). The vertices of these regions, as defined for this study, are given in Table 1. Observed Pearson correlation coefficient values were computed between each SSTA index and the individual grid cells within the region using all months in the period 1990–2012. The SSTA indices had high Pearson correlation coefficient values for all grid cells in the corresponding region (r > 0.7 for all grid cells in all regions; Fig. 1b–d) and so were deemed representative of local SSTA conditions. A fourth SSTA index was also created for an area of the Tasman Sea (see Fig. 1, Table 1), chosen as that most severely impacted during the 2018/2019 MHW, to assess the relationship between the wider region and the three aquaculture subregions.

Table 1.  Vertices defining each the regions of interest for this study; Hauraki Gulf, Western Cook Strait and Foveaux Strait, and a segment of the Tasman Sea.

ACCESS-S1 was the Australian Bureau of Meteorology’s operational seasonal ensemble prediction system until October 2021 (Hudson et al. 2017), when it was upgraded to ACCESS-S2. ACCESS-S1 was developed in collaboration with the UK Met Office (UKMO). It utilises the UKMO Global Coupled model 2.0, which consists of the UKMO Global Atmosphere 6.0 model (Williams et al. 2015); the European ocean model, Nucleus for European Modelling of the Ocean (Madec 2008) and the Los Alamos sea-ice model (Hunke and Lipscomb 2004); together with the Joint UK Land Environment Simulator land surface model (Walters et al. 2017). The ocean model has a grid resolution of approximately 25 × 25 km in the horizontal and 75 vertical layers to a total depth of 6000 m, starting at 1-m thickness at the surface.

The ACCESS-S1 system included a set of retrospective forecasts (hindcasts) over a 23-year period (1990–2012), comprising of an 11-member ensemble on the 1st, 9th, 17th and 26th of every month. In this study, the hindcasts generated on the 1st of each month were used. This hindcast set was used to bias correct forecasts and remove any model drift with lead time as well as to assess model skill in predicting past events. The model was run 11 times for each start date in the hindcast period to give an ensemble of 11 forecasts, using a set of perturbed initial conditions (Hudson et al. 2017). Ensemble forecasts seek to account for uncertainty in initial conditions and model error by producing a distribution of possible forecast outcomes that encapsulate these sources of error (Hudson et al. 2017). A small ensemble spread suggests a confident forecast, whereas a large spread indicates a greater uncertainty in the forecast, and it should be interpreted with more caution. The ensemble mean was computed by averaging the 11 ensemble members for each start date and lead time. Monthly SST values and climatologies were then created by averaging over the 1990–2012 hindcast period for each start month. For each start date the appropriate monthly climatology was subtracted from the corresponding hindcast data to produce SSTAs. Area-averaged SSTA indices were calculated for Hauraki Gulf, Western Cook Strait and Foveaux Strait regions (see Fig. 1).

The real-time ACCESS-S1 prediction system was made operational in late 2018 (Smith and Spillman 2020) and ran in real-time until late 2021. For seasonal timescale forecasts, ACCESS-S1 generated 11 ensemble members out to 6 months every day, which were then combined with forecasts from the previous 8 days to produce a larger, time-lagged 99-member ensemble (see Hudson et al. 2017). An ensemble mean was produced by averaging the 99 ensemble members, and data were averaged along lead times to give monthly means. The time-lagged 99-member ensembles starting on 1 January, 1 February and 1 March 2019 were used in this study.

SSTA indices for the Hauraki Gulf, Western Cook Strait and Foveaux Strait regions were created by area averaging the values in all grid cells within each region. Pearson correlation coefficients were then computed between the model hindcast ensemble mean monthly SSTA index and the observed monthly SSTA index for each region for 1990–2012. Correlations were calculated between persistence SSTA forecasts and observed regional SSTA in the same manner, providing a minimum skill benchmark (Spillman and Alves 2009). For these correlations, values below −0.42 and above 0.42 are statistically significant (two-tailed t-test, n = 23). Additionally, the standard deviation of the area-averaged SSTA in each region was calculated for all model ensemble members. The standard deviation across the ensemble was then obtained by averaging the variances and taking the square root of the result. Corresponding observed SSTA standard deviations were also calculated.

In order to achieve a large enough sample size, it was necessary to include all months and ocean grid cells in each region for the probabilistic skill assessment. This is because when categorising extreme events in a distribution, if the sample taken from the distribution is small, it is possible that few or no instances of some extreme event categories will occur in the sample. In all three regions, the correlations between the grid cells and the regional SSTA index is high (>0.84) (Fig. 1b–d). The model SSTA 90th percentile (top decile) threshold was calculated for each year/start date/lead time combination in the hindcast using a cross-validation method: for a given start date and lead time (0–5 months), SSTA values from all hindcast ensemble members and years were pooled, omitting the year/start date combination for which the 90th percentile was being computed. The 90th percentile was then computed for this data by ranking all samples and determining the SSTA value for which 90% of the samples were below this value and 10% were above it. This was repeated for every year/start date/lead time combination, always leaving out the start date and year the percentile was being computed for comparison with. The 90th percentile threshold for observed SSTA values was calculated for each month and year in the hindcast period using the same cross-validation method.

A multi-category contingency table was then calculated to assess the skill in predicting above the 90th percentile in Hauraki Gulf, Western Cook Strait and Foveaux Strait using the following method: If X percent or more forecast ensemble members correctly predicted SSTA values above (below) the 90th percentile, then this was counted as a hit (correct negative). If X percent or more forecast ensemble members incorrectly predicted SSTA values above (below) the 90th percentile, then this was counted as a false alarm (miss). Here X is a probability threshold, and this process was repeated for values of X ranging between 0 and 100%. The multi-category contingency table was then used to calculate Receiver Operating Characteristic (ROC) curves and ROC areas (A). Reliability Diagrams and Brier Skill Scores (BSS) were also computed.

The ROC curve determines the ability of the forecast to discriminate between events and non-events and gives an indication of how often an event occurs, given it has been forecast, for a specific forecast probability (Mason and Graham 1999). In other words, these curves answer the question, ‘for a given forecast probability threshold at which I say an event will occur, how often will I be correct and how often will it be a false alarm?’. ROC curves are calculated by comparing the hit rate (i.e. the fraction of observed events that were correctly forecast) against the false alarm rate (i.e. the fraction of observed non-events that were incorrectly forecast as events) for a range of probabilistic thresholds (Mason and Graham 1999). The area under the ROC curve (A) is used as a summary statistic, with a larger area indicating more skilful probabilistic forecasts (see Table 2). A perfect forecast (100% hit rate, 0% false alarm rate, A = 1) sits in the upper left corner of the axes. A forecast that sits along the 1:1 line is as good as a random guess of whether the event will occur or not (A = 0.5), and any line to the left of the 1:1 line indicates skill relative to a guess.

Table 2.  Summary scores for SSTA forecasts for Hauraki Gulf, Western Cook Strait and Foveaux Strait regions for all months in 1990–2012.

Reliability diagrams for events where SSTA >90th percentile were generated by binning forecast probabilities and plotting against the observed rate of occurrence of the event in each bin (Wilks 2006). This indicates how reliably a given forecast probability will reflect the true probability of the occurrence of an event (Wilks 2006). In other words, reliability diagrams answer the question ‘when an event is forecast to occur with probability X, what is the observed frequency with which that event actually occurs?’. The model has skill over climatology, for a given forecast probability and lead time, when the corresponding point sits within the grey shaded area of the reliability diagrams shown. The vertical boundary of the grey area represents a fixed climatological forecast, and the sloped boundary represents the division between a positive or negative BSS relative to climatology. The BSS is a measure of the model’s skill to predict the probability of an event occurring compared to a reference forecast (Mason and Stephenson 2008). Here, both the climatological frequency of the event (e.g. climatologically, there is a 10% chance of exceeding the 90th percentile) and a persistence forecast are used. The BSS values of one and zero indicate a perfect score and no skill relative to the reference respectively, and a negative value indicates a forecast worse than the reference (Wilks 2006).

Correlations between observed and model ensemble mean monthly SSTA indices for Hauraki Gulf, Western Cook Strait and Foveaux Strait for each start date and lead time are shown in Fig. 2. Correlations are generally higher at shorter lead times for all model start dates and all locations. The model has higher skill than persistence forecasts for the majority of start and target month combinations, with the notable exceptions of lead 0 forecasts starting in April–June and October for Hauraki Gulf and May, July and October for Western Cook Strait. Model skill is often higher for Hauraki Gulf at longer leads than the other two locations, but it is also lower than persistence skill on more occasions (31% of skill values, Fig. 2a). The lowest correlations are for forecasts of summer months January, February and March at all three locations, with correlations generally not statistically significant beyond a 1-month lead. Conversely, the highest skill for all three locations is for forecasts of winter months (June–August), with the exception of July forecast from 1 July for Western Cook Strait (Fig. 2b) and June and July forecasts from 1 May for Foveaux Strait (Fig. 2c). When forecasting winter months from summer (e.g. December and January start dates), correlations drop below the significance threshold initially for some months before increasing again for winter months at longer lead times.

Fig. 2.  Correlations for regional SSTA indices for (a) Hauraki Gulf (b) Western Cook Strait and (c) Foveaux Strait between model hindcasts and observations for 1990–2012. The x-axes show the start date for the hindcasts runs (i.e. hindcast runs beginning on the 1st of each month were assessed), and the y-axes show the lead times of the hindcast runs (i.e. lead times of 0–5 months were assessed). Hatching indicates where skill is not significantly significant (r < 0.42, n = 23, two-tailed t-test). Grey dots indicate where model skill is greater than persistence skill. The text in each square indicates which month is being forecast and the colour is the correlation value.

The model is able discriminate between heat events (greater than the 90th percentile) and non-events and has skill over random chance for all regions and lead times (Fig. 3). Skill reduces with increasing lead time for all three locations. The areas under the ROC curves (A) range from 0.77 to 0.88, 0.75 to 0.89 and 0.7 to 0.83 for lead times of 2, 1 and 0 months for Hauraki Gulf, Western Cook Strait and Foveaux Strait, respectively (Table 2).

Fig. 3.  ROC curves for forecasts of SSTA above the 90th percentile for (a) Hauraki Gulf, (b) Western Cook Strait and (c) Foveaux Strait for all start dates in 1990–2012. A perfect forecast (100% hit rate, 0% false alarm rate) sits in the upper left corner of the axes. A climatological forecast sits along the 1:1 dashed line, and any line to the upper left of this is skilful relative to climatology. Probability thresholds values at which an event is forecast to happen are noted in the circles along the individual curves. See Table 2 for area under the curve values (A).

The model has skill over climatology for high and low probability forecasts for SSTA >90th percentile at all lead times shown for the Hauraki Gulf and Western Cook Strait regions (Fig. 4a, b). This is indicated for a given forecast probability and lead time when the corresponding point sits within the grey shaded area of the reliability diagram. The exception is the uppermost forecast probability for lead 2 forecasts for both locations, most likely due to an insufficient sample size (see inset bar plots). For Foveaux Strait, the model is skilful only at shorter lead times for high and low probability forecasts (Fig. 4c). The forecast is overconfident for all regions and lead times, as indicated by the curves for each lead time being above the 1:1 line for low probabilities and below the line for high probabilities (Fig. 4a–c). An over-forecasting bias is also present for all lead times and regions, as indicated by the downward shift of the curves relative to the 1:1 line. BSS values referenced to persistence for lead times of 2 months to 0 months range from 0.17 to 0.37 for Hauraki Gulf, from −0.08 to 0.22 for Western Cook Strait and from 0.22 to 0.4 for Foveaux (Table 2). BSS values referenced against climatology are generally lower (though positive) for all locations, due to strong variability, at all lead times.

Fig. 4.  Reliability diagram for SSTA forecasts in the 90th percentile for all start dates in 1990–2012 for (a) Hauraki Gulf, (b) Western Cook Strait and (c) Foveaux Strait. Forecast probability distributions for each location are in shown in the inset bar plots.

Model SSTA indices over time for the three regions, for the period 1990–2012 and for lead times 0–2 months are shown with observations in Fig. 5. Observed standard deviations are similar across the three regions (0.56–0.62, Table 2). In the model, mean ensemble standard deviations for the three indexes are all higher than those of the observed values, with the highest values occurring at a 2-month lead (0.71–0.86, Table 2). Correlations between observed and model indexes are high for a lead of 0 months for all three locations (0.80–0.83, Table 2), and they decrease with lead while maintaining correlations exceeding 0.5 at a lead of 2 months. The majority of observed data points are captured by the ensemble spread (grey shading) at all locations and leads (Fig. 5).

Fig. 5.  Time series of observed and ensemble mean monthly SSTA indices for (a) Haruki Gulf, (b) Western Cook Strait and (c) Foveaux Strait for 1990–2012. The top, middle and bottom plots of each panel are for lead times of 0, 1 and 2 months, respectively. Ensemble spread is indicated by the grey shading.

Fig. 6 shows the forecast probabilities from the real-time system for SSTA >90th percentile for the period from January to March 2019 in the Tasman Sea and around NZ. At lead time of 0 months (first plot in Fig. 6a, second in Fig. 6b, third in Fig. 6c), the model captured the spatial extent of the extreme SSTs reasonably well and with high probabilities. The exception to this was around the South Island for February at a lead time of 0 months, where the model incorrectly forecast to have a high chance of SSTA >90th percentile. Forecasts issued on 1 January underestimated the extent of the event in February (1-month lead) and March (2-month lead), with a marked reduction in areas with high forecast probabilities of SSTA >90th percentile compared to that forecast for January at a lead time of 0 months. Similarly, the area of SSTA >90th percentile forecast from 1 February for March (1-month lead) is reduced compared to that forecast for February (0-month lead) and is less than was observed in March.

Fig. 6.  Forecast probabilities for monthly SSTA exceeding the historical 90th percentile for January, February and March 2019 for model forecasts issued on (row a) 1 January 2019, (row b) 1 February 2019 and (row c) 1 March 2019. The final row (d) shows where observed monthly values exceeded the historical 90th percentile (yes – red, no – grey). The historical reference period is 1990–2012.

The timescales of the teleconnections between the Tasman Sea and coastal regions of NZ provide insight into drivers of heat transport between these locations and, thus, the link between Tasman Sea MHWs and extreme marine heat events around NZ. This has implications for the ability to successfully forecast the impacts of these events, since atmospherically driven phenomena tend to have limited predictability at longer lead times compared to ocean-driven phenomena (Cowan et al. 2019). To understand these timescales, Pearson correlation coefficients were computed between a lagged observed Tasman Sea regionally averaged SSTA index (Tasman index) and the observed SSTA in each grid cell around NZ for January, February and March 1990–2012. The instantaneous correlation (Fig. 7, first column) between the Tasman index and the grid cells contributing to that index is high (mostly >0.6 for all months shown), indicating that this index is representative of conditions in this part of the Tasman. Overall, correlations between the Tasman index and the whole domain decrease with increasing lag. However, February Tasman index correlations remain over 0.6 for the majority of the Tasman index region and a section following the path of the Tasman Front around the northern tip of NZ in March. This is not unexpected, as the Tasman Sea supplies the water flowing around the North Cape of NZ and into the East Auckland Current (Chiswell et al. 2015). The low correlations to the southeast of NZ can be attributed to the incursion of Subantarctic Water into this region (Chiswell et al. 2015; de Burgh-Day et al. 2019). The most notable overall decreases in correlation for the three subregions are from the Tasman index in December and SSTA in January (Fig. 7b) to the index in December and the SSTA in February (Fig. 7f). A similar decrease in correlation is apparent from the Tasman index in January with SSTA in January (Fig. 7a) to the index in January and the SSTA in February (Fig. 7e). This suggests there is a change in SST conditions around NZ between December and February, which is not driven by the Tasman Sea.

Fig. 7.  Correlations between the observed Tasman Sea SSTA index and all observed SSTA values in the region for January, February and March 1990–2012 at lags of 0–2 months. The green lines denote the Tasman Sea index area as well the index regions for Hauraki Gulf, Western Cook Strait and Foveaux Strait.

Of the three key aquaculture regions assessed here, forecasts were most skilful for Hauraki Gulf in the north, and the least skilful for Foveaux Strait in the south. Model ensemble mean forecasts are more skilful than persistence forecasts for most start dates and lead times. ACCESS-S1 also shows promising performance for probabilistic predictions. Historically February is the warmest month in terms of SSTs around NZ (see, for example, fig. 3B of Chiswell and O’Callaghan 2021), and this is therefore when a MWH is most likely to raise SSTs above critical thresholds for aquaculture industries. Skill in summer is lower than winter for all three regions, which is consistent with the findings of de Burgh-Day et al. (2019). They attributed this to increased SST variability in summer compared to winter, which in turn may be related to the shallower mixed layer depth associated with warmer summer SSTs facilitating greater responsiveness to atmospheric variability. Despite this, model skill still exceeds that of persistence (Fig. 2) for the majority of summer month forecasts.

The Hauraki Gulf has statistically significant correlations for 75% of the start date-target date combinations shown in Fig. 2a. Considering only the summer months (December–March), 79% of the start date-target date combinations have skill over a persistence forecast (Fig. 2a). The ability of the model to discriminate between events and non-events is better than climatology for SSTA >90th percentile at all lead times considered. The area under the ROC for Hauraki Gulf is comparable to that of Western Cook Strait, suggesting that the ability of the model to discriminate between events and non-events is similar in the two regions. Despite this, the Hauraki Gulf region has the best reliability of the three regions considered with only the forecast probability of 100% at a lead time of 2 months falling outside the skilful region. The Hauraki Gulf region is fed primarily by subtropical waters from the Tasman Front, which forms the East Auckland Current as it passes around Cape Reinga of the North Island (Chiswell et al. 2015). This makes the higher relative correlation in Hauraki Gulf unsurprising, since ACCESS-S1 has higher skill with these subtropical waters to the west of NZ than any other nearby body of water (de Burgh-Day et al. 2019).

Of the three regions considered, Western Cook Strait has the highest overall ensemble mean correlation, with 82% of the start date/lead time combinations being statistically significant (Fig. 2b). Of the summer months (December–March), 88% start date/target date combinations have skill over persistence. All points fall within the skilful region of the reliability diagram for Cook Strait, except for forecast probabilities of 40 and 100% at a lead time of 2 months (Fig. 4b). Cook Strait is usually fed by the Tasman Sea via the d’Urville and Westland Currents (Fig. 1). ACCESS-S1 is relatively skilful in predicting SST in the Tasman Sea (de Burgh-Day et al. 2019), which in turn likely drives some model skill in the western end of the Western Cook Strait via the d’Urville Current. The lowest ensemble mean correlation is for February forecasts issued on 1 November. Stevens (2014) found that the January–March period is when there is the most north–south variability in flow through Cook Strait, which makes prediction more challenging. ACCESS-S1 has lower skill for the southeastern end of Cook Strait (de Burgh-Day et al. 2019), so change in the strength of flow through Cook Strait could reduce the skill for the western end. Chiswell et al. (2017) have shown that changes in SST in the Western Cook Strait are related to upwelling around Farewell Spit, upstream of the strait, which in turn provides nutrients to mussel farms in Pelorus Sound (Zeldis et al. 2013). Skilful weekly and seasonal forecasts of SST, as well as extreme heat events, in the western half of the Cook Strait could therefore be a beneficial source of guidance for mussel farms and other aquaculture industries in Pelorus Sound and nearby regions.

Foveaux Strait has the lowest overall skill of the three regions. Of the start date/target date combinations, 77% were above the statistical significance threshold (Fig. 2c), and the area under the ROC curve was the lowest of the three regions (Table 2). However, all summer (December–March) start date/target date combinations have skill over persistence forecasts for this region. The model is able to discriminate between events and non-events better than climatology for all lead times considered (Fig. 3), and all forecast probabilities at a lead time of 0 months fall within the skilful region of the reliability diagram (Fig. 4). Foveaux Strait had the highest 0-month lead time BSS referenced to persistence and the lowest referenced to climatology, which is reflective of the seasonally variable nature of SSTA in this region. Foveaux Strait is fed from the west by the Subtropical Water of the Fiordland Current (Chandler et al. 2021), which introduces some influence on skill from Tasman Sea processes (Forcén-Vázquez et al. 2021). ACCESS-S1 has relatively low SST skill to the south of NZ due to a number of influences, including the increased complexity of the currents and the impact of the Southern Ocean warm bias (de Burgh-Day et al. 2019). While water from the Southern Ocean largely follows the boundary of the Snares Shelf where Stewart Island Rakiura and Foveaux Strait are situated, there are variations in supply from the west (Chandler et al. 2021), which potentially explains the poor skill of a persistence forecast.

Regional indices for each of the three regions showed good skill, with forecast SSTA indices matching observed SSTA indices well at 0-month lead times (Fig. 5, Table 2), despite some clear forecast ‘misses’ that emerge at longer lead times (e.g. the forecast failing to capture the two cold anomalies in 2004/2005; Fig. 5). Model skill values for these regional SSTA indices (Table 2), combined with the high correlation between the grid cells in each region and their SSTA index (Fig. 1b–d), indicate that the regional SSTA indices are both representative of the region and have utility for predicting extreme marine heat events.

The data used in this analysis are monthly means, and therefore, any shorter duration extreme heat events will be smoothed out. We apply the general premise of SST values exceeding the 90th percentile to define an extreme heat event, but we do not strictly adhere to the current definition of a MHW as described by Hobday et al. (2016a). This work represents a preliminary analysis of extreme heat event prediction skill in ACCESS-S1 and provides an indication of potential prediction skill for MHWs on the sub-monthly timescale, using the definition of Hobday et al. (2016a). A more detailed assessment of MHW prediction using daily forecasts is planned for the future.

As noted in the Introduction, the regions defined for this analysis are intended to strike a compromise between the limitations of the resolution of the ocean component of ACCESS-S1 (~25 km) and being small enough in extent to provide useful indicators of the SSTA forecast skill of ACCESS-S1 for local aquaculture industries. The skill of ACCESS-S1 forecasts may be sensitive to the extent of the regions chosen, however; so an analysis of the sensitivity of forecast skill to region size and location would be recommended for future work.

While wintertime marine heat events can still affect marine industries, e.g. winter spawning and stocking periods, summertime MHWs are more likely to have more detrimental impacts. Our assessment of probabilistic skill for monthly SSTA >90th percentile skill here used all seasons, due a limited hindcast set, and this overall skill is likely to be higher than that for summer months only, due to the contribution of more skilfully predicted winter months (Fig. 2). To better understand the ability of the model to make useful predictions of summertime extreme heat events in an operational setting, we assessed real-time model forecasts of the 2018/2019 Tasman Sea MHW. The spatial extent of the areas observed to have SSTA >90th percentile was well reproduced at a lead time of 0 months throughout the 2018/2019 event (Fig. 6). However, at longer lead times, there was a marked decline in spatial extent where SSTA >90th percentile. This is particularly true for forecasts from 1 January for February and March. This raises the question of what was driving the 2018/2019 MHW and heat transport from the Tasman Sea to coastal regions around NZ, and whether this can go some way to explaining the decline in model accuracy beyond a 0-month lead time .

Instantaneous correlations between SSTAs in the eastern Tasman Sea and the three key aquaculture regions around NZ are high, but generally drop off rapidly at lags of 1 and 2 months (Fig. 7, middle and right-hand columns). The exception to this at longer leads is Hauraki Gulf, which may be due to the influence of the East Auckland Current, which facilitates the transport of water from the Tasman sea into the Hauraki Gulf region. Holbrook et al. (2020) discuss the physical mechanisms that can drive MHWs, and note that coupled air–sea interactions and atmospheric preconditioning have been the cause of many extratropical MHWs, including the 2017/2018 Tasman Sea MHW. This is consistent with the findings of Elzahaby et al. (2021) that air–sea heatflux-driven MHWs are shallower, and occur predominantly in summer. Bowen et al. (2017) investigated the causes of interannual variability of SST in the Southwest Pacific, and found no evidence for a single mechanism dominating the heat balance. Furthermore, they found that the horizontal advection of heat by boundary currents was unlikely to explain the interannual variability of temperature around NZ. They did however note that air–sea heatflux plays a role in the interannual SST variability in the region. The large drop off in observed SSTA correlations between the Tasman Sea index and the three key regions after the first month suggest that at least some part of the teleconnection between these regions is also atmospherically driven (and thus operating on shorter timescales). This is consistent with the rapid decrease in the spatial extent of areas with a high predicted probability of SSTA >90th percentile at lead times of 1 and 2 months (Fig. 6), as atmospherically driven phenomena tend to have little predictability beyond 1 month (Cowan et al. 2019). In cases where MHWs are ocean-driven, the eastward drift of Tasman waters is likely to have a stronger influence on SSTAs around NZ and is likely to have greater predictability on timescales greater than 1 month (Behrens et al. 2020).

De Burgh-Day et al. (2019) showed that ACCESS-S1 is relatively skilful in predicting SSTAs in the Tasman to the west of NZ, most likely due to the relatively dynamically simple and slow-moving nature of this body of water. This, coupled with skilful forecasts of the Tasman Sea monthly mean index (0.86, 0.67, 0.59 at lead times of 0, 1 and 2 months, respectively), indicates that ocean-driven extreme SSTAs around NZ can be expected to have better predictability than atmospherically driven extreme SSTAs.

Despite the short lead time with which ACCESS-S1 was able to predict the full extent and severity of the 2018/2019 Tasman Sea MHW and its impacts on NZ waters, it did provide indications of increased SSTAs at longer lead times, with higher than climatological odds (i.e. >10% chance) of SSTA >90th percentile in most inshore regions at lead times of 1 and 2 months (Fig. 6). Additionally, the reasonable probabilistic skill of ACCESS-S1 for Hauraki Gulf, Western Cook Strait and Foveaux Strait (see Table 2) gives confidence that future instances of extreme monthly mean SSTAs (>90th percentile) are predictable in those locations at greater than 1-month lead time.

The aquaculture industry in NZ has plans to expand, cognisant that, as it does so, it will be in a changing climate (NZGAS 2019). Based on the current definition by Hobday et al. (2016a), MHWs can occur any time of year, and so impacts will be dependent on industry activities and vulnerabilities to increased ocean temperatures at a given time of year. February is the hottest month for the waters around NZ, and as summertime MHWs generally result in the most extreme absolute temperatures, it is very useful to have skilful forecasts for this time of year. Advanced warning of large marine heat events has significant practical application for aquaculture (e.g. Spillman and Hobday 2014; Hobday et al. 2016b). A high probability of SST in the coming months falling in the warmest 10% of hindcast years suggests increased likelihood of a MHW. Similarly, a high chance of extreme temperatures indicates an upcoming period in which aquaculture industries could experience stock losses. The improved information provided by probabilistic forecasts over a deterministic forecast highlights the strengths of an ensemble prediction system, which aims to sample the distribution of possible outcomes. The benefits to industry of such a system lie in being able to predict shifts in this distribution to higher or lower odds of a given outcome. Importantly, this requires industry managers to develop mitigation strategies that account for uncertainty (Tommasi et al. 2017).

While ACCESS-S1 monthly SST forecasts have been shown to be skilful around NZ and in the Tasman Sea, there is room for future improvement. ACCESS-S1 has a Southern Ocean warm bias, which reduces the mixed layer depth and increases the model SST variability around NZ (de Burgh-Day et al. 2019), which then degrades forecast skill. Another limitation of the ACCESS-S1 system is the relatively short hindcast period of 23 years. This allows fewer samples of more slowly varying modes of climate variability, such as ENSO, which could influence SSTs around NZ, and limits possible skill analyses that can be done due to sample sizes. These shortcomings will be addressed at least in part by the Bureau’s recently-upgraded seasonal prediction system, ACCESS-S2, which has a longer hindcast period (1981–2018) and a weakly coupled ocean–atmosphere initialisation scheme. An extension of this analysis to ACCESS-S2 is left for future work.

This study represents an initial investigation into the potential to produce useful and skilful advance warning of extreme monthly SSTs in the Tasman and around NZ. Our ability to skilfully forecast MHWs that impact NZ coastal regions will depend on our ability to accurately represent processes at several spatial and temporal scales. This initial work has highlighted the complex nature of the teleconnections between the Tasman Sea and NZ. More work is required to further investigate the drivers of prediction skill for MHWs in these regions, including understanding the processes which drive deep and shallow water heat extremes in the region, and how these extremes are transported from the Tasman to NZ coastal waters. An extension of this work in future will be to develop MHW forecast products using daily ACCESS-S forecasts and a formal MHW definition (e.g. Hobday et al. 2016a).

Extreme marine heat events are projected to increase in frequency and severity under climate change, and the Tasman Sea has been identified as a climate hotspot (Hobday and Pecl 2014). Oliver et al. (2017) analysed the 2015/2016 Tasman Sea MHW and found that events of this severity were 6.5 times more likely due to the influence of anthropogenic climate change. The Tasman Sea is projected to continue warming due to an increasing anthropogenic influence (Oliver et al. 2014), further increasing the likelihood of MHWs. Skilful forecasts of ocean heat extremes in regional areas of importance to aquaculture will be a valuable tool for marine operators in this region to mitigate losses due to MHWs, especially in a warming climate.

The ACCESS-S1 Hindcast data used in this study is currently publicly available for research purposes via an OPeNDAP archive, or on the National Computational infrastructure . For further information on access see http://poama.bom.gov.au/general/hindcast_data.html. ACCESS-S1 forecast data for the 2018/2019 Tasman Sea case study period will be shared upon reasonable request to the corresponding author. Reynolds Optimum Interpolation Sea Surface Temperature V2 and V2.1 data were obtained from the US National Oceanic and Atmospheric Administration (NOAA; https://www.ncei.noaa.gov/products/optimum-interpolation-sst) by permission/licence. Data will be shared upon reasonable request to the corresponding author with permission from NOAA, or can be obtained directly.

The authors declare no conflicts of interest.

This work was funded by the National Institute of Water and Atmospheric Research Strategic Science Investment Fund Climate, Aquaculture and Coasts and Oceans Centres.