1.

## Introduction

Currently ocean forecasts with a state-of-the-art model provide short-term forecasts (e.g. daily or shorter time scale) and longer-term forecasts (e.g. seasonal, multi-year to decadal predictions) (e.g. Martin, 2011; MacLachlan et al., 2014; Meehl et al., 2014; Hernandez et al., 2015; Hudson et al., 2017; Schiller et al., 2020). Global ocean prediction became feasible from 1997 due to the Global Ocean Data Assimilation Experiment (GODAE) (Bell et al., 2009). Ocean data assimilation plays an important role for all model-based ocean forecasts (Balmaseda et al., 2013; Fujii et al., 2019; Storto et al., 2019).

A model-based dynamical ocean prediction system can provide forecasts for ocean temperature, salinity, sea surface height and currents. Based on these results, the key ocean phenomena (e.g. meandering and fronts, eddies, surface mixed layer, thermocline depth, equatorial and coastally trapped waves, upwelling of cold water, Rossby waves, storm surge, etc.) can be derived. Therefore, ocean forecasts can provide information to services to locate and map historical shipwrecks and hazardous areas; support search and rescue efforts; assist with ship navigation, structural design, operational safety, oil spill management for oil and gas industry and marine environmental protection; manage fisheries resources and so on.

Short-term ocean prediction has become more and more important owing to more frequent oceanic activities and social-economics developments. The short-term ocean prediction systems are generally based on a high-resolution eddy-resolving or eddy-permitting ocean models. The models can be of a regional or global domain, a stand-alone ocean model or coupled weather/climate model (Usui et al., 2006; Wilkin and Hunter, 2013; Blockley et al., 2014). The forecasting system with a stand-alone ocean model ignores two-way atmosphere-ocean interactions. It needs to be provided with external atmospheric forcing to drive the model and external data for discharge of river runoff during the forecast period. In principle, a coupled ocean-atmosphere model-based forecasting system should be better than the stand-alone ocean model, since it can integrate forward after initialization and river runoff can be produced in the model itself. However, it generally requires more computer resources. For coastal forecast systems, the model domain can be regional or global. For the global domain model, it generally has a nested or multi-nested grid in which the model resolution becomes finer for coastal domains of interest (Brassington et al., 2012). The open boundary conditions in the regional model need to be provided using an interpolation of results of an external and larger domain model possibly with a coarser grid configuration, or an oceanic reanalysis (Wilkin and Hunter, 2013).

A complete list of the operational ocean forecasting system operated worldwide can be found in the web page of GODAE Ocean Predict (https://www.godae-oceanview.org/science/ocean-forecasting-systems/). This website also provides the information on system descriptions, ocean model, assimilation characteristics, routine set-up and national reports. At present stage, most ocean forecast systems only provide daily forecasts with 7 forecast days or so. This study is the first to investigate the predictability of long-range daily forecasts up to 28 days in the Indo-Pacific oceans based on the Australian Community Climate Earth-System Simulator (ACCESS-S1). ACCESS-S1 is the current version of the Australian Bureau of Meteorology’s seasonal forecast system and can provide real-time predictions from daily to multi-week to seasonal timescales (Hudson et al., 2017). It has much higher resolution models than previous one (Wang et al., 2011). The horizontal resolution in previous ocean model is only 2° and it increases to 1/4° in ACCESS-S1. Even though the current ocean model is not fully eddy resolving, it is eddy permitting and does simulate the larger eddies e.g. eddies shedding off the East Australia Current. The focus of this paper is on the ability of eddy-permitting global models (using the ACCESS-S1 model) to predict daily ocean variability beyond the typical 7 or so day forecasts produced by the current high resolution ocean forecast systems.

This study is organised as follows. Section 2 will introduce the model, initialization scheme and ensemble generation method. The model results such as sea surface temperature (SST), sea surface height (SSH) and eddy predictions in the Tasman Sea will be presented in Section 3. The last section will provide a discussion and conclusions.

2.

## ACCESS-S1 system

2.1.

### Model descriptions

The ACCESS-S1 forecast system uses the Met Office Global Coupled model 2.0 (GC2) (Williams et al., 2015). GC2 is comprised of components such as Global Atmosphere 6.0 (GA6.0), Global Land 6.0 (GL6.0), Global Ocean 5.0 (GO5.0) and Global Sea Ice 6.0 (GSI6.0). The atmosphere and land surface horizontal resolution is 0.833°×0.556° (N216) and GA6.0 has 85 vertical levels. GA6.0 and GL6.0 are fully documented by Walters et al. (2017), and GSI6.0 by Rae et al. (2015).

The ocean model GO5.0 is based on version 3.4 (v3.4) of NEMO (Nucleus for European Models of the Ocean) (Madec and the NEMO team, 2016), and is described by Megann et al. (2014). The horizontal ocean grid, known as ORCA025, has 1/4° resolution (1442 × 1021 grid points) at global scale decreasing polewards (an isotropic Mercator grid in the Southern Hemisphere, matched to a quasi-isotropic bipolar grid in the Northern Hemisphere with poles at 107°W and 73°E). This has a 28 km grid spacing near the equator reducing to 6 km grid spacing near the Arctic. The model has 75 vertical levels where the level thickness is a double tanh function of depth such that the level spacing increases from 1 m near the surface to 200 m at a depth of 6000 m. The concentration of levels near the sea surface aims to resolve shallow mixed layers and to give the potential to capture diurnal variability. The ocean model has 35 levels over the top 300 m.

2.2.

### Initialization and ensemble generations

Reanalyses from the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA-Interim project are used to initialize atmosphere and land surface temperatures in the hindcast members (MacLachlan et al. 2014). For the land surface moisture climatological values are used. The oceanic data assimilation uses NEMOVAR, which is an incremental first guess at an appropriate time three-dimensional variational (3DVAR) data assimilation scheme. NEMOVAR assimilates observations of sea-surface temperature (SST), sea-surface height (SSH), in situ temperature and salinity profiles and sea ice concentration (Martin et al., 2007; Donlon et al., 2012; Waters et al., 2015). Note that the oceanic data assimilation here is implemented in the ocean model forced with surface fields from the Met Office numerical weather prediction (NWP) system at a three-hourly frequency. The fluxes are calculated online in the NEMO model using the Coordinated Ocean Research Experiments (CORE) bulk formula formulation (Large and Yeager, 2009) and interpolated in time to each model time step. Sea ice concentration data from Special Sensor Microwave Imager (SSM/I) satellites provided by EUMETSAT OSISAF (http://osisaf.met.no) are also assimilated into the sea ice component.

Two types of perturbation are implemented in the atmospheric model. One uses a Stochastic Kinetic Energy Backscatter (SKEB2) scheme that represents model uncertainty (Bowler et al., 2009). SKEB2 is designed to represent unresolved processes and provide grid-scale perturbations during the model integration. Each ensemble member evolves differently due to these grid-scale perturbations. The other perturbation represents model’s initial condition uncertainty with scaled, randomly chosen 7-day differences of reanalysed atmospheric states from the period 1981–2010. The differences for winds, temperature, specific humidity and surface pressure are multiplied by a scaling factor so that the perturbations have magnitude equal to analysis uncertainty. The scaled perturbations are then added and subtracted from the unperturbed initial atmospheric state each day so that pairs of perturbed members are created (Hudson et al., 2017). The ocean model, sea-ice model and land surface model are not perturbed.

2.3.

### Observed data for model's skill assessment

Hindcast results are verified with observed SST analyses from the NOAA optimum interpolation 1/4° daily SST analysis (OISST), which are available on http://www.esrl.noaa.gov/psd/data/gridded/data.noaa.oisst.v2.html (Reynolds et al., 2007). In addition, we also compare to observed 1/4° daily SSH anomaly (SSHA) altimeter data, which are downloaded from the website: http://marine.copernicus.eu/. They are gridded data from multi altimeter satellites computed with respect to a twenty-year mean reference period (1993–2012). This product is computed with an optimal and centered computation time window. All the missions are homogenized with respect to a reference mission, which is currently Jason-2 (Pujol et al., 2016).

3.

## Results

The hindcasts of ACCESS-S1 initialized on the 1st of February, May, August and November for the period 1990–2012 are investigated here. There are 11 members at each start date and the ensemble mean results are generally used here except for the investigation of eddy predictions where we use individual members in the Tasman Sea.

3.1.

### Model SST bias

The model climatology is the daily mean for model hindcast results over the period 1990–2012. This study investigates the long range daily forecast up to 28 days. Each forecast day has a corresponding climatology, e.g. there are a total of 28 climatological means initialised from 1st February. Figure 1 displays the model SST bias. Left panels represent the model biases at the 7th, 14th, 21st and 28th forecast days starting on the 1st of all months. Right panels are the model biases at the same 21st forecast day but for different months. The SST bias in 7th and 28th seems similar probably due to the high persistence of SST. However, the absolute errors averaged over the Indo-Pacific Oceans (30°E-70°W, 50°S-60°N) continuously grow from 0.17 °C at lead day 1 to 0.29 °C at day 28. In the first 28 days of the forecast, SST biases are in the range of ±0.3 °C in most regions, except in some coastal areas and some regions with high variability such as the eastern equatorial Pacific and Kuroshio. The Pacific cold tongue bias, which is common in many climate models, reaches up to about −0.8 °C on the 28th lead day and it will further increase beyond this lead time and finally is up to around −2 °C after a few-month forecast time (Zhou et al., 2015). The Pacific Ocean has a small warm bias in most areas except the northeast Pacific between 40°N-60°N and eastern equatorial Pacific. Conversely the North Indian Ocean and western Australian coastal regions tend to have a cold bias which increases with lead time. The spatial patterns of SST bias are quite different in different months. The model has a warm bias in most of the South Pacific in February (Fig. 1e), but the warm bias is in the North Pacific in May and August (Fig. 1f,g). This may be related to too much heat remaining in the upper ocean due to ocean stratification forming in summer. In August, the model has a relatively large cold bias in the eastern Indian Ocean. Our further analysis suggests that the cold bias is mainly due to too strong surface southeasterly winds in the Indian monsoon season which results in too strong oceanic upwelling in this region.

Fig. 1.

SST (°C) bias at the forecast lead time (a) 7th day, (b) 14th day, (c) 21st day and (d) 28th day. (e) The SST bias at the 21st forecast day initialised from 1st February. (f), (g) and (h) are the same as (e) except for May, August and November, respectively.

3.2.

### SST anomaly (SSTA) prediction skill

The model SSTA is computed by removing the model climatology which is defined in 3.1. No seasonal signals are in the model SSTA. Figure 2 shows model SSTA correlation skill (left panels), persistence skill (middle panels) calculated with observed daily data and their differences (right panels) at a 7-day interval over the first 28-day lead. The model's correlation skill is the correlation between model predicted results and the corresponding observed ones at each lead time over the whole hindcast period. The persistence skill is commonly used as a reference for the model's prediction skill. The persistence correlation is the autocorrelation from the initial conditions calculated from the OISST data. Persistence SSTA forecasts have very high skill in the equatorial Pacific (>0.9). The spatial patterns of model correlation skill are similar to those of the persistence skill at the same forecast lead time. A high/low predicted SSTA correlation skill in a specific area coincides with a high/low persistence skill. Model skill beats persistence skill in most regions at two weeks lead time and beyond. The pattern of SSTA persistence skill is similar to that of a regression onto the monthly Niño 3.4 SSTA index of global SSTA (Zhou et al. 2015). The high persistence skill in the low and middle latitude Pacific occur in the regions affected by the El Niño and Southern Oscillation (ENSO). It suggests that the SSTA daily prediction skill is likely linked to the interannual time scale ENSO events. The North Indian Ocean has a lower persistence skill than the South Indian Ocean. This could be because the South Indian Ocean is largely affected by the western Pacific through the Indonesian Throughflow.

Fig. 2.

Left panels: model’s predicted SSTA correlation at forecast time (a) 7th day, (b) 14th day, (c) 21st day and (d) 28th day starting from the 1st of February, May, August and November over the period 1990–2012. Middle panels (e, f, g and h) are the same as the left panels but for persistence skill. Right panels (i, j, k and l) are the differences between correlation of the model and persistence.

Figure 3 shows the model SSTA correlation skill at the 14th forecast day starting from 1st February and the differences of correlation between February and the other months such as May, August and November. The averaged correlation coefficients over the Indo-Pacific oceans in all months are similar to each other. They are about 0.7 at forecast day 1 and decrease to about 0.5 at day 28. However, the spatial structures of correlation skill in different months are different. The model forecasts in February have higher correlation skill in the most North Pacific than those in August and May, whereas lower prediction skill in the seas of north and east of Australia and east of South America compared to the other months.

Fig. 3.

(a) Model correlation skill for predicted SSTA at 14th forecast day initialised from 1st February. (b), (c) and (d) are the difference between February and May, August and November, respectively.

Figure 4 exhibits the predicted SSTA root mean square error (RMSE, left panels) and RMSE normalised using the observed SSTA standard deviation (right panels). Large RMSE values occur in the equatorial Pacific, Kuroshio and ACC regions where there is strong variability. The magnitude and spatial pattern of RMSE at 7 days lead seem to have similar pattern to that at 28 days lead. It suggests that the model's forecast errors increase rapidly in the first week or two.

Fig. 4.

Left panels: model’s predicted SSTA (°C) root-mean-squared errors (RMSEs) at the forecast lead time (a) 7th day, (b) 14th day, (c) 21st day and (d) 28th day over the hindcast period 1990–2012. Right panels are the same as the left panel but for normalized RMSE using observed SSTA standard deviations.

3.3.

### Prediction skill of SSTA indices

Eight SSTA indices, which represent different domains in the Indo-Pacific, are defined in Fig. 5. The SSTA averaged in the area of 180°E-70°W and 20S°-20°N stands for eastern tropical Pacific (EP); 120°E-160°E and 20°S-20°N for western tropical Pacific (WP); 120°E-120°W and 20°N-45°N for the North Pacific (NP); 120°E-70°W and 20°S-45°S for the South Pacific (SP); 30°E-80°E and 0–25°N for the North Indian Ocean (NI); 30°E-110°E and 0–40°S for the South Indian Ocean (SI); 110°E-120°E and 15°S-30°S for the Leeuwin Current region (LC) and 150°E-160°E and 20°S-40°S for East Australian Current region (EAC).

Fig. 5.

Different domains used for calculating SSTA indices.

Figure 6 displays the prediction correlation skill of the above eight SSTA indices. The model prediction skill (red solid line) is lower than persistence (black dash line) in the first few days, but higher at longer leads out to 28 days. This initial poorer correlation is related to the model initialization method. Firstly, the ocean model initialization from NEMOVAR reanalysis is implemented with the ocean model alone forced with external atmospheric forcing rather than the coupled GC2 model (Waters et al., 2015). This leads to a possible initial shock during the forecast period. Secondly, the OISST reanalysis is regarded as ‘truth’ here and used to measure the model's prediction skill. Therefore it is likely that the NEMOVAR reanalysis will disagree with OISST, since the OISST reanalysis are obtained with objective optimal interpolation of the observed satellite SST (Reynolds et al., 2007), but the SST in NEMOVAR is derived from assimilating the observed satellite and in-situ SST data into the model (Waters et al., 2015). Therefore, the SST data in OISST and NEMOVAR reanalysis are, in fact, obtained from the almost same observed SST data but with a different method. It should not mean that SST in OISST is better and more accurate than that in NEMOVAR reanalysis. We checked the correlation of NEMOVAR reanalysis and Reynolds observed daily SSTA on the 1 February over the period 1990–2012. It is below 0.5 in some areas of the western Pacific, northern Indian Ocean and ACC (not shown). Therefore, it is not surprising that the model has a lower correlation skill than the persistence correlation in the first few days.

Fig. 6.

Model predicted correlations (red solid) and persistence correlations (black dash) for SSTA indices in (a) east Pacific, (b) west Pacific, (c) north Pacific, (d) south Pacific, (e) north Indian Ocean, (f) south Indian Ocean, (g) Leeuwin Current and (h) East Australian Current.

In general, the model skill exceeds persistence skill beyond one-week forecast time. For EP, EAC and SP indices, the model’s correlation skill is just a little higher than the persistence. Note that the persistence skill for the EP index is very high and even above 0.9 at the lead 28th day probably due to ENSO being active in this region. For SI and LC indices, the model skill is significantly higher than the persistence (Fig. 6). For the other three indices, the model skill is relatively higher than the persistence. The LC and EAC indices represent the Australian western and eastern coastal regions, respectively. The model seems to have a better performance on the Australian western coast than the eastern coast.

The model's RMSE (red solid line), RMSE of persistence forecasts (black dashed line) and average ensemble spread (green dashed line) are displayed in Fig. 7. During the first few days, the model RMSE is larger than the persistence. It is consistent with a lower correlation skill than the persistence in the same lead time. We had discussions for this issue before. Regarding the EP index, the model's RMSE is almost the same as persistence beyond a few days. For the indices EP, NP, SP, the model only has a smaller improvement in RMSE skill, whereas the model has more improvement in the Indian Ocean. In particular, the model has a significant improvement for the LC and EAC indices, which represent the Australian coastal regions. In general, the model has less room to improve the RMSE skill in basin scale than coastal regions, since the RMSE in basin scale is much smaller than that in coastal areas.

Fig. 7.

Same as Fig. 6 but for RMSE and ensemble spread (green) (unit: °C).

If the model is given perfect perturbations to generate ensemble members, the model's RMSE and average ensemble spread should be approximately equal based on the formula (Fortin et al., 2014):

$\mathit{RMSE}\approx \sqrt{\left(\frac{R+1}{R}\right)\frac{1}{T}\sum _{t=1}^{T}{s}_{t}^{2}}=\sqrt{\left(\frac{R+1}{R}\right)}{\left(\overline{{s}_{t}^{2}}\right)}^{1/2}$
where $R$ is ensemble member and ${\mathrm{}s}_{t}^{2}\mathrm{}$ is ensemble variance at time t. RHS of the above equation stands for the average ensemble spread. Note that it is not equal to the average ensemble standard deviation. For a sufficiently large ensemble size, the correction factor $\sqrt{\left(\frac{R+1}{R}\right)}$ that depends on $R$ vanishes. Generally, the model average ensemble spread is smaller than the corresponding model's RMSE, particular in the Indian Ocean It indicates that the model's perturbation is too small to sufficiently account for the model's uncertainty. This problem is likely due to the fact that the model's ocean, sea ice and land surface models are not perturbed and only the atmospheric model is perturbed here. It is interesting for the EAC index that the ensemble spread increases quickly and is almost comparable to RMSE after the 12-day forecast time.

3.4.

### Sea surface height anomaly (SSHA) prediction skill

Variations in SSHA reflect the variations of subsurface heat content via the thermosteric component. For a good prediction for SSHA it is useful to know the status of both sea surface and subsurface in advance. The SSHA model correlation and RMSE skill is displayed in Figs. 8 and 9, respectively, at lead times 7, 14, 21 and 28 days. SSHA has a relatively low persistence skill in the central and eastern equatorial Pacific, which may be associated with tropical instability waves, but a high persistence skill in the western tropical Pacific possibly relating to the sea water accumulation in this region forced by the easterly trade winds. In contrast, SSTA has a high persistence in the central and eastern persistence skills but a relatively low skill in the western Pacific.

Fig. 8.

Left panels: model’s predicted SSHA correlation at forecast time (a) 7th day, (b) 14th day, (c) 21st day and (d) 28th day starting from the 1st of February, May, August and November over the period 1990–2012. Middle panels (e, f, g and h) are the same as the left panels but for persistence skill. Right panels (i, j, k and l) are the differences between correlation of the model and persistence.

Fig. 9.

Left panels: model’s predicted SSHA (°C) root-mean-squared errors (RMSEs) at the forecast lead time (a) 7th day, (b) 14th day, (c) 21st day and (d) 28th day over the hindcast period 1990–2012. Right panels are the same as the left panel but for normalized RMSE using observed SSHA standard deviations.

Persistence correlation skill is higher than or comparable to the model prediction skill for the 7 days lead time in most regions except the Nino3.4 area (Fig. 8e). It indicates that the model needs to improve the assimilation of satellite altimeter data. However, with as lead time increases, the model shows its superiority and beats persistence in the most parts of the tropical Indo-Pacific. Note that the model prediction skill is always lower than the persistence almost everywhere poleward of 30oN or 35oS and in the tropical western Pacific. The model's ability to beat persistence is lower for SSHA than SSTA, likely because of the high persistence of SSHA compared to SSTA.

The model SSHA forecasts have large RMSE in regions such as Kuroshio, ACC and EAC that have a large variability and are dominated by eddies. The spatial pattern is similar to that of SSTA RMSE. For the normalized RMSE with the observed SSHA standard deviation, the ACC and Kuroshio regions the errors still remain high.

3.5.

### Forecasting mesoscale eddies in the Tasman Sea

ACCESS-S1 is able to predict not only basin-scale oceans, but also coastal regions. The prediction skill for two SSTA indices LC and EAC, defined in Section 3.3 to represent Australian coastal regions, have been examined. Here we will further investigate the predictability for mesoscale eddies in the Tasman Sea with ACCESS-S1. Woodham et al. (2015) explored eddy forecasts in this area with a high-resolution regional ocean model BLUElink. The model, with a horizontal resolution of 0.1° in Australian coastal regions, can only provide 7-day daily ocean forecasts. Their results showed that the model lost ability to forecast some eddies even after 72 hours.

The EAC is a strong western boundary current in the Australian coastal region. It flows southward from Coral Sea around 15°S and tends to be intensified till its separation point around 32°S. The EAC produces relatively large mesoscale eddies which have significant impacts on the physical, chemical, and biological properties of the Tasman Sea. The Tasman Sea is situated around the EAC’s separation point and the EAC extension region and this area is eddy-rich and of high variability. It is a marginal sea of the South Pacific Ocean, situated between Australia and New Zealand. The area of 32°S to 37°S and 150°E to 155°E in the Tasman Sea is analysed here to evaluate the predictability of mesoscale eddies. The relatively large size of the eddies that separate from the EAC are represented by the eddy permitting model used here.

Two cases for model's eddy prediction in the Tasman Sea are investigated here. One shows a good prediction starting from 1 November, 2009 (Fig. 10) and the other is a poor prediction issued 1 August 2003 (Fig. 11). Figures 10 and 11 display the time-evolving eddies (based on SSHA variable) in the Tasman Sea with a three-day interval up to the 28th forecast day for the model predictions (left two columns) and the corresponding observations (right two columns). Generally, a cold-core eddy (CCE) upwells the water in the centre of the eddy and a warm-core eddy (WCE) downwells the water. The upwelled water comes from deeper in the ocean and, hence, is cooler and contains nutrients. The downwelled water results in warmer temperature anomalies and less surface nutrients.

Fig. 10.

Left two panels (a) are model’s predicted SSHA in the Tasman Sea for 1st to 28th November 2009 with a 3-day interval, issued 1 November 2009. Right two panels (b) are the corresponding observed SSHA (unit: m).

Fig. 11.

Same as Fig. 10 but for mode forecasts starting 1 August 2003.

Figure 10 shows the forecasts starting on 1st November 2009. A strong CCE and a strong WCE are present in the Tasman Sea. They move southward along the coast. The model can capture the tracks, centre locations and strength of the two eddies up to one month in advance. However, the CCE in the model becomes a little weaker as forecast lead time increases, whereas the observed CCE remains stable till the 28th November. The observations show that a second WCE located to the north of the CCE moves southwards and intensifies (Fig. 10b). However, the model predicted centre of the second WCE is too far east. Furthermore, a second CCE which follows the second WCE moves southward along coastline earlier than the observations.

Figure 11 shows the SSHA ensemble mean forecast starting on 1st August 2003, together with observations (Fig. 11b). The model fails to predict the eddy evolution beyond 7 days forecast lead time. For the 1st forecast day, the model’s eddy types, locations and strength are close to the observation (Fig. 11b), but the CCE strength in the north Tasman Sea is a little underestimated. The CCE in the model become weaker and then almost dissipates as forecast lead increases. A new weak CCE forms in the west of the WCE and the WCE shifts eastward (Fig. 11a). However, the observations show that one strong CCE in the north and one weak CCE in the east of the strong WCE on 1 August 2003. The WCE moves southwards continuously whereas both CCEs seem to stay in place and the south CCE become stronger with time.

Figure 12 shows all 11 forecast members for eddies and sea surface wind stresses issued on 1 August 2003 at the 28th forecast day. The eddies in terms of location and strength look quite similar to each other, suggesting that the model ensemble spread is too small. However, the external atmospheric forcing such as wind stress and total heat flux (not shown) are quite different in each member due to perturbation of both initial conditions and model physics in the atmospheric model. This may indicate that the ocean internal dynamics plays the major role for formation and movement of these eddies. The ocean model physics and ocean initial conditions need to be perturbed in the future to increase the ensemble spread for ocean variables and improve eddy predictability.

Fig. 12.

Model individual SSHA forecasts (unit: m) and sea surface wind stress (unit: N/m2) in the Tasman Sea on 28 August 2003 (28 day lead), issued 1 August 2003.

Figure 13a displays the SSHA predicted/observed spatial correlation skill for each forecast case (ensemble mean) at forecast time 1st day (black), 7th day (red) and 14th (green) day over the period 1993–2012. The correlation for the 1st day is mostly between 0.6 and 0.9. The 7th and 14th day prediction correlations are highly related to the 1st day correlation. If the 1st day prediction correlation is high/low, the prediction skill of 7th and 14th day is generally high/low. This indicates that the ocean model initialization plays a major role in determining forecast skill. However, the prediction skill for some cases drop very quickly with increasing forecast lead time, even if the model states are quite closed to observations at the initial time. It may be caused by the initialisations. The ocean eddies are usually hundreds of meter deep or even deeper. Although the model's eddies are closed to the observed ones at the sea surface due to the data assimilation satellite altimeter data, the model's eddies in the subsurface may not be corrected properly due to less observed data.

Fig. 13.

(a) The spatial correlation between the model prediction and observed SSHA in the Tasman Sea (the domain as per Figs. 10–12) at forecast day 1 (black), 7 (red) and 14 (green). (b) The model averaged spatial correlation (red line) at forecast time from day 1 to 28 over the period 1993–2013. Dots indicate forecast correlations for model start dates. There is a total of 80 cases at the same lead time. (c) Same as (b) but for RMSE (unit: m).

Figure 13b shows the model spatial correlation with observations with forecast lead time for each model start date. The black dot indicates the correlation skill for each start date. The red line stands for the averaged correlation for all start dates at a given lead time. The difference of the highest and lowest correlation skill becomes bigger with increasing lead time. The model averaged correlation skill is around 0.76 from 1 day lead and drops to about 0.5 after a couple of weeks. Figure 13c is the same as Fig. 13b but for RMSE. The RMSE increases from 0.12 m at the 1st day to about 0.2 m at the 28th lead day. This is consistent with the drop of prediction correlation skill with increasing lead time (Fig. 13b).

4.

## Conclusions and discussions

This study takes advantage of the availability of eddy permitting global seasonal prediction models to explore the potential to produce daily ocean forecasts to forecast lead times beyond those available for global ocean forecast systems. We explore the potential for extended long range daily ocean forecast up to 28 days into the future. This is the first time this has been done with an eddy permitting seasonal forecast coupled model, ACCESS-S1 in this study. Our results indicate that ACCESS-S1 is able to provide accurate and useful daily predictive information beyond two-week lead time. This suggests that: (i) the lead time of around 7-days used in most operational daily ocean forecast systems can be extended by using a coupled model; (ii) the forecast system with a coupled climate model has greater advantage in accuracy and no limit in forecast time than that with ocean model alone which is largely restricted by external NWP forcing during the forecast period. Although the horizontal resolution of the ocean model used here is 1/4° and it is only eddy-permitting rather than eddy-resolving model, the national report in 2019 from the Met Office points out that there is generally similar accuracy in the short-range forecast for the NEMO model with 1/4° and 1/12° horizontal resolutions, in spite of some of the dynamics being better represented in the higher resolution model (https://www.godae-oceanview.org/documents/q/category/govst/system-reports/national-systems-reports-2019/).

The prediction skills for daily SSTA and SSHA in the Indo-Pacific over the period 1990–2012 have been assessed. The daily SST bias up to the 28th forecast day in ACCESS-S1 is within 0.3 °C in most regions except some areas with strong currents or strong temperature variability, e.g. Kuroshio, ACC and eastern equatorial Pacific. The model SSTA prediction skill is higher than persistence skill beyond the first few days of forecast lead time. Eight SSTA indices representing different domains are defined here. The model has much higher prediction skill in the northern Indian Ocean, southern Indian Ocean, western Pacific and Leeuwin Current than persistence skill. However, for the regions in the eastern Pacific, southern Pacific and EAC, the model skill is only slightly higher or comparable to persistence skill. It is noted that the SSTA in eastern tropical Pacific has high persistence skill (correlation >0.9 at the 28th forecast time). It is likely related to interannual ENSO events.

The model’s SSHA prediction skill is comparable or slightly lower than persistence skill in the 7th day of the forecast in the tropical Indo-Pacific. However, it exceeds the persistence from the second week in the most regions of the tropical Indo-Pacific, except in the western Pacific where persistence skill is high till the fourth week of the forecast. Eddy prediction in the Tasman Sea is also explored in detail. The averaged spatial correlation (total 80 cases for a forecast time over the period 1992–2012) between the model prediction and observations is about 0.5 at the 14-day lead. Two example cases are presented: a good prediction and a poor prediction. For the good prediction case, the model can predict three eddies' strength and location a few weeks in advance. However, the model can only predict well in the first week for the poor case. After one week, the model’s predicted eddies for the poor case gradually show different tracks to the observed situation. In general, ACCESS-S1 has demonstrated an ability to predict the large-scale eddies in the Tasman Sea, which is an achievement given that these eddies are very challenging to predict due to their highly variability.

Based on the above evaluation, the ACCESS-S forecasts can be improved further in two ways. One is through better data assimilation and the other is via improving the ensemble generation method. Currently, ocean data assimilation is implemented in a standalone ocean model which is forced with external atmospheric forcing (Waters et al., 2015). The atmospheric initial conditions for the forecast are then simply interpolated onto the model grids from the reanalysis of ECMW ERA-Interim. A better way would be to apply a coupled data assimilation method. It could potentially produce a balanced state for all components and each variable in the climate model and avoid an initial shock during the forecast period. The model ensemble generation also needs to be improved, since the ensemble spread is much smaller than the corresponding RMSE. They should be comparable to each other for a reliable ensemble. Currently, the ocean model is not perturbed and only the atmospheric model is perturbed. All model components for both model physics and initial conditions should be perturbed in order to improve the predicted reliability.

This study is an initial validation of daily prediction skill of sea surface temperature and height out to 28 days lead. Further work is underway including: (i) evaluation of other variables in the upper subsurface, such as ocean currents, mixer layer depth, thermocline depth; (ii) investigation of multi-week (weekly averaged) ocean forecasts which is a bridge between daily and seasonal forecasts; (iii) validation of more hindcast experiments from ACCESS-S1 over the period 1990–2012 starting from four different dates (1st, 9th, 17th and 25th) in all calendar months. It will produce more statistically robust results on predictability due to more case studies.