The Australian wave forecast system with statistical wind correction

2.1. Atmospheric model

The UK Met Offices’ unified model is the basis for the APS weather model. The fourth upgrade (APS4) contains the UK Met Offices’ Parallel Suite 45 code, designed to run with a dynamic ocean model. However, the current implementation uses static sea surface temperature fields. Coupling to a dynamic ocean is planned at a later stage.

Fig. 1 and Table 1 show verification metrics for various weather models relative to marine advanced scatterometer (ASCAT) winds. A couple of things stand out. First, APS4 shows a reduced wind bias relative to APS3, from −0.27 to −0.12 m s⁻¹. Second, the root-mean-square error (RMSE; RMSE=[n−1∑in( m i− o i)2]1/2) in wind speed improved from 1.27 to 1.26 m s⁻¹. Third, there is a better representation of the probability density function (central panel) in APS4. However, the probability of exceedance P(X > x) shows lower skill in APS4, meaning that high wind speeds in the range of 18–35 m s⁻¹ are less accurate. Interestingly, APS3 is slightly better for P(X > x) than APS4. For comparison, Durrant et al. (2013) found a negative bias (B; B=[n−1∑inmi]−[n−1∑inoi]=mean(m)−mean(o)) of −0.52 m s⁻¹ with RMSE of 1.57 m s⁻¹ in the APS1 wind field and applied bias correction. The NWP reanalysis, for example, ERA5 (Hersbach et al. 2020) is known for biases at high wind speed; hence, statistical correction for extreme winds is valid. The bottom-left panel in Fig. 1 shows that the probability of exceedance of ERA5 surface wind speed is under-predicted, starting at ~17 m s⁻¹, also reflected in a lower value for slope and a lower score in histogram fit (I; I=∑inmin(mi,oi)÷∑inoi; Swain and Ballard 2002) (see legend in the central panel). Nonetheless, the ERA5 reanalysis features the best skill regarding RMSE, scatter index (SI; Cardone et al. 1996; s.d.(m – o) ÷ mean(o)) and Pearson’s correlation coefficient (R). It is therefore considered the best product in the world.

Fig. 1.

Satellite verification statistics for marine winds from ACCESS-G3 (top), ACCESS-G4 (centre) and ERA5 (bottom). Statistics were averaged globally from July 2023 to August 2024. Panels show the probability of exceedance (left), probability density function and scatter density (right). Statistics in the legend are the number of observations (N), Pearson’s correlation coefficient (R), root-mean-square error (RMSE), bias (B), scatter index (SI), least squares slope forced through the origin (fit₀) and histogram fit (I).

Fig. 2 shows the spatial distribution of mean wind bias from July 2023 to August 2024. It confirms the reduction in wind speed bias with fewer blueish regions visible in APS4 overall than in APS3. The most noticeable improvements are along the Antarctic ice edge, the equatorial Pacific and the eastern Indian Ocean. The western Indian Ocean remains biased low. Spatial patterns in ERA5 reanalysis bias are like those in ASP4, perhaps except for the Southern Ocean, where more frequent underestimation of surface winds occurs.

Fig. 2.

Spatial mean wind speed bias from (top to bottom) ACCESS-G3, ACCESS-G4 and ERA5. Observations refer to ASCAT satellite data from July 2023 to August 2024.

The negative bias in wind speed originates from a prescribed static wave-boundary layer. The NWP operates as a standalone model. Two-way coupling can reduce biases at high wind speeds and decrease the scatter (Wiese et al. 2020). Janssen and Bidlot (2018) illustrated the importance of two-way coupling and stress, stating that the principal component is the transfer of momentum and energy to the surface. Breivik et al. (2022) showed that reducing the Charnock parameters at high wind speeds can increase wind speed in the wave-boundary layer.

2.2. Wave model

AUSWAVE-G is based on a 12-km multiple-resolution SMC grid, with refinement around sub-grid scales features and shelves at ~6 km and features the observation-based wind input and dissipation terms (Zieger and Greenslade 2021). The legacy AUSWAVE-G APS3 uses the GEBCO 2020 (GEBCO Bathymetric Compilation Group 2020) gridded bathymetry data. Bathymetric depths were evaluated against measured depths derived from wave buoy deployments around the Australian coastline. The APS4 upgrade did not change the wave model grid and underlying bathymetry.

Surface gravity waves grow due to surface winds. Wave equilibrium assumes the dominant terms (i.e. wind input, resonant wave–wave interaction and wave breaking) on the right-hand side in balance and the left-hand side of the governing equation to be zero at first order (Thomson et al. 2013). Equilibrium wave theory shows that significant wave height H_s scales with wind speed U102:Hs∝U102 (Janssen 1998). This means that the error in the wind field is related to the error in the wave field. A detailed skill assessment of the global wind field is essential for the predictability of the wave field. Satellite verification yields an error of 1.26 m s⁻¹ in the surface wind field (global averages), which would result in an error of 0.41 m in significant wave height. The negative bias in the ACCESS wind field (APS1) yields a negative bias in the wave field with a H_s bias of −0.22 m and RMSE of 0.49 m (Durrant et al. 2013). Durrant and co-authors concluded that the bias in H_s is highly sensitive to the bias in U₁₀, overall and in spatial variation.

The global wave model AUSWAVE-G APS3 increased the skill on average by 13% compared to the previous version (Zieger and Greenslade 2021). The improvement comes from increasing the flux to the waves by 5% to compensate for inherent biases in surface wind speed. Setting the FLX4 namelist parameter CDFAC to 1.05 increases the fluxes from the atmosphere to the waves by 5%. Initial trials of the legacy AUSWAVE-G3 model with ACCESS-G4 surface winds were unsatisfactory. They resulted in a large overprediction of the wave field with less skilful bulk scatter statistics (i.e. bias, RMSE, SI and R), as shown in Table 2. Scaling the drag coefficient would change the momentum flux across the entire range of surface winds, leading to an overestimate in the low range. The wave model then integrates these overestimates, producing less skilful wave predictions. As a result, to compensate for biases in the wind field, a bias correction will be adopted.

Verification of concatenated 12-h forecasts of ACCESS-G4 reveals a slight bias in wind speed (see Fig. 1). A close look at the probability of exceedance plot (left panel) indicates that the surface winds are in excellent agreement to ~15 m s⁻¹. An adjustment of only a fraction of the winds is needed. To do so, we apply a threshold-based bias correction to target the underprediction of extreme surface winds with wind magnitudes that exceed the threshold w₁ are scaled with w₂ as defined in Eqn 1.

Table 2 lists configurations with variations in wind speed correction and the resulting altimeter verification for the wave model simulation. Since the surface wind speed skill is remarkably close between the APS3 and APS4 versions of the NWP model, one would expect a similar skill between the AUSWAVE versions. The change in wind speed bias from −0.27 to −0.14 m s⁻¹ will overestimate the sea state everywhere with a mean bias of 0.07 m in significant wave height and 2% increase in SI, as listed in Table 2. The threshold-based approach applied here can alleviate this overestimate. A closer look at the metrics in Table 2 reveals that the bias correction improves the legacy model G3 in SI, R and bias. Fig. 3 shows that the G4 wave model matches the probability density function of the observations and represents extreme wave heights exceptionally well to the point that it can beat ERA5 in skill. One can conclude that open ocean fields of H_s are verified well by altimeter observations.

Fig. 3.

Satellite verification for significant wave height. Statistics are for the AUSWAVE-G3 operational model (top), the configuration for AUSWAVE-G4 (centre) and ERA5 (bottom). Statistics are averaged between 55°N and 55°S from 1 August 2023 to 1 January 2023 and 22 February to 1 August 2024. Panels show the probability of exceedance (left), probability density function and scatter density (right). Statistics in legend are the number of observations (N), Pearson’s correlation coefficient (R), root-mean-square error (RMSE), bias (B), scatter index (SI), least squares slope forced through the origin (fit₀) and histogram fit (I).

Fig. 4 shows the spatial distribution of mean H_s observed by altimeters. ASP3 and APS4 are similar in bias patterns. At high latitudes, the APS4 bias (middle panel) is slightly larger than that of APS3, which may result from the somewhat more energetic surface winds (see Fig. 2). By contrast, these regions are prominent for long-period swells, indicating a lack of understanding of swell dissipation. The most noticeable difference in ERA5 reanalysis is a more prominent negative bias in H_s in the Southern Ocean, North Pacific and North Atlantic. The G4 model features a more negative bias in the west-central equatorial Pacific, while a positive bias persists east of Drake Passage in the South Atlantic.

Fig. 4.

Spatial mean bias in significant wave height. Panels show from top to bottom: AUSWAVE-G3, AUSWAVE-G4 and ERA5. Observations refer to altimeter data with statistics averaged between 55°N and 55°S from 1 August 2023 to 1 January 2023 and 22 February to 1 August 2024.

The third-generation spectral wave model WW3 explicitly contains source terms representing swell decay and negative input. The wave model in the ERA5 reanalysis does not contain an explicit swell decay term. Instead, the level of white-capping dissipation was retuned to mimic the slow attenuation of long-period swell (Hersbach et al. 2020). The negative input in WW3 represents non-breaking dissipation due to turbulence in the atmosphere–wave boundary layer. Although the swell decay rate has been observed in field experiments, one cannot distinguish background interaction on either side of the sea surface, i.e. in the air or water (Young et al. 2013). Pathirana et al. (2023) found that swell waves are attenuated too strongly for swells propagating long distances and recommended a lower value for negative wind input in observation-based wave physics (ST6). Zieger et al. (2015) and later Liu et al. (2019) fully describe the ST6 physics package. An earlier version of the configuration with a negative input of 5% and a constant swell decay coefficient is selected (Zieger et al. 2015), consistent with the conclusion of Pathirana et al. (2023). In addition, the background turbulence swell dissipation rate is constant b₁ = 0.22 × 10⁻³. The rest of the model configuration follows the configuration described in Zieger and Greenslade (2021).

3.1. Forecast verification

Results presented in Section 2 demonstrated the skill of the new model configuration compared to the legacy operational wave model and ERA5 reanalysis for open ocean and deep-water wave conditions from altimeter observations. The wave hindcasts were forced with concatenated 12-h NWP forecasts, including analysis fields and short-range forecasts, to minimise the accumulation of forecast errors. Hindcasts are free simulations with no explicit use of data assimilation. Under the assumption that forecast errors grow with lead time, these short-range forecasts will have the most negligible errors associated with them.

We cycled the model every 12 h to assess the forecast skill and produced 7-day forecasts for March–November 2024 to generate a database for forecast verification. Fig. 5 shows the skill for H_s as a function of forecast lead time and demonstrates that the proposed G4 configuration outperforms or matches the legacy global system (G3) entirely. The RMSE (top panel in Fig. 5) does not change significantly over the first +2 to +3 days and then grows at a higher rate during days 3–7. Interestingly, the forecast error is persistently lower than the error of the legacy model (G3) and accumulates at a lower rate with forecast lead time. A similar pattern is evident when looking at R in significant wave height (second panel from the top in Fig. 5). R is close to the legacy system up to the 3-day forecast range. Beyond 3 days, the new model features less R loss. AUSWAVE-G3 showed an increase in bias with a lead time of up to 5 days (bottom panel in Fig. 5). The G4 model configuration reduced the model bias with lead time.

Fig. 5.

Model skill for significant wave height plotted as a function of forecast time (hours, horizontal axis). Panels show verification metrics (top to bottom): root-mean-square error, Pearson’s correlation (R), scatter index (SI) and bias (B). Line colours indicate different models, as indicated in the legend.

Zieger and Greenslade (2021) showed that the forecast errors in APS3 surface winds grow linearly across the full forecast range. The wave model takes ~3 days in lead time to accumulate errors in new wind-sea from the NWP model and propagate these uncertainties across the sea as swell. Swell waves differ from local wind-sea and originate from remote storms much earlier. The data assimilation in NWP models ensures high precision of storm location, allowing for accurate swell propagation at short time scales away from the storm centres (Komen et al. 1994). Swell waves are lower in magnitude than wind waves; therefore, errors associated with wind-sea waves are more significant than those for swells. For example, Southern Ocean swells dominate Australia’s North-West Shelf, and wave buoy verification shows smaller RMSE and mean absolute error (MAE; MAE=n−1∑in|mi−oi|) than all of Australia (Table 3). Pathirana et al. (2023) quantified swell presence in over 75% of the oceans, which can explain the lag in lead time in forecast error. The rate of error accumulation is higher in the NWP model than in the wave model.

3.2. Site-based verification

Altimeter verification presented in previous sections demonstrated skill improvements when applying bias correction in the new model configuration for open-ocean sea states. However, most users of operational wave forecasts are interested in the coast, and most wave buoys around Australia are located near the coast, as shown in Fig. 6. The Southern Ocean Flux Station, SOFS (Schulz et al. 2012), is omitted in Fig. 6 due to its location ~640 km south-south-west of Tasmania. We included US wave buoys from the North Pacific to verify the model skill in the open ocean. Open ocean refers to sites with depths greater than half the mean wavelength and a few kilometres distant from the coast. In deep-water locations, the wave field is unaffected by coastline geometry and the processes associated with limited depth and shallow water.

Fig. 6.

Location of in situ wave buoys used for verification. There are 69 sites around Australia (left panel) and 14 in the North Pacific (right panel). Solid blue shading indicates the width of the continental shelf (i.e. areas of less than 75-m depth).

Wave hindcast skill across all in situ wave buoys listed in Table 3 and Fig. 7. The results are close to identical for G3 and G4. A closer look at the metrics in Table 3 shows almost no difference, and the results agree with the altimeter verification presented earlier. A global RMSE of 0.29 m is a slight improvement from that of G3 for the same period. Compared to G3, G4 improves in the tail of the probability of exceedance – a small but noteworthy improvement. The ERA5 reanalysis compares well in terms of the shape of the probability density function, as shown in the second column panels in Fig. 7. The ERA5 shows an increase in scatter in H_s in the range up to ~2.5 m. The frequency of extreme H_s appears to be underrepresented in ERA5 and scatter density plots indicate that the bulk of scatter above 2.5 m falls below the line of best fit with a regression coefficient of the least squares slope forced through the origin (fit₀) of 0.972. Table 3 also indicates that the skill in ERA5 wave reanalysis is not far behind that of the G4 model for wave buoys located in the North Pacific and the North-West Shelf of Australia.

Fig. 7.

Wave buoy verification for significant wave height. Statistics are for the AUSWAVE-G3 operational model (top), the configuration for AUSWAVE-G4 (centre) and ERA5 (bottom). Statistics are averaged between 55°N and 55°S from 1 August 2023 to 1 January 2023 and 22 February to 1 August 2024. Panels show the probability of exceedance (left), probability density function and scatter density (right). Statistics in legend are the number of observations (N), Pearson’s correlation coefficient (R), root-mean-square error (RMSE), bias (B), scatter index (SI) and least squares slope forced through the origin (fit₀).

One interesting feature in Table 3 is that the ERA5 reanalysis features a negative bias at high latitudes, i.e. North Pacific wave buoys. By contrast, the opposite is true at low latitudes, i.e. Australian wave buoys. The bottom panel in Fig. 4 confirms this feature and shows the spatial distribution of biases. One can only speculate on the source of these biases. The spatial winds bias shown in Fig. 2 does not support such a trend.

The wave model (AUSWAVE) and wave reanalysis have different resolutions. The spatial resolution of AUSWAVE is 0.125°, with refinements of 0.0625° around islands and continental shelves. The ERA5 wave reanalysis fields are archived at 0.50°, whereas the wave model resolution is 0.36° (Hersbach et al. 2020). Bulk statistics in Table 3 omit wave buoy locations that fell within the land mask of the reanalysis and exclude locations mapped too far offshore. Offshore locations tend to increase the variability, i.e. standard deviation, and introduce a representation bias in the model. As a result, there is a divergence in the number of data points used in buoy verification.

Bulk metrics in Table 3, such as RMSE and R, do not adequately reflect the model’s performance at the distribution’s tail, as shown in the probability of exceedance in Fig. 7. For illustrative purposes, we plotted wave buoy time series for two stations with frequent extreme waves. Figures 8 and 9 show the time series of H_s and wave period for SOFS (Schulz et al. 2012), located south-south-west of Tasmania (Southern Hemisphere), and NOAA wave buoy 46070, located in the south Bering Sea (Northern Hemisphere). The agreement between the models and observations is remarkable for H_s. The SOFS observations show more scatter than the NOAA buoy. The wave model skill at the SOFS location is the lowest in the verification, highlighting the challenges in wave modelling in the Southern Hemisphere with an ocean subject to infinite fetch. The onset of peaks in the time series matches very well for both models, resulting in a R greater than 0.93. Nonetheless, the peaks in H_s are higher in the G4 model. Surprisingly, this is not just the case for H_s above 10 m, as indicated in the probability of exceedance plots. A similar pattern for H_s exists below 10 m, where G4 peaks exceed those of ERA5. Local minima in the time series match between models, concluding that the G4 model has greater variability in H_s than the ERA5 wave reanalysis, and Table 3 shows that the standard deviation of G4 is persistently above the standard deviation of ERA5.

Fig. 8.

Time series of significant wave height (H_s) observed by wave buoys and predicted by AUSWAVE-G4 and ERA5. The panel shows (from top) the Southern Ocean Flux Station (141.3°E 47.2°S) and NOAA wave buoy 46070 (south-west Bering Sea; 175.3°W 55.1°N).

Fig. 9.

Time series of wave period observed by wave buoys and predicted by AUSWAVE-G4 and ERA5. The panel shows the mean zero-crossing periods (T_z) at the Southern Ocean Flux Station (141.3°E 47.2°S) and the peak wave periods (T_p) at NOAA wave buoy 46070 (south-west Bering Sea; 175.3°W 55.1°N).

Verification for the wave periods is more challenging than for significant wave height. Owing to data processing, there is inherently more scattering in the buoy observations. Fig. 9 shows an increase in scatter in the observations for zero-crossing and peak wave periods compared to H_s (Fig. 8). There is better agreement between the models than between the model and buoy observations. An increase in scatter yields a reduction in the R between the model and observations for the wave periods, with values of ~0.7–0.8. The similarity within the wave models is unsurprising since the models share a standard spectral resolution with frequencies logarithmically spaced, which differs to the frequency discretisation used in buoy data processing. Another shared aspect is that the wave models use the same methodology to calculate the resonant four-wave interaction. This method is critical for the down-shift of the peak of the wave spectrum, effectively allowing the wave field to outrun the storm and transition from wind waves to swell (Hasselmann 1960).