Wave forecast investigations on downscaling, source terms, and tides for Aotearoa New Zealand

This study evaluates the effects of downscaling, source terms, and tidal interactions on numerical wave forecasts in Aotearoa New Zealand. We utilised a set of three nested domains (from global to regional scale) to examine significant wave height (Hs), first-order mean period (Tm01), and peak wave direction at two coastal locations, Banks Peninsula and Baring Head. Downscaling markedly improved forecast accuracy at Baring Head, a tidally constricted region, reducing Hs forecast error by 28 %. However, improvements at Banks Peninsula were minimal, likely due to its open-coast characteristics, which are adequately represented even by lower-resolution models. This variability was also evident in the Tm01 predictions, with notable improvements in bias reduction through model downscaling, particularly at Baring Head. Using default source term 6 (ST6) parameters generally improved Hs predictions on the west coast but worsened them on the east, indicating a geographical dependency in model performance. Tidal forcing had a small impact on the overall forecast skill, and its impact was mostly noticed at Baring Head, where tides force large variability. However, the tidally driven wave model showed smaller 12 h variability compared to observations. The study underscores the importance of tailored modelling approaches that consider local geographical and hydrodynamic conditions to optimise wave forecasting.

1 Introduction

Understanding local wave variability is an important aspect of most coastal engineering projects (e.g. Camus et al., 2011; Kamphuis, 2020; Kroon et al., 2020). Beyond these studies, there are numerous ecological (e.g. Coppin et al., 2020), infrastructure and logistics (e.g. Camus et al., 2019; Lucio et al., 2024), and coastal safety (e.g. de Vos and Rautenbach, 2019; Altomare et al., 2020) questions that require a thorough description of offshore and coastal waves. Forecasting ocean waves is important for oceanic traffic safety, on- and offshore industrial operations, and recreational activities. Extreme event forecasting (and the associated accuracy in magnitude and timing) has obvious implications for coastal safety and infrastructure. Increased sea level and cyclone activity due to climate change (e.g. Vitousek et al., 2017; Diamond and Renwick, 2015) has implications for infrastructure design, coastal adaptation, and conservation (Toimil et al., 2020). Examples of these include coastal areas having increased susceptibility to wave impact due to rising sea levels (Hannah, 2004; Hannah and Bell, 2012; Hauer et al., 2016) and the potential of increased large-wave events (Albuquerque et al., 2024).

Figure 1 Nested model domains and respective forecasted wave fields for 06:00 UTC 29 June 2021. Colour shade represents significant wave height from each domain. Global model outputs are shown as the most external colour shade. Basin-scale model outputs are shown inside the largest black contour area. The innermost contour highlights the area of the high-resolution model around the mainland of Aotearoa New Zealand. The white arrows show wave height and peak wave direction plotted for every 40th grid point. The red circles are the location of in situ measurements from wave buoys at Banks Peninsula and Baring Head (the northernmost red dot in this figure). The 200 m isobath is shown as a dark red contour.

Currently, several global wave forecasts are freely available, e.g. the Global Forecast System by the National Centers for Environmental Prediction (NCEP) (Tolman et al., 2002), ERA5 by the European Centre for Medium-Range Weather Forecasts (ECMWF) (Hersbach et al., 2018), and WAVERYS by Copernicus Marine Environment Monitoring Service (CMEMS) (Law-Chune et al., 2021). However, these wave forecast systems use low-resolution models (>20 km) that cannot represent complex bathymetry and sometimes fail at correctly simulating large coastal wave events (Fanti et al., 2023). Therefore, regional and local wave forecasts are important.

In Aotearoa New Zealand, the National Institute of Water and Atmospheric Research (NIWA) is tasked to conduct research and develop tools aimed at enhancing the countries' resilience to wave-related and other environmental hazards. Being an island nation, Aotearoa New Zealand has an extensive exclusive economic zone (EEZ) and coastline to manage. NIWA has developed a platform called EcoConnect that generates and disseminates tailored environmental information services in near-real time (Webster et al., 2008; Moore et al., 2022), including operational forecasts of atmospheric, hydrological, storm-tide, and wave conditions, as well as data obtained from various observational installations.

EcoConnect forecasts waves routinely using Wave Watch III^® (WW3) which is a third-generation spectral wave model that uses an energy-based approach to describe the physical processes of wave growth and transformation at oceanic scales (WW3DG, 2019). WW3 is implemented within EcoConnect using a set of three nested grids to appropriately resolve spatial scales in the wave field surrounding the country (Fig. 1). An adequate model resolution is needed for different regions of Aotearoa New Zealand. Open ocean regions, such as the west and east coasts of the South Island, may be adequately simulated using low-resolution models. Conversely, regions with complex coastlines, such as the Cook Strait and Hauraki Gulf, could require higher-resolution models (<4 km) (Fig. 1). In both cases, validation is necessary to identify regions that need increased model resolution and vice versa.

The WW3 models used in EcoConnect all use the same parameter settings derived from a calibration study (Gorman and Oliver, 2018) in which a global wave model was auto-calibrated against satellite altimetry data using an iterative process that finds values for each model parameter which minimises the root mean square difference between model and observed significant wave heights over the global model domain. Initial short-term calibration tests were conducted using two alternative input–dissipation source term packages, namely, ST2 and ST4, in which ST4 performed better and was selected for a 1-year calibration. ST4 physics include developments in ST1 and ST3 and the latter source terms were not tested by Gorman and Oliver (2018). However, ST6 physics (Rogers et al., 2012; Zieger et al., 2015) has not been tested in EcoConnect's WW3 implementations yet. WW3 ST6 is similar to Simulating Waves in the Nearshore (SWAN, Holthuijsen et al., 1993; Ris et al., 1995) model physics, which is meant for coastal applications. ST6 models have shown improved/similar results compared to ST4 experiments during storm passage and ambient conditions (e.g. Kalourazi et al., 2021; Zou et al., 2023; Meucci et al., 2023b) and has also been used for global wave climate simulations (e.g. Meucci et al., 2023a).

Ocean currents also create important variability in the wave field, and they have been an active topic of research in the last few decades (e.g. Zhang et al., 2022). Wave–current interaction studies have shown variability ranging from large-scale oceanic currents (e.g. Barnes and Rautenbach, 2020) to small-scale tidally dominated regions (e.g. Vincent, 1979; Ris et al., 1995). For instance, a wave train propagating against opposing currents tends to increase wave steepness, which in turn can generate wave breaking and dissipation (Holthuijsen, 2007). Opposing currents can also reduce the wave period/length and have been called “wave straining” by Holthuijsen and Tolman (1991). These effects have been observed and simulated in different regions of the globe using short-term time series (e.g. Ardhuin et al., 2012; Rapizo et al., 2017; Halsne et al., 2024). However, a long-term (≥1-year) evaluation of the importance of tidal currents on the wave forecast is still needed.

The aim of the present study is to analyse the impact of wave model downscaling (including storm tide forcing) on wave forecast for Aotearoa New Zealand. The impact of downscaling is assessed using a set of three nested model domains with increasing resolution but with similar grid (i and j) dimensions, thus implying the same computational cost. A comparison between source term physics (ST4 and ST6) is made using the intermediate model grid. The impact of tides and storm surge forcing is analysed using the model with the highest resolution in which two experiments are compared: with and without tidal and storm surge forcing. Both in situ and remotely sensed validation/investigations are performed, and the physical dynamics of Aotearoa New Zealand's wave mean fields are discussed, especially in the Cook Strait. In this research article, we acknowledge the existence of names of geographical locations in Te Reo Māori. An exception to that is Aotearoa New Zealand, which is widely used in that form.

2 Methods

2.1 The operational forecasting system

The EcoConnect platform ingests data and runs numerical models to deliver forecasts for a variety of natural hazards. EcoConnect operates autonomously via the Cylc workflow meta scheduler (Oliver et al., 2018) and starts by downloading the United Kingdom Met Office's global model atmospheric forecast (UKMO). This model is a configuration of the Unified Model (UM, Maher and Earnshaw, 2022), and provides lateral boundary conditions for the New Zealand Limited Area Model (NZLAM) atmospheric model, itself a local configuration of the UM, whose domain extends from eastern Australia to the Chatham Islands (external black contour in Fig. 1). NZLAM provides the initial and lateral boundary conditions for the New Zealand Convective-Scale Model (NZCSM) numerical model, a convection-permitting configuration of the UM covering just the Aotearoa New Zealand's landmass and its coastal waters (innermost domain in Fig. 1). The three domains of the wave models cover the same area as UKMO (global), NZLAM, and NZCSM atmospheric models; more details are provided in the next section.

Both NZLAM and NZCSM are configured to provide input data for a suite of different hazards models. These include a hydrological river flow model, TopNet; forecast streamflow for just under 50 000 river reaches around Aotearoa New Zealand (Cattoën et al., 2022); and a hierarchy of wave and current forecast models based on WW3 and the River and Coastal Ocean Model (RiCOM) (more details about these models can be found in Sect. 2.2 and 2.4). Observation datasets collected and disseminated within EcoConnect include satellite imagery, surface weather station, river gauge, and wave buoy data. These model outputs are generated, processed, compiled, and archived by bespoke tasks in the EcoConnect workflow, all orchestrated by Cylc.

2.2 Wave model

The WW3 version 6.07.1 (WW3DG, 2019) is used for wave forecasts in EcoConnect. The model represents the sea state by the two-dimensional ocean wave spectrum $F(\mathit{k},\mathit{x},t)$, which gives the energy density of the wave field as a function of wavenumber $\mathit{k}=(k_x,k_y)$ at each position $\mathit{x}=(x,y)$ in the model grid and time t of the simulation. The spectrum evolves subject to the radiative transfer equation

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial N}{\partial t}}+{\mathbf{\nabla }}_{(x,y)}\cdot \left(\dot{\mathit{x}}N\right)+{\displaystyle \frac{\partial }{\partial k}}\left(\dot{k}N\right)+{\displaystyle \frac{\partial }{\partial \mathit{\theta }}}\left(\dot{\mathit{\theta }}N\right)={\displaystyle \frac{S}{\mathit{\sigma }}}\end{array}$$

for the wave action $N(k,\mathit{\theta },\mathit{x},t)=F(\mathit{k},\mathit{x},t)/\mathit{\sigma }\left(k\right)$, where the dots represent time derivatives, θ is the propagation direction, and σ=2πf is the relative (radian) frequency associated with waves of wavenumber magnitude k through the linear dispersion relation

$$\begin{array}{}\text{(2)}& \mathit{\sigma }{}^{\mathrm{2}}=gk\tanh kd.\end{array}$$

The frequency σ of waves propagating subject to gravity g and water depth d is observed relative to a frame of reference moving with a mean current U, providing a Doppler shift from the absolute (radian) frequency:

$$\begin{array}{}\text{(3)}& \mathit{\omega }=\mathit{\sigma }+\mathit{k}\cdot \mathit{U}.\end{array}$$

While WW3 internally computes the wavenumber spectrum $F(\mathit{k},\mathit{x},t)$ due to its invariance properties, for output purposes, this is converted to the traditional frequency–direction spectrum

$$\begin{array}{}\text{(4)}& F(f,\mathit{\theta },\mathit{x},t)=\left({\displaystyle \frac{\mathrm{2}\mathit{\pi }}{c_{\mathrm{g}}}}\right)F(\mathit{k},\mathit{x},t).\end{array}$$

The terms on the left-hand side of Eq. (1) represent spatial advection and the shifts in wavenumber magnitude and direction due to refraction by currents and varying water depth. The source term S on the right-hand side of Eq. (1) represents all other processes that transfer energy to and from wave spectral components. It can be expressed as

$$\begin{array}{}\text{(5)}& S=S_{\text{in}}+S_{\text{ds}}+S_{\text{nl}}+S_{\text{ice}}+\mathrm{\dots }\end{array}$$

to include contributions from wind forcing (S_in), energy dissipation (S_ds), weakly nonlinear four-wave interactions (S_nl), scattering and wave–ice interactions (S_ice), and other terms.

2.2.1 Wave forecast implementations

The wave component of EcoConnect consists of three nested implementations of WW3 (WW3DG, 2019):

1. The global model (GLOBALWAVE) implemented on a regular latitude–longitude grid covering longitudes from 0° to 360° at 0.234375° (<26 km) resolution and latitudes from −81.25° to +81.25° at 0.15625° (∼17 km) resolution (Gorman and Oliver, 2018). Atmospheric inputs are provided by the UKMO (Maher and Earnshaw, 2022). GLOBALWAVE forecast operates twice daily (00:00 and 12:00 UTC) to provide 6 d forecasts. Daily ice concentration fields are also sourced from the UKMO. No current or sea level inputs are used. A global time step of 900 s is used, with minimum time steps of 180 s for spatial advection and source term integration, and 900 s for refraction.

2. The regional model (NZWAVE) is implemented on a regular latitude–longitude grid covering longitudes from 143.3203125° to 184.5703125° at 0.05859375° (<6 km) resolution and latitudes from −54.45° to −20.85° at 0.0390625° (∼4 km) resolution (Fig. 1). Atmospheric inputs are provided by NIWA's deterministic forecast implementation of the Unified Model for the Tasman Sea and Aotearoa New Zealand (NZLAM; see Sect. 2.3 for details), with both models running four times daily (00:00, 06:00, 12:00, and 18:00 UT) 72 h forecasts. No ice, current, or sea level inputs are used. Spectral boundary conditions for NZWAVE are sourced by GLOBALWAVE. A global time step of 300 s is used, with minimum time steps of 60 s for spatial advection and source term integration and 300 s for refraction.

3. The mainland Aotearoa New Zealand model (NZWAVE-HR) is implemented on a regular latitude–longitude grid covering longitudes from 163.21° to 181.67° at 0.029296875° (<3.25 km) resolution and latitudes from −48.54° to −30.84° at 0.019531250° (∼2 km) resolution (Fig. 1). Atmospheric inputs are provided by NIWA's deterministic convection-resolving forecast implementation of the Unified Model (NZCSM; see Sect. 2.3 for details), with both models running four times daily (00:00, 06:00, 12:00, and 18:00 UT) 48 h forecasts. Storm surge forecasts of sea level and current fields (NZSURGE-HR, forced by NZCSM) are combined with tidal sea levels and currents derived from the NZTIDE harmonic tidal model to provide input fields for NZWAVE-HR (more details about NZSURGE-HR and NZTIDE in Sect. 2.4). Spectral lateral conditions for NZWAVE-HR are provided by NZWAVE. A global time step of 180 s is used, with minimum time steps of 60 s for spatial advection, 30 s for source term integration, and 180 s for refraction.

2.3 Atmospheric forcing

Atmospheric forcing for the wave models are derived from UKMO (global model, Maher and Earnshaw, 2022) and NIWA's family of numerical weather prediction (NWP) models in EcoConnect. UKMO provides 10 m winds and sea ice concentration for GLOBALWAVE, whereas NZWAVE and NZWAVE-HR are forced by near-surface winds from the NIWA's regional models. NIWA operates two limited area NWP models: i) the NZLAM, which featured a 12 km horizontal resolution between 2007 and late 2019 and since then a 4.4 km horizontal resolution, and ii) the NZCSM, which runs with a 1.5 km horizontal resolution. NZLAM's domain is represented by the larger black outline in Fig. 1, and NZCSM's domain is the smaller outline in Fig. 1, covering Aotearoa New Zealand's main landmass and its coastal waters.

Both models are based on the UM, a non-hydrostatic, fully compressible, deep-atmosphere model whose dynamical core, ENDGame (Even Newer Dynamics for General atmospheric modelling of the environment), solves the equations of motion using mass-conservation, semi-implicit, semi-Lagrangian, time-integration methods (Wood et al., 2014). However, they feature differing science configurations and workflow setups. NZLAM's workflow includes a three-dimensional variational (3D–Var) data assimilation method which ingests observations from satellites, aircraft, ships, buoys, and land surface synoptic weather stations and is configured to use the Met Office GA6 Global model science settings (Walters et al., 2017). NZLAM forces NZWAVE with outputs of 1 h temporal resolution, and the forecast extends 72 h into the future. It runs four times a day, at 00:00, 06:00, 12:00, and 18:00 UTC, generating analyses at each cycle.

NZCSM is a convection-permitting model with minor changes to NZLAM configuration that improve model stability to run over Aotearoa New Zealand's complex terrain with high resolution in a smaller domain. The scientific configuration of NZCSM is equivalent to the setup described in Bush et al. (2020). Like NZLAM, NZCSM is warm-cycled, it restarts from an output of the previous forecast. However, NZCSM does not perform its own data assimilation. Instead, at the start of each forecast cycle, the larger-scale analysis from NZLAM, which has benefited from its data assimilation, is merged with the forecast from the previous NZCSM cycle to give an improved atmospheric state from which another forecast can begin. NZCSM's lateral boundary conditions are derived from NZLAM with a 20 min update interval and operates four times daily, on the 00:00, 06:00, 12:00, and 18:00 UTC analysis cycle, forecasting 48 h ahead. Forecast outputs from NZCSM, including the driving data for the downstream wave models and other components in EcoConnect, are made available at 30 min temporal resolution. Previous studies have found NZCSM wind speeds to have errors in the range of 2 %–15 % of observed values in the Cook Strait region of Aotearoa New Zealand. The highest errors typically occur during more extreme events (Yang et al., 2017). The version of NZCSM used in this has improved results over the simulations evaluated in Yang et al. (2017) (not shown).

2.4 Water level and velocity forcing

Tidal and storm surge (infra-inertial) variability of water level and depth-averaged velocity are predicted using RiCOM (River and Coastal Ocean Model), a semi-implicit semi-Lagrangian finite-element model based on an unstructured triangular grid that can be run as both a harmonic solver and a time-stepping hydrodynamic solver (Walters, 1992, 2005, 2006). Tides (NZTIDE) and storm surge (NZSURGE and NZSURGE-HR) are calculated separately within EcoConnect for increased computation efficiency. Simulating barotropic tides would require smaller time steps; instead, they are resolved harmonically, saving computational time. These physical components are then summed to give total water level and velocity.

The hydrodynamic part of the RiCOM solves the Reynolds-averaged Navier–Stokes equations assuming the incompressibility condition. A semi-implicit time-stepping scheme is used to solve advection using a semi-Lagrangian algorithm (Staniforth and Côté, 1991), including a power series particle tracking method (Walters et al., 2007). The Coriolis term is added explicitly using a third-order Adams–Bashforth scheme (Walters et al., 2009).

Two storm surge forecasts (NZSURGE and NZSURGE-HR) are included in EcoConnect, both running four times daily using the time-stepping RiCOM model forced by 10 m winds and mean sea level pressure but sourcing these inputs from the 72 h NZLAM and 48 h NZCSM forecasts, respectively. The new atmospheric forcing file is smoothed from the previous file over the first 6 h of the forecast to account for changes in the atmospheric variables due to data assimilation. Wind forcing is included in the model as surface stress using a quadratic formulation based on wind velocity with a drag coefficient specified by Wu (1982). Surface pressure is used to calculate an inverse barometer surface level, which is applied as a loading term. Relaxing to an inverse barometer water level within a radiation boundary condition is applied as the lateral open boundary conditions. Further details on the NZSURGE operational model can be found in Lane et al. (2009) and verification of the results in Lane and Walters (2009).

The tidal model (NZTIDE) is calculated using a harmonic tidal model formulation of RiCOM (Walters et al., 2001) based on a harmonic decomposition in time and a finite-element approximation in space. This model provides amplitudes and phases for the eight largest tidal constituents around Aotearoa New Zealand: M2, N2, S2, K2, K1, O1, P1, and Q1. Dependent variables are expressed in terms of harmonic expansions of these constituents with the nonlinear bottom friction included as a series expansion (Walters, 1992). Equilibrium tide and self-attraction/loading tide are also included in the formulation (Goring and Walters, 2002). The results of the tide model are interpolated from the unstructured grid onto a regular Cartesian grid for simplified viewing and outputs. For each forecast period, the tidal sea surface and currents are reconstructed from these tidal constituents. Tides in Aotearoa New Zealand are dominated by the M2 constituent, which proceeds around the two main islands in an anticlockwise direction. Tidal flows through Cook Strait (Fig. 1), the narrow strait between the North Island / Te Ika A Maui and South Island / Te Waipounamu, can be especially strong because the tide levels are close to 180° out of phase at either end. NZTIDE and NZSURGE-HR water levels and currents are summed and interpolated onto the NZWAVE-HR domain to provide total water level and velocity forcing for that wave model.

A new model of tides and storm surges has been developed using TELEMAC. This new finite-element model simulates tides and storm surge simultaneously using a time-stepping approach (for technical details about the model please refer to Hervouet, 2000; Moulinec et al., 2011). This new model is forced by the UKMO global atmospheric model. For comparison against RiCOM results (NZTIDE and NZSURGE-HR), we extracted water level and velocity output from this model at the Baring Head station for wave–current interaction analysis. This new model is planned to replace RiCOM as the operational water level and velocity forecasting model in the near future. In this work, we provide initial comparisons which might highlight the importance of the migrating to this new model.

2.5 Observations

2.5.1 Satellite data

Near-real-time, gridded, satellite derived daily average significant wave heights from CMEMS were used to validate the forecasts between 1 January and 31 December 2021. This gridded product includes daily mean and maximum significant wave height data from different altimeter missions. This product merges along-track data from Jason-3, Sentinel-3A, Sentinel-3B, SARAL/AltiKa, Cryosat-2, CFOSAT, and HaiYang-2B missions and delivers daily data at 2° spatial resolution with an uncertainty ranging from 0.12 to 0.44 m (CMEMS, 2024, last access: 11 April 2024).

2.5.2 Buoy data

Two in situ sites with near-real-time wave observations are used to validate the predicted wave height, mean period, and peak direction. The Banks Peninsula wave buoy is a directional wave rider moored approximately 17 km east of Banks Peninsula in the South Island / Te Waiponamou at a latitude of 43°45^′ S and longitude of 173°20^′ E at approximately 80 m of water depth (Fig. 2a). The mooring location has been continuously maintained since 1999 with data gaps limited to buoy or mooring failures. The location is exposed to a wide range of swells from the northeast to the south. The mean significant wave height over the observation period is 2.1 m with an average mean period (Tm02) of 6.5 s (Walsh, 2017).

Figure 2 Model bathymetries at Banks Peninsula (a) and Baring Head (b) showing locations of wave buoy measurements (red dots). NZWAVE-HR bathymetry is displayed using the colour shading and green contour. The 30 and 80 m isobaths from NZWAVE (blue) and GLOBALWAVE (black) are also shown for comparison.

The Baring Head station has been instrumented since 1995 with a non-directional wave rider buoy at the beginning and with directional wave rider buoys since 2014. The site is located at approximately 44 m depth, 2 km west of Baring Head lighthouse and 15 km southeast of Wellington / Te Whanganui-a-Tara (Fig. 2b). The site is located inside the Cook Strait and is sheltered from many wave directions by both South Island / Te Waiponamou and North Island / Te Ika A Maui. Most swells come from the south, but the location is also exposed to a narrow swell window from the west/northwest (Allis et al., 2021).

Observations from 1 January 2021 to 31 December 2021 from the two stations are used in this study. Before this period, both stations switched from the second (Tm02) to the first-order mean period (Tm01). These mean periods represent ratios with regard to the zeroth spectral moment and provide different estimates of the average wave period depending on the spectrum's shape. Time series of significant wave height, first-order mean period (Tm01), and peak direction are used for model validation. Wave observations are available every 15 min and are filtered using a 1 h moving average filter. Significant wave height and peak direction are decomposed into u and v components, filtered (1 h moving mean) and then recomposed to remove spiky variability in both time series. The locations of wave buoy measurements are shown in Fig. 2, as well as NZWAVE-HR bathymetry and GLOBALWAVE and NZWAVE 30 and 80 m isobaths. Here, one can expect that the bathymetry resolution will generate differences in how wave refraction will occur in the nearshore, especially for long-period swells.

2.6 Experiment design and evaluation metrics

Five numerical wave forecasts are evaluated in this study. Three simulations are part of EcoConnect's forecast: GLOBALWAVE, NZWAVE, and NZWAVE-HR. The latter one includes tidal and storm surge forcing via water level and current inputs to the model – hereinafter tides (Table 1). All models use a calibrated version of ST4 physics (Gorman and Oliver, 2018). The operational forecasts have run in the current cluster since 2018 and are analysed in this study.

Table 1 Wave forecast experiments, their domains, forcings, lateral conditions and source term configurations.

Two additional forecasts are run to address the questions raised in this study. In one simulation, we use ST6 physics with default parameters (Table 2.8 in WW3DG, 2019) applied to the intermediate domain, and the forecast experiment is called NZWAVE-ST6 (Table 1). This experiment sheds light on ST6 performance around Aotearoa New Zealand waters for later consideration of these source terms in EcoConnect's forecasts. The last experiment removed tides and storm surge forcing from the highest-resolution wave model and is called NZWAVE-HR-NOTIDES (Table 1). Comparisons between NZWAVE-HR-NOTIDES and NZWAVE-HR show the importance of ocean currents as forcing to generate wave variability in coastal regions of Aotearoa New Zealand. The two additional forecasts start on 1 December 2020 from rest and they have a spin-up period of 1 month – similarly to Gorman and Oliver (2018).

We use 24 h forecasts starting daily at 00:00 UTC. Model results are interpolated to the satellite observational field, and the closest model grid point is used for comparison against in situ data. Forecast outputs are written half-hourly for the highest-resolution models (NZWAVE-HR and NZWAVE-HR-NOTIDES) and hourly for the rest of the models. Model outputs are linearly interpolated to the observational frequency before model–data comparison. Model significant wave height daily averages are validated using satellite data. Model time series of significant wave height (Hs), peak period (Dp), and first-order mean period (Tm01) are also validated at the Banks Peninsula and Baring Head coastal stations (Fig. 2). The forecasts are objectively validated using the root mean square error (RMSE) given by

$$\begin{array}{}\text{(6)}& \text{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{n}}{\mathrm{\Sigma }}_{i=\mathrm{1}}^n(x_i-y_i{)}^{\mathrm{2}}}\end{array}$$

and linear correlation (r) given by

$$\begin{array}{}\text{(7)}& r={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n(x_i-\stackrel{\mathrm{‾}}{x})(y_i-\stackrel{\mathrm{‾}}{y})}{\sqrt{{\sum }_{i=\mathrm{1}}^n(x_i-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}}\sqrt{{\sum }_{i=\mathrm{1}}^n(y_i-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}}}}\end{array}$$

between observed (x) and predicted (y) results, where $i=\mathrm{1},\mathrm{2},\mathrm{\dots },n$ are the observation times or locations and the averages (⁻) are applied in time or space. Spatial fields of significant wave height RMSE and bias (model − observations) are computed using daily averaged model results interpolated to satellite observations horizontal grid.

The role of currents in wave refraction is further analysed using the vertical component of the relative vorticity (ζ), which is defined as $\mathit{\zeta }=\partial v/\partial x-\partial u/\partial y$, where x and y (u and v) are the zonal and meridional components of space (velocity). According to Dysthe (2001), the curvature of a wave ray (χ) travelling over variable currents can be defined as $\mathit{\chi }=\mathit{\zeta }/v_{\mathrm{g}}$, where v_g is the wave group velocity. Using this equation, one expects positive relative vorticity to refract waves (or bend wave rays) in one direction and negative relative vorticity to shift waves in the opposite direction.

Harmonic analysis is applied to in situ observations and NZWAVE-HR forecast time series of significant wave height, mean period, and peak direction using t_tide (Pawlowicz et al., 2002). This allows us to identify tidal oscillations in these wave parameters.

3 Results

This section analyses the impacts of downscaling, source terms, and tides on forecasting significant wave height, wave peak direction, and mean period (Tm01). Initially, the historical 24 h forecasts are evaluated against satellite data and in situ observations from Banks Peninsula (open coast) and Baring Head (constricted region) using RMSE, linear correlation analysis, and bias metrics. This section is concluded by highlighting the impacts of the tides and storm surge on the wave parameters, especially at Baring Head.

3.1 Mean spatial fields

Mean bias and RMSE analyses between altimeter data and the wave forecasts are given in Figs. 3 and 4. The 1-year average of satellite observations indicates an area with large significant wave height (∼4.5 m ) below 48° S, which penetrates the Tasman Sea and offshore regions of southeast Aotearoa New Zealand (Fig. 3a). The northeast region of Te Ika a Maui / North Island is marked with the smallest average waves due to the landmass barrier that protects the region from south/southwest swells. All forecasts using calibrated source term 4 (ST4) (Gorman and Oliver, 2018) show a relatively larger positive mean bias of significant wave height on the west side of Aotearoa New Zealand, whereas the bias is small on the eastern side of the country (Fig. 3b, c, e, and f). The spatial mean bias in the smallest domain area (black contour in Fig. 3) varies from 0.10 to 0.13 m between all ST4 forecasts. Applying ST6 physics (NZWAVE-ST6) with default parameters reduces the positive bias on the western side of Aotearoa New Zealand compared to the same model using calibrated ST4 (NZWAVE) (Fig. 3c and d). However, it also produces a negative bias on the eastern side of the country, generating a mean bias of −0.07 m computed for the smallest domain. Removing tides from the forecast (NZWAVE-HR-NOTIDES) shows no significant change in the mean bias field analysis (Fig. 3e and f), implying that tides and storm surge forcing do not significantly affect the temporal mean values of these wave parameters.

Figure 3 Maps of mean observed significant wave height (m) (a) and forecast bias (m) from GLOBALWAVE (b), NZWAVE (c), NZWAVE-ST6 (d), NZWAVE-HR-NOTIDES (e), and NZWAVE-HR (f) models. Statistics are computed between 1 January and 31 December 2021. The spatial means are shown in the figure title and highlighted over defined regions.

Figure 4 Maps of observed significant wave height standard deviation (a) in metres and forecast RMSE (m) from GLOBALWAVE (b), NZWAVE (c), NZWAVE-ST6 (d), NZWAVE-HR-NOTIDES (e), and NZWAVE-HR (f) models. Statistics are computed between 1 January and 31 December 2021. The spatial means are shown in the figure title and highlighted over defined regions in metres.

In Fig. 4, 1-year significant wave height standard deviation from altimeter data is shown alongside the RMSE comparisons with the various historical forecasts described in Sect. 2.2. There is considerable variability in wave height south of Tasmania and South Island / Te Waipounamu of Aotearoa New Zealand in which large standard deviations of significant wave height (∼1.5 m) are observed (Fig. 4a). This region of large variability fades towards the north and reaches minima values (∼0.75 m) above 30°S and on the eastern side of Aotearoa New Zealand. GLOBALWAVE and NZWAVE show a similar RMSE in the whole analysed region and for the smallest domain (Fig. 4b and c). This is associated with the same ST4 parameters used in both models. Model resolution does not significantly impact these results, and GLOBALWAVE, NZWAVE, and NZWAVE-HR-NOTIDES have similar mean RMSEs. However, satellite observations have low spatial resolution (2°), and some coastal points show larger RMSEs, which are associated with localised larger uncertainty in the satellite observations.

Applying ST6 default parameters (NZWAVE-ST6) reduces model RMSE south of Tasmania and minimises the mean RMSE from 0.43 to 0.41 m compared to NZWAVE for their whole domain. NZWAVE-ST6 also decreases the significant wave height RMSE northwest of the North Island but increases it on the east side of the South Island. In general, ST6 physics reduces the mean RMSE in the NZWAVE-HR region by 0.01 m (Fig. 4c and d). Including tides as forcing in the wave model (NZWAVE-HR) slightly increased (0.01 m) the average RMSE of significant wave height.

3.2 In situ significant wave height (Hs)

Model evaluation against in situ significant wave height (Hs) shows a similar RMSE (0.30–0.32 m) and correlation coefficient (0.92) between the ST4 sensitivity experiments at Banks Peninsula – an open-coast station (Fig. 5). However, relatively large differences are found when analysing the Hs bias. A reduction in the negative bias from −0.18 m (GLOBALWAVE) to −0.08 m (NZWAVE) and an even further reduction (−0.04 m) are found in NZWAVE-HR-NOTIDES (Fig. 5a, b, and d). Including tides in the simulation (NZWAVE-HR) does not generate a marked impact on evaluation metrics (RMSE, bias, and r) and shows a slight degradation in bias (Fig. 5d and e). Using standard ST6 parameters (NZWAVE-ST6) degraded the forecast bias (−0.20 m) and RMSE (0.38 m) compared to the calibrated ST4 forecast (NZWAVE) (Fig. 5b and c).

Figure 5 Scatter plots of observed and forecasted significant wave heights from GLOBALWAVE (a), NZWAVE (b), NZWAVE-ST6 (c), NZWAVE-HR-NOTIDES (d), and NZWAVE-HR (e) simulations at the Banks Peninsula station. Root mean square error (RMSE), mean bias, and coefficient of correlation (r) are shown in the top-left corner of each panel. The scatter colour represents the probability density estimate.

Significant wave heights smaller than 2 m have larger probability density estimates in the forecasts and observations. Wave heights between 2 and 6 m tend to have an even distribution with regard to over- and underestimation of the observed values. Some events with a wave height between 4 and 6 m are overestimated by the forecasts and fall within the 6 to 8 m category. Within the observed 6 to 8 m category, most forecasted waves are within the same range, except for a few data points that fall into the 4 to 6 m category. The largest waves (∼8 m) were underestimated by GLOBALWAVE and NZWAVE-ST6.

At Baring Head, a constricted coastal region, larger improvements in the evaluation metrics are generated by the downscaling in comparison to Banks Peninsula (Fig. 6). A decrease of 25 % (0.09 m) in the significant wave height RMSE is shown from the lowest-resolution (GLOBALWAVE) to the highest-resolution model (NZWAVE-HR-NOTIDES) (Fig. 6a and d). The intermediate grid generates a 17 % (0.06 m) decrease in RMSE compared to GLOBALWAVE (Fig. 6a and b). However, applying standard ST6 parameters to the intermediate model domain (NZWAVE-ST6) reduced the forecast skill (RMSE and r) back to the level of the GLOBALWAVE forecast (Fig. 6a and c). Including tides and storm surge forcing into the highest-resolution model (NZWAVE-HR) generates a slight improvement in the forecast RMSE and correlation coefficient (Fig. 6d and e).

Figure 6 Scatter plots of observed and forecasted significant wave heights from GLOBALWAVE (a), NZWAVE (b), NZWAVE-ST6 (c), NZWAVE-HR-NOTIDES (d), and NZWAVE-HR (e) runs at Baring Head station. Root mean square error (RMSE), mean bias, and correlation coefficient (r) are shown in the top-left corner of each panel. The scatter colour represents the probability density estimate.

Significant wave height forecast bias switches from positive (GLOBALWAVE =0.13 m) to negative (NZWAVE $=-\mathrm{0.12}$ m) in the first downscale exercise, but it is improved (−0.05 m) in the highest-resolution model (NZWAVE-HR-NOTIDES). The largest positive bias in GLOBALWAVE might be related to the reduced levels of wave height dissipation in the region due to the poor representation of the coastline and bathymetry near Baring Head (Fig. 2b). The negative wave height bias in NZWAVE (−0.12 m) is further reduced in NZWAVE-HR-NOTIDES (−0.05 m), which better represents the underwater topography and coastline.

Wave heights tend to be smaller (∼1 m) at Baring Head compared to Banks Peninsula (∼2 m). Its probability density estimate is the largest near 1 m wave height in all forecasts. The region with the highest probability density is located around the 1 : 1 line for most forecasts except GLOBALWAVE. The latter forecast has a larger density estimate above the 1 : 1 line which represents overestimation for most of those small-wave events.

3.3 In situ mean period (Tm01)

Model evaluation analysis of the first-order mean period (Tm01) shows root mean square error (RMSE) and absolute bias smaller than 0.8 and 0.3 s, respectively, for all forecasts at the Banks Peninsula station (Table 2). All models had similar and high correlation coefficients with observations of Tm01 (>0.8). The most marked changes occur between NZWAVE-ST6 and NZWAVE. A small degradation of the forecast RMSE and correlation coefficient is found when comparing NZWAVE-ST6 (RMSE=0.77 s, r=0.85) and NZWAVE (RMSE=0.65 s, r=0.89). However, a smaller absolute bias is simulated by the former model (−0.04 s) compared to the latter forecast (−0.14 s). A Tm01 gradual bias reduction is found from the lowest-resolution (GLOBALWAVE, $\text{bias}=-\mathrm{0.28}$ s) to the highest-resolution forecast without tides (NZWAVE-HR-NOTIDES, $\text{bias}=-\mathrm{0.10}$ s). Including tides and storm surge forcing slightly degraded the Tm01 bias (NZWAVE-HR, $\text{bias}=-\mathrm{0.11}$ s) compared to the same model without a varying water level and current forcing (NZWAVE-HR-NOTIDES, $\text{bias}=-\mathrm{0.10}$ s).

Table 2 First-order mean wave period (Tm01) skill evaluation of different forecasting wave models at the Banks Peninsula and Baring Head stations. Tm01 evaluation analysis include root mean square error (RMSE), bias (model – obs.), and correlation coefficient.

Model evaluation statistics of the mean period show a larger error at Baring Head compared to Banks Peninsula. Maximum Tm01 forecast RMSE (NZWAVE-ST6 =1.75 s) more than doubled when comparing the two stations (Table 2). GLOBALWAVE has the smallest Tm01 forecast RMSE (1.46 s) and the largest correlation coefficient (r=0.85). The downscaled models (NZWAVE and NZWAVE-HR-NOTIDES) increased the RMSE and reduced the correlation coefficient compared to GLOBALWAVE. However, a gradual improvement in the Tm01 bias from GLOBALWAVE (−0.78 s) to NZWAVE (−0.67 s) and NZWAVE-HR-NOTIDES (−0.35 s) is seen. Moreover, adding tidal and storm surge forcing in NZWAVE-HR slightly reduced the RMSE (1.63 s) and absolute bias (−0.34 s).

3.4 In situ peak wave direction

At the Banks Peninsula station, the observed peak wave direction is largely dominated by the south component (∼30 % of occurrence), with the significant wave height reaching more than 6 m (Fig. 7a). Large waves (>6 m) also arrive from the SSW quadrant; however, this direction has a smaller percentage of occurrence (∼5 %). The second-most-frequent direction is the SSE (∼15 %) and is followed by an eastern component (∼13 %). This eastern component, however, is dominated by smaller waves.

Figure 7 Directional histograms (wave roses) of peak wave direction and significant wave height from buoy measurements (a), GLOBALWAVE (b), NZWAVE (c), NZWAVE-ST6 (d), NZWAVE-HR-NOTIDES (e), and NZWAVE-HR (f) runs at the Banks Peninsula station.

The numerical simulations show a wider spread in larger waves arriving from the S and SSW quadrants at the Banks Peninsula station (Fig. 7). That differs from the observations which have the largest frequency of occurrence more concentrated on the S component. GLOBALWAVE and NZWAVE-ST6 have an even distribution between those two directions, whereas NZWAVE shows a slightly larger preferential distribution towards the S component. The highest-resolution models, however, show larger occurrences of the SSW direction component. All simulations show a marked northeast component, with the percentage of occurrence ranging from ∼16 % (GLOBALWAVE) to ∼21 % (NZWAVE) and wave height reaching values above 4 m. This large-wave (>4 m) northeastern component is also observed in the buoy measurements, however, with a smaller percentage of occurrence (∼12 %).

Peak wave direction observations from the SSE have the largest frequency of occurrence at the Baring Head station (Fig. 8a). This might be associated with strong tidal currents that occur in the region with a marked southeastward residual component and its associated relative vorticity (Fig. 7 of  Walters et al., 2010), which steer incoming southerly swells towards a similar current direction (more details on the process are discussed in Sect. 3.5). Waves from the south quadrant are the second-most-frequent ones. Both the SSE and S direction can have the largest significant wave heights (>6 m). In contrast, the NW (third-most-frequent) component is marked by smaller waves (<1 m) locally generated by the strong NW winds that often happen in the region (Reid, 1996).

Figure 8 Directional histograms (wave roses) of significant wave height from buoy measurements (a), GLOBALWAVE (b), NZWAVE (c), NZWAVE-ST6 (d), NZWAVE-HR-NOTIDES (e), and NZWAVE-HR (f) runs at Baring Head station.

All forecasting models show a prevalence of waves coming from the south (Fig. 8). The percentage of occurrence ranges from 52 % (GLOBALWAVE) to 58 % (NZWAVE-HR). The second-most-common peak wave direction is the SSW. Both S and SSW components can have waves larger than 4 m, but waves between 1 and 2 m are the most frequent. NZWAVE is the only model that has its third-largest component associated with the NW direction – similarly to observations. GLOBALWAVE and NZWAVE-ST6 have NNW as their third-most-frequent direction. The highest-resolution models (NZWAVE-HR and NZWAVE-HR-NOTIDES), however, show the SE component to be the third most frequent. The lack of a NW/NNW component in those models might be associated with weaker northwesterlies in the atmospheric forcing (NZCSM), which are not able to generate the observed smaller (<1 m) and frequent waves in the region or the model configuration itself.

3.5 Tidal influence on wave height, period, and direction

A close look at intra-daily variability shows high-frequency oscillations in the observed significant wave height and mean period (Tm01), which are often matched by ∼12 h peaks in predicted wave height and period by NZWAVE-HR (green and blue dots in Fig. 9a and b). This 12 h variability simulated by NZWAVE-HR seems to be generated by the interaction between waves and tidal currents. The inclusion of tides as forcing, however, generated only a small reduction (increase) in the forecast RMSE (correlation coefficient) with respect to Hs observations compared to the simulation without tidal forcing (NZWAVE-HR-NOTIDES) for the period shown in Fig. 9. This wave–current interaction has been observed and simulated in different regions of the globe (e.g. Ardhuin et al., 2012; Rapizo et al., 2017; Barnes and Rautenbach, 2020; Halsne et al., 2024). Wave buoy measurements have additional high-frequency oscillations (green dots in Fig. 9a and b) not accounted for by NZWAVE-HR due to its limitations in model approximations and/or in its forcings. In this study, the forecast model without tidal and storm surge forcing (NZWAVE-HR-NOTIDES) does not show any 12 h variability in significant wave height and mean period (Tm01) and showed a slightly lower skill on forecasting Hs compared to NZWAVE-HR (Figs. 9 and 11).

Figure 9 Time series of significant wave height (a), first-order mean period (Tm01) (b), and depth-averaged meridional currents (c) at Banks Peninsula from buoy measurements (green), NZWAVE-HR-NOTIDES (magenta), and NZWAVE-HR (blue). The RMSE, correlation coefficient, and bias are presented next to each experiment label in the legend. Those statistical metrics were computed for the period shown in the figure.

At Banks Peninsula, NZWAVE-HR simulates this 12 h variability that occurs due to the interaction between southerly waves and the anticlockwise propagation of the tidal wave around Aotearoa New Zealand's continental shelf. Around midday and midnight on 21 October 2021, southward (negative) tidal currents flow against a swell propagating northeastward (20°) (red rectangles in Fig. 9). These counter-currents increase wave height and reduce the wave period every tidal cycle near low tide, which propagates on the continental shelf as a progressive wave – peaks/troughs in water level match peaks/troughs in tidal currents (Walters et al., 2001). This growth in significant wave height and shortening in wave period while facing opposing currents is well predicted in theory and observed and simulated in different regions (Phillips, 1977; Vincent, 1979; Ardhuin et al., 2012; Rapizo et al., 2017; Barnes and Rautenbach, 2020; Halsne et al., 2024). The reduction in wave period/length is called “wave straining” by Holthuijsen and Tolman (1991). It is a combination of the “concertina effect” (Ardhuin et al., 2017; Wang and Sheng, 2018), related to changes in wavelength, and “energy bunching” (Baschek, 2005). Wave peak direction has values near 200°, and tidal variability is virtually absent between 20 and 23 October 2021 (not shown).

Hourly maps of the ocean current field reveal tidally driven variability in significant wave height at the Banks Peninsula station (Fig. 10). Large significant wave height from buoy observations (2.09 m) and the NZWAVE-HR forecast (2.19 m) coincides with opposing southward currents at noon on 21 October 2021 (Fig. 10a). The swell decays throughout the rest of the day but also oscillates with the currents as shown in Fig. 9. Observed and predicted significant wave height reached local minima of 1.90 and 1.89 m at 21 h on 21 October 2021, when currents are the weakest (Fig. 10j). Measured and simulated significant wave height increase again with opposing currents at 22 h on 21 October 2024 and increase to local maxima at 23 h of the same day. Wave peak direction does not show a tidal signal in either observations or NZWAVE-HR. This can be explained by the small change in relative vorticity near the wave buoy throughout the tidal cycle (red and blue shade in Fig. 10). According to Dysthe (2001), changes in relative vorticity would be responsible for changes in the wave direction or curvature of the wave ray.

Figure 10 Hourly maps of ocean currents in m s⁻¹ (black arrows) and relative vorticity (s⁻¹, red and blue shading) from NZTIDE. Observed (magenta arrow) and NZWAVE-HR-predicted (blue arrow) significant wave height and peak direction at different times on 21 October 2021 at Banks Peninsula. The values of observed and predicted (NZWAVE-HR) significant wave height for each time frame are written for reference. Model wave variables are extracted from the grid point closest to the observations.

Tidal forcing plays a larger role at Baring Head due to its location in a region of constricted circulation with strong tidal currents through Cook Strait and the measurement site being located near the coastline. Large variability is observed in peak wave direction and significant wave height (Fig. 11). Around midnight and noon on 13 January 2021, southeastward (counter) currents (red rectangles in Fig. 11) and its associated relative vorticity increase significant wave height (by around 10 cm) and steer their direction to arrive from 140°. The opposite happens when the currents flow northwestward and the waves decrease their height and shift direction to 220°, generating tidally driven oscillations in the wave direction with an amplitude of 40°. NZWAVE-HR simulates this tidal variability but with smaller amplitudes of significant wave height (5 cm) and peak wave direction (6°). NZWAVE-HR-NOTIDES, however, does not show any 12 h variability in significant wave height or peak direction. Despite the smaller tidal amplitudes in wave height and direction compared to observations, they explain the slight improvements in forecasted significant wave height and mean period (Tm01) at Baring Head (Fig. 6 and Table 2)

Figure 11 Time series of significant wave height (a), peak direction (b), and depth-averaged meridional currents (c) at Baring Head from buoy measurements (green), NZWAVE-HR-NOTIDES (magenta), and NZWAVE-HR (blue). In (c), black dots represent depth averaged currents from TELEMAC. The RMSE, correlation coefficient, and bias are presented next to each experiment label in the legend. Those statistical metrics were computed for the period shown in the figure.

Tidal variability is also observed in the mean period (Tm01), which has a tidal amplitude of about 1.5 s, whereas the model has an amplitude of 0.5 s (not shown). The smaller tidal variability in the model can be explained by the smaller current speeds in RiCOM (NZTIDE+NZSURGE-HR). TELEMAC, a newly developed storm tide forecast, shows twice the current velocity generated by RiCOM at this location, which would create larger tidal variability in NZWAVE-HR (black dots in Fig. 11c). For instance, an analytical analysis conducted by Barnes and Rautenbach (2020) based on a simple model for wave refraction from zero-velocity water into a steady current derived by Johnson (1947) shows that waves with a period of 12 s being acted on by a current field of 0.8 m s⁻¹ at a 45° angle could shift another 30° due to wave refraction. This is a similar case to the observed waves at Baring Head.

Hourly maps of the velocity, relative vorticity, and observed wave fields show a large oscillation (up to 80°) in peak direction that occurs during the tidal cycle at Baring Head. At 03:00 UTC on 13 January 2021, the observed wave peak direction has its smallest value (142.8°) when local currents and their relative vorticity approach their peak (Fig. 12c). NZWAVE-HR shows small levels of steering, and the peak direction is 188.4°. When currents are flowing northwestward and relative vorticity is negative at 07:00 UTC on 13 January 2021, observed peak direction shows its daily maximum steering and tends to be orthogonal to the flow (Figs. 11 and 12g). The predicted wave peak direction switches about 5°, and waves arrive from an angle of 193°. This smaller amplitude in the tidally forced variability in peak direction is attributed to the hydrodynamic model's low resolution and/or harmonic solver, which is not able to capture spatial (relative vorticity) and temporal variability in current speed.

Figure 12 Hourly maps of ocean currents in m s⁻¹ (black arrows) and relative vorticity (s⁻¹, red and blue shading) from NZTIDE. Observed (magenta arrow) and NZWAVE-HR-predicted (blue arrow) significant wave height and peak direction at different times on 13 January 2021 at Baring Head. Observed and predicted (NZWAVE-HR) peak direction for each time frame is written for reference. Model wave variables are extracted from the grid point closest to the observations.

A harmonic analysis of wave parameters reveals a marked influence of the main semi-diurnal tidal constituent (M2) on observed and forecast significant wave height and peak direction (Table 3). The same analysis is applied to mean period (Tm01), but the amplitudes are smaller than 1 s and are not shown. At Banks Peninsula, the observed significant wave height (Hs) and peak direction show M2 amplitudes of 2.1 cm and 2.9°, which are well represented by NZWAVE-HR (2.0 cm and 3.0°). The analysed wave parameters show larger tidally driven amplitudes at Baring Head compared to Banks Peninsula. This happens due to Baring Head's geography, which generates larger tidal currents and hence a greater influence on the wave parameters. The observed M2-forced variation in significant wave height is about 5 cm and is closely simulated by NZWAVE-HR (3.8 cm). The observed peak wave direction shows a large (17.1°) M2-forced oscillation which is not reproduced by NZWAVE-HR. This might be related to NZWAVE-HR's relatively low resolution (2 km), weaker currents in the ocean forcing, and/or lack of a two-way wave–current coupled system.

Table 3 Amplitudes of M2 tidal constituent for significant wave height (Hs) and peak direction from observations/NZWAVE-HR.

4 Discussion and conclusions

We used a set of numerical wave forecasts to evaluate the importance of downscaling, source terms, and tides on the wave forecast around Aotearoa New Zealand. The focus was mainly on two sites on the southeast coast of the country, where we have long-term time series from wave buoys. At these locations, a thorough validation of the significant wave height, mean period, and peak direction was conducted, and the impact of tides on the forecast was analysed.

Wave model downscaling showed a marked impact at a coastal scale, especially at Baring Head – a constricted coastal region. A reduction of 25 % in the significant wave height forecast error was achieved by downscaling from the low-resolution global model (GLOBALWAVE, RMSE=0.36 m) to the highest-resolution model (NZWAVE-HR-NOTIDES, RMSE=0.27 m) at Baring Head. At Banks Peninsula, however, model downscaling did not show a large impact on the wave forecast. This might be related to the station's geography, which is located in an open-coast continental shelf region, and its wave conditions can be well simulated by low-resolution models. Model downscaling generates a reduction in the bias of the absolute mean period (Tm01) at the two stations. At Banks Peninsula, mean absolute bias was reduced from −0.28 s (GLOBALWAVE) to −0.10 s (NZWAVE-HR-NOTIDES). At Baring Head, a larger reduction in the bias was found when comparing the GLOBALWAVE (−0.78 s) to the NZWAVE-HR-NOTIDES (−0.35 s) mean period (Tm01). Nevertheless, model downscaling did not show a marked impact on the mean period RMSE. No marked difference in the wave peak period was found when comparing the different downscaled models.

The use of ST6 default parameters improved the significant wave height forecast on the west coast of Aotearoa New Zealand but degraded it on its east coast. Satellite comparisons show that a forecast using ST6 default parameters reduced the significant wave height RMSE to 0.02 m (0.01 m) in the southwest Pacific (Aotearoa New Zealand) region in comparison to a model forecast using calibrated ST4 parameters. It also decreased the positive bias generated by the calibrated ST4 model from 0.13 to −0.04 m in the southwest Pacific region but created a larger region of negative bias on the eastern side of Aotearoa New Zealand. This explains the degraded significant wave height forecast at Banks Peninsula and Baring Head when comparing NZWAVE-ST6 to NZWAVE. Gorman and Oliver (2018) used altimeter wave data to find the set of ST4 parameters that minimises significant wave height RMSE via an iterative process. The same process can be applied to ST6 parameters to further reduce satellite significant wave height RMSE to a point in which NZWAVE-ST6 might improve its results against in situ observations. In addition, one may try to minimise the RMSE of in situ and satellite significant wave height altogether since altimeter data have a larger observational uncertainty (0.12–0.44 m) compared to buoy measurements (<0.10 m).

Tides and storm surge forcing showed a marked impact on the wave variability at Baring Head. This is explained by the site location, which is in a large tidally constricted region – Cook Strait. A southerly swell interacting with opposing tidal currents coming from the north with positive relative vorticity generates an increase in wave height by around 10 cm, a decrease in mean period of around 1.5 s, and a shift in wave direction of around 40°. When the tidal current is flowing northward and generating a field with negative relative vorticity, the wave parameters are affected in the opposite sense. The high-resolution tidally driven model (NZWAVE-HR) shows amplitudes of 5 cm (significant wave height), 0.5 s (mean period), and 6° (peak direction). The improvement in skill is, however, small when looking at the whole-year statistics. The RMSE of significant wave height only reduced from 0.27 to 0.26 m, and a small increase in the correlation coefficient from 0.94 (NZWAVE-HR-NOTIDES) to 0.95 (NZWAVE-HR) was observed. These improvements are smaller in comparison to the decrease in error generated by the downscaling. This can be explained by the long time series used in the study, which mixes periods where tides are more and less important – spring–neap cycle. Another factor is the weaker currents generated by the RiCOM tide and storm surge models, which might underestimate the impact of currents on the wave field. At Banks Peninsula, a similar wave–current interaction process occurs, but the inclusion of tides in the simulation did not show any impact on the wave forecast skill. This might be related to the station's geography, which is located in an open-coast continental shelf region, further from the coastline, and its wave conditions were not largely affected by changes in ocean currents and relative vorticity during the tidal cycle.

The results found by our 2 km resolution tidally driven wave forecast model (NZWAVE-HR) are comparable to the 1 km resolution simulation generated by Albuquerque et al. (2021). Those authors ran a wave hindcast for Aotearoa New Zealand using SWAN and forced by the Climate Forecast System Reanalysis (CFSR, Saha et al., 2010, 2014). At Baring Head, they found significant wave height hindcast bias, correlation coefficient, and RMSE values to be −0.09, 0.87, and 0.16 m compared to −0.05, 0.95, and 0.26 m (NZWAVE-HR). The increased correlation in NZWAVE-HR might be related to the high-spatial-and-temporal-resolution atmospheric forcing as well as the additional tidal input into the forecast. Nevertheless, the larger RMSE might be attributed to the lower spatial resolution (2 km) in NZWAVE-HR compared to the 1 km grid spacing in the hindcast. At Banks Peninsula, similar results were also found when comparing NZWAVE-HR to the wave hindcast generated by Albuquerque et al. (2021). Hindcast significant wave height bias, correlation coefficient, and RMSE were 0.05, 0.88, and 0.14 m compared to −0.06, 0.92, and 0.30 m (NZWAVE-HR). This comparison suggests that model calibration using in situ data and/or increased model resolution are still needed to further reduce the model significant wave height RMSE at Banks Peninsula. Mean period comparisons between the forecast and hindcast models show similar results. At Baring Head, the mean bias, correlation coefficient, and RMSE were −0.31, 0.78, and 1.36 s, respectively, in the SWAN hindcast and −0.34, 0.75, and 1.63 s in the forecast. Improved statistics from both models were found for Banks Peninsula in comparison to Baring Head. At Banks Peninsula, mean period bias, correlation coefficient, and RMSE were 0.09 s, 0.73 and 0.89, in the hindcast and −0.11, 0.89, and 0.65 s in the forecast (NZWAVE-HR).

Recent developments in WW3 include the implementation of the spherical multiple-cell (SMC) grid system (Li, 2022). This approach allows for adjustable spatial resolutions within a single model domain, offering the potential for more precise simulations in coastal regions. It also reduces the model resolution in less complex regions, which saves computational time. Therefore, continued validation efforts are crucial as they help identify these specific regions where high-resolution grids are the most beneficial, thereby enhancing our forecasting capabilities. Moreover, investigations towards a two-way coupled wave–current forecasting system (e.g. Couvelard et al., 2019; Fragkou et al., 2023) should be conducted.