Long range daily ocean forecasts in the Indo-Pacific oceans with ACCESS-S1

2.1. Model descriptions

The ACCESS-S1 forecast system uses the Met Office Global Coupled model 2.0 (GC2) (Williams et al., Citation2015). 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. (Citation2017), and GSI6.0 by Rae et al. (Citation2015).

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. (Citation2014). 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. Citation2014). 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., Citation2007; Donlon et al., Citation2012; Waters et al., Citation2015). 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, Citation2009) 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., Citation2009). 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., Citation2017). 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., Citation2007). 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., Citation2016).

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., Citation2015). 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.

Display full size

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.

Display full size

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.

Display full size

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.

Display full size

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.

Display full size

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., Citation2015). 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., Citation2007), but the SST in NEMOVAR is derived from assimilating the observed satellite and in-situ SST data into the model (Waters et al., Citation2015). 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.

Display full size

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).

Display full size

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., Citation2014): $$\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.

Display full size

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.

Display full size

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 30^oN or 35^oS 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. (Citation2015) 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).

Display full size

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

Display full size

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/m²) in the Tasman Sea on 28 August 2003 (28 day lead), issued 1 August 2003.

Display full size

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).

Display full size

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).