Introduction

This page documents the workflow behind the “From Buoys to CCHaPS: Tracking the 40-year trend of extreme waves on Australia’s east coast” analysis prepared for the 4th Wind Wave Symposium in Perth 2025. It links long-term hourly wave observations from the Queensland and New South Wales coastal buoy networks with the Australian Climate Service wave hindcasts WHACS and CCHaPS to examine how extreme significant wave heights have behaved from the mid-1980s through to 2024.

The main objectives are to (i) derive consistent estimates of average recurrence intervals (ARIs) and 100-year significant wave heights from the buoy records, (ii) compare these with extremes simulated by WHACS and CCHaPS at co-located grid points, and (iii) assess whether mixed-climate behaviour associated with tropical-cyclone and non-tropical-cyclone events can be detected in both observations and model hindcasts. Extreme value analyses are carried out using peaks-over-threshold and annual-maximum approaches, fitting standard Gumbel and GEV models alongside a mixed-climate “two-Gumbel” model implemented in the {loopevd} R package.

We also explore non-stationary trends in extremes by allowing the Gumbel location parameter to vary with time, summarising results as 20-year changes in 100-year Hs and associated multiplication factors in annual exceedance probability (AEP). The code below reproduces the key figures and tables used in the symposium talk and provides a reproducible template for extending the analysis to other sites, time periods, or wave model products.

Previous extreme value analysis for NSW wave buoys is explained in a technical report by Baird. To our knowledge, no such analysis has been done for all of the QLD buoys.

Data sources

New South Wales wave buoys (MHL)

New South Wales wave data are provided by Manly Hydraulics Laboratory (MHL) on behalf of the NSW Government, through the long-term offshore Waverider buoy network and associated reporting program (for example the NSW Offshore Wave Heights dataset and annual wave-climate summaries). The analysis here uses hourly significant wave height (Hs) time series from the offshore buoys, following the configuration described in the Baird report on NSW offshore wave climate and the NSW Nearshore Wave Tool update. In line with that work, we take a consistent hourly data window starting in 1985 for the NSW sites used in this notebook.

The station metadata and tables in the accompanying presentation and report show that the underlying buoy archives are often longer than 1985, with some sites having earlier, non-hourly or intermittent records. Those earlier segments are not used in the extreme-value analysis here: we restrict to the continuous hourly records from 1985 onwards so that the NSW buoys are treated consistently in time and can be cleanly compared with the WHACS and CCHaPS hindcasts.

We rely entirely on the quality control and processing carried out by MHL and documented in their data products and reports; no additional QA/QC is applied in this notebook.

Queensland wave buoys (Queensland Government)

Queensland wave data are sourced from the offshore buoy network operated by the Queensland Government as part of the Coastal Data System – historical wave data and related wave climate annual summary reports. We use the processed hourly Hs records from the selected offshore Waverider sites, taking a consistent hourly data window starting in 1995, as documented in the Queensland Government wave-climate report referenced in the slides.

As with NSW, the metadata tables and original reports make clear that some Queensland buoys have records extending back before 1995 (for example, shorter bursts or earlier deployments at the same location). For the purposes of this analysis we only use the continuous hourly records from 1995 onwards, both to match the period used in the Queensland wave-climate assessment and to maintain a common baseline when comparing with model hindcasts.

Again, we work directly with the custodians’ processed hourly time series and do not perform any extra QA/QC beyond what is built into the official Queensland Government products.

Data used in this notebook

The code in this notebook:

  • reads the hourly Hs time series for the chosen NSW (MHL) and Queensland (QGov) buoys over the 1985–present (NSW) and 1995–present (Qld) analysis windows,
  • extracts annual maxima and peaks-over-threshold samples from these hourly series, and
  • passes those samples into the extreme-value models used in the “Buoys to CCHaPS” presentation.

All site choices, record lengths and start dates for the hourly data follow the configuration given in the Baird NSW report and the Queensland Government wave-climate reports, and are summarised in the tables in the slides and in the linked documentation above.

Model data extraction at buoy locations

For each buoy site, model time series were taken from the nearest model grid point to the buoy coordinates.

For the Coupled Coastal Hazard Prediction System (CCHaPS) Hindcast for Australia
(CSIRO Data Access Portal, collection csiro:65669), we use the native unstructured SCHISM–WWM mesh. CCHaPS is a 2D hydrodynamic–wave model run on a triangular grid with more than 1.4 million coastal nodes, and provides hourly water levels, depth-averaged currents and wave variables at each node.:contentReferenceoaicite:0
For each buoy we:

  • take the nominal buoy latitude–longitude from the MHL / Queensland Government metadata;
  • search the full CCHaPS node list and identify the nearest wet node by great-circle distance; and
  • extract the full model significant wave height (Hs) time series at that node over the buoy analysis window.

No additional spatial interpolation is applied: all CCHaPS–buoy comparisons use the raw Hs from this single nearest model node.

For the WHACS: Wave Hindcast for the Australian Climate Service
(DOI: 10.25919/83kj-kp95), we use the native WaveWatch III spherical multi-cell (SMC) grid. WHACS is a global wave hindcast running from 1979 to near-present, with coastal resolutions of order 6–7 km in shelf regions and coarser resolution offshore, and archives bulk and spectral wave parameters at sub-daily (3-hourly) time steps for extremes analysis.:contentReferenceoaicite:1
For each buoy we:

  • locate the nearest WHACS SMC grid cell centre to the buoy position;
  • ensure this grid point is ocean-covered (wet cell); and
  • extract the model Hs time series from that single nearest grid point over the same calendar period used for the buoy (1985–present for NSW, 1995–present for Queensland).

As with CCHaPS, we do not perform spatial averaging or interpolation between WHACS grid points; instead we use a simple nearest-neighbour approach so that extremes are compared between the buoy and a single, clearly defined model location. The extracted CCHaPS and WHACS Hs series are then treated in exactly the same way as the buoy series (annual maxima and peaks-over-threshold) in the subsequent extreme-value analysis.

Below is the R code explaining the analysis and results.

suppressPackageStartupMessages(require(terra))
## Warning: package 'terra' was built under R version 4.4.3
suppressPackageStartupMessages(require(loopevd))
## Warning: package 'loopevd' was built under R version 4.4.3
OZ_buoys = vect("data/OZ_buoys.shp")

QLD_am = as.data.frame( read.table("data/QLD_annMax.csv",header=TRUE,sep=","))
names(QLD_am) = gsub("\\.","",names(QLD_am))
names(QLD_am) = gsub("-","",names(QLD_am))

NSW_am = as.data.frame(read.table("data/NSW_annMax.csv",header=TRUE,sep=","))
QLD_years = QLD_am$year
NSW_years = NSW_am$year

QLD_am[QLD_years <= 1995,] = NA
NSW_am[NSW_years <= 1985,] = NA

QLD_am[QLD_years > 2024,] = NA
NSW_am[NSW_years > 2024,] = NA


usenames = c("cairns","townsville","mackay","haypoint","emupark","mooloolaba","brisbane",names(NSW_am))
usenames = usenames[usenames!= "X" & usenames!= "year"]

CCHaPS_am = as.data.frame(read.table("data/CCHaPS_annMax.csv",header=TRUE,sep=","))
WHACS_am = as.data.frame(read.table("data/WHACS_annMax.csv",header=TRUE,sep=","))

EC_buoys = vect(lapply(usenames,function(x) OZ_buoys[which(OZ_buoys$shortName == x)]))
EC_buoys = EC_buoys[order(EC_buoys$Lat,decreasing = TRUE)]
usenames = EC_buoys$shortName


EC_am2 = cbind(as.data.frame(sapply(usenames[1:7], function(x) QLD_am[[x]][QLD_years >= NSW_years[1]])),
              as.data.frame(sapply(usenames[8:14], function(x) NSW_am[[x]])))

years = NSW_am$year
years_CCHaPS = CCHaPS_am$year
years_WHACS = WHACS_am$year

CCHaPS_am2 = as.data.frame(sapply(usenames, function(x) CCHaPS_am[[x]]))
WHACS_am2 = as.data.frame(sapply(usenames, function(x) WHACS_am[[x]]))

#cbind(EC_buoys$shortName,names(EC_am2),names(CCHaPS_am2),names(WHACS_am2))

#EC_am$year = NSW_am$year
plot(OZ_buoys,ext = ext(buffer(EC_buoys,500000)),col = "darkgrey")
points(EC_buoys,col = "blue")
maps::map(add=TRUE)
text(EC_buoys,"shortName",pos=1:4)

Figure 1. Locations of the 14 long-term East Coast offshore wave buoys used in this analysis (blue points), shown relative to the wider Australian buoy network (grey) and coastline. Labelled sites run from Cairns in the north to Eden in the south and are used to compare observed extremes with CCHaPS and WHACS model hindcasts.

plot_empirical <- function (x, xns = NULL, unitz = "-", add=FALSE,...) 
{
    ndat = length(x)
    empi_AEP = -1/log((1:ndat)/(ndat + 1))
    if(!add) plot(empi_AEP, sort(x), col = 1, pch = 1, log = "x", xlab = "Average Recurrence Interval [years]", 
        ylab = unitz, ...)
    if(add) graphics::points(empi_AEP, sort(x),...)
    if (!is.null(xns)) 
        graphics::points(empi_AEP, sort(xns), col = "grey")
    graphics::grid()
}
i = 1
par(mfrow = c(5,3),las=1,mar = c(2,4.5,2,1)+0.1)
plot(1,1,col = NA,axes=FALSE,xlab = "",ylab = ""); legend("topleft",c("Buoy","CCHaPS","WHACS"),col = c("black","blue","red"),pch = 1)
for(i in 1:14){
  x = EC_am2[,i]
  x = x[!is.na(x)]
  plot_empirical(x,xlim = c(.4,100),main = paste0(EC_buoys$shortName[i]," n=",length(x)," depth=",round(EC_buoys$Depth_m[i])," [",round(EC_buoys$CCHaPSdept[i]),"m]"),ylim = c(1,9))
  try(plot_empirical(CCHaPS_am2[,i],add=TRUE,col="blue"))
  try(plot_empirical(WHACS_am2[,i],add=TRUE,col="red"))
  
}

Figure 2. Empirical return-interval plots of annual maximum significant wave height (Hs) at each buoy site, comparing observations (black) with co-located CCHaPS (blue) and WHACS (red) hindcast grid points. The x-axis shows average recurrence interval on a logarithmic scale; the y-axis shows annual maximum Hs, with panels ordered from north to south along the coast.

nlocs = ncol(EC_am2)
set.seed(2)
obs_gev = as.data.frame(t(sapply(1:nlocs,function(x) loopevd::evd_params(EC_am2[,x][!is.na(EC_am2[,x])],"fgev"))))
obs_gumbel = as.data.frame(t(sapply(1:nlocs,function(x) loopevd::evd_params(EC_am2[,x][!is.na(EC_am2[,x])],"fgumbel"))))
obs_gumbelx = as.data.frame(t(sapply(1:nlocs,function(x) loopevd::evd_params(EC_am2[,x][!is.na(EC_am2[,x])],"fgumbelx"))))
## Error in evd::fgumbelx(x, method = "BFGS", warn.inf = !silent) : 
##   observed information matrix is singular; use std.err = FALSE
## Error in evd::fgumbelx(x, method = "BFGS", warn.inf = !silent) : 
##   observed information matrix is singular; use std.err = FALSE
## Error in evd::fgumbelx(x, method = "BFGS", warn.inf = !silent) : 
##   observed information matrix is singular; use std.err = FALSE
## Error in evd::fgumbelx(x, method = "BFGS", warn.inf = !silent) : 
##   observed information matrix is singular; use std.err = FALSE
## Error in evd::fgumbelx(x, method = "Nelder-Mead", warn.inf = !silent) : 
##   observed information matrix is singular; use std.err = FALSE
## Error in evd::fgumbelx(x, method = "BFGS", warn.inf = !silent) : 
##   observed information matrix is singular; use std.err = FALSE
## Error in evd::fgumbelx(x, method = "Nelder-Mead", warn.inf = !silent) : 
##   observed information matrix is singular; use std.err = FALSE
CCHaPS_gev = as.data.frame(t(sapply(1:nlocs,function(x) loopevd::evd_params(CCHaPS_am2[,x][!is.na(CCHaPS_am2[,x])],"fgev"))))
CCHaPS_gumbel = as.data.frame(t(sapply(1:nlocs,function(x) loopevd::evd_params(CCHaPS_am2[,x][!is.na(CCHaPS_am2[,x])],"fgumbel"))))
CCHaPS_gumbelx = as.data.frame(t(sapply(1:nlocs,function(x) loopevd::evd_params(CCHaPS_am2[,x][!is.na(CCHaPS_am2[,x])],"fgumbelx"))))
## Error in evd::fgumbelx(x, method = "BFGS", warn.inf = !silent) : 
##   observed information matrix is singular; use std.err = FALSE
## Error in evd::fgumbelx(x, method = "BFGS", warn.inf = !silent) : 
##   observed information matrix is singular; use std.err = FALSE
## Error in evd::fgumbelx(x, method = "BFGS", warn.inf = !silent) : 
##   observed information matrix is singular; use std.err = FALSE
## Error in evd::fgumbelx(x, method = "BFGS", warn.inf = !silent) : 
##   observed information matrix is singular; use std.err = FALSE
## Error in evd::fgumbelx(x, method = "Nelder-Mead", warn.inf = !silent) : 
##   observed information matrix is singular; use std.err = FALSE
## Error in evd::fgumbelx(x, method = "BFGS", warn.inf = !silent) : 
##   observed information matrix is singular; use std.err = FALSE
## Error in evd::fgumbelx(x, method = "Nelder-Mead", warn.inf = !silent) : 
##   observed information matrix is singular; use std.err = FALSE
## Warning in evd::fgumbelx(x, method = "CG", warn.inf = !silent): optimization
## may not have succeeded
## Error in evd::fgumbelx(x, method = "BFGS", warn.inf = !silent) : 
##   observed information matrix is singular; use std.err = FALSE
## Error in evd::fgumbelx(x, method = "Nelder-Mead", warn.inf = !silent) : 
##   observed information matrix is singular; use std.err = FALSE
## Error in evd::fgumbelx(x, method = "BFGS", warn.inf = !silent) : 
##   observed information matrix is singular; use std.err = FALSE
## Error in evd::fgumbelx(x, method = "Nelder-Mead", warn.inf = !silent) : 
##   observed information matrix is singular; use std.err = FALSE
WHACS_gev = as.data.frame(t(sapply(1:nlocs,function(x) loopevd::evd_params(WHACS_am2[,x][!is.na(CCHaPS_am2[,x])],"fgev"))))
WHACS_gumbel = as.data.frame(t(sapply(1:nlocs,function(x) loopevd::evd_params(WHACS_am2[,x][!is.na(CCHaPS_am2[,x])],"fgumbel"))))
WHACS_gumbelx = as.data.frame(t(sapply(1:nlocs,function(x) loopevd::evd_params(WHACS_am2[,x][!is.na(CCHaPS_am2[,x])],"fgumbelx"))))
## Error in evd::fgumbelx(x, method = "BFGS", warn.inf = !silent) : 
##   observed information matrix is singular; use std.err = FALSE
obs_gev$AIC2 = obs_gev$AIC;obs_gev$AIC2[obs_gev$shape <= 0] = Inf
CCHaPS_gev$AIC2 = CCHaPS_gev$AIC;CCHaPS_gev$AIC2[CCHaPS_gev$shape <= 0] = Inf
WHACS_gev$AIC2  = WHACS_gev$AIC; WHACS_gev$AIC2[WHACS_gev$shape <= 0] = Inf


obs_ACS = apply(cbind(obs_gumbel$AIC,obs_gev$AIC2,obs_gumbelx$AIC),1,which.min)
CCHaPS_ACS = apply(cbind(CCHaPS_gumbel$AIC,CCHaPS_gev$AIC2,CCHaPS_gumbelx$AIC),1,which.min)
WHACS_ACS = apply(cbind(WHACS_gumbel$AIC,WHACS_gev$AIC2,WHACS_gumbelx$AIC),1,which.min)

obs_ss = sign(obs_gev$shape)+2
CCHaPS_ss = sign(CCHaPS_gev$shape)+2
WHACS_ss = sign(WHACS_gev$shape)+2


evds = c("Gumbel","GEV","Mixed Gumbel")
signs = c("-","","+")
data.frame(location = usenames,Buoy = paste(signs[obs_ss],evds[obs_ACS]), 
           CCHaPS = paste(signs[CCHaPS_ss],evds[CCHaPS_ACS]),
           WHACS = paste(signs[WHACS_ss],evds[WHACS_ACS]))
##        location           Buoy         CCHaPS          WHACS
## 1        cairns + Mixed Gumbel + Mixed Gumbel       + Gumbel
## 2    townsville          + GEV          + GEV          + GEV
## 3        mackay          + GEV          + GEV          + GEV
## 4      haypoint          + GEV + Mixed Gumbel + Mixed Gumbel
## 5       emupark + Mixed Gumbel          + GEV       + Gumbel
## 6    mooloolaba       - Gumbel       - Gumbel       - Gumbel
## 7      brisbane       - Gumbel       - Gumbel       - Gumbel
## 8      byronbay       - Gumbel       - Gumbel       - Gumbel
## 9  coffsharbour       - Gumbel       + Gumbel       - Gumbel
## 10   crowdyhead       - Gumbel       - Gumbel       - Gumbel
## 11       sydney       - Gumbel       - Gumbel       - Gumbel
## 12   portkembla       - Gumbel       - Gumbel       - Gumbel
## 13  batemansbay       - Gumbel       - Gumbel       - Gumbel
## 14         eden       - Gumbel       - Gumbel       - Gumbel

Table 1. Preferred stationary extreme value model for annual maximum Hs at each East Coast buoy site, for observations, CCHaPS and WHACS. Entries indicate whether a Gumbel, GEV or mixed-climate two-Gumbel model (“Mixed Gumbel”) gives the lowest AIC, and the sign of the GEV shape parameter (−, 0, +) where applicable.

plot_EVD = function(i){
  x = EC_am2[,i]
  x = x[!is.na(x)]
  plot_empirical(x,xlim = c(.4,100),main = paste0(EC_buoys$shortName[i]," n=",length(x)," depth=",round(EC_buoys$Depth_m[i])," [",round(EC_buoys$CCHaPSdept[i]),"m]"),ylim = c(1,9))
  #try(plot_empirical(CCHaPS_am2[,i],add=TRUE,col="blue"))


  mod_ARI  = -1/log((1:100)/(100 + 1))
  # convert to AEP (Annual Exceedance Probability)
  mod_AEP = 1 - exp(-1/mod_ARI)
  
  obs_gumbel_RL  = loopevd::qevd_vector(x = obs_gumbel[i, ],  p = 1-mod_AEP, evd_mod_str = "fgumbel")
  obs_gev_RL  = loopevd::qevd_vector(x = obs_gev[i, ],  p = 1-mod_AEP, evd_mod_str = "fgev")
  obs_gumbelx_RL = loopevd::qevd_vector(x = obs_gumbelx[i, ], p = 1-mod_AEP, evd_mod_str = "fgumbelx", interval = c(0,12))
  
  CCHaPS_gumbel_RL  = loopevd::qevd_vector(x = CCHaPS_gumbel[i, ],  p = 1-mod_AEP, evd_mod_str = "fgumbel")
  CCHaPS_gev_RL     = loopevd::qevd_vector(x = CCHaPS_gev[i, ],  p = 1-mod_AEP, evd_mod_str = "fgev")
  CCHaPS_gumbelx_RL = loopevd::qevd_vector(x = CCHaPS_gumbelx[i, ], p = 1-mod_AEP, evd_mod_str = "fgumbelx", interval = c(0,12))
  
  lines(mod_ARI,obs_gumbel_RL,lty=3,col="black")
  lines(mod_ARI,obs_gev_RL,lty=2,col="black")
  lines(mod_ARI,obs_gumbelx_RL,lty=1,col="black")
  
  #lines(mod_ARI,CCHaPS_gumbel_RL,lty=3,col="blue")
  #lines(mod_ARI,CCHaPS_gev_RL,lty=2,col="blue")
  #lines(mod_ARI,CCHaPS_gumbelx_RL,lty=1,col="blue")
  return(max(obs_gev_RL))
}
i = 1
par(mfrow = c(5,3),las=1,mar = c(2,4.5,2,1)+0.1)
plot(1,1,col = NA,axes=FALSE,xlab = "",ylab = ""); 
legend("topleft",c("Gumbel","GEV","Mixed Gumbel"),col = "black",lty = 3:1)

for(i in 1:14) plot_EVD(i)

Figure 3. Fitted stationary extreme value models for annual maximum Hs at each buoy site. Points show the empirical annual maxima, and black lines show Gumbel (dotted), GEV (dashed) and mixed-climate two-Gumbel (solid) return-level curves as a function of average recurrence interval, illustrating how well each model represents the upper tail at each location.

par(mfrow = c(2,2),las=1,mar = c(2,4.5,2,1)+0.1)
plot(1,1,col = NA,axes=FALSE,xlab = "",ylab = ""); 
legend("topleft",c("Gumbel","GEV","Mixed Gumbel"),col = "black",lty = 3:1)

rl100 = 1:3
for(i in c(1,3,11)) rl100[i] = plot_EVD(i)

Figure 4. Example return-level plots for three representative buoy sites (Cairns, Mackay and Sydney). As in Figure 3, points show empirical annual maxima and lines show fitted Gumbel, GEV and mixed-climate two-Gumbel curves, highlighting differences in tail behaviour between tropical (Cairns, Mackay) and mid-latitude (Sydney) locations.

obs_gumbelns = as.data.frame(t(sapply(1:nlocs,
                  function(x) loopevd::evd_params(EC_am2[,x][!is.na(EC_am2[,x])],
                  evd_mod_str = "fgumbel",
                  nsloc = data.frame(year = years[!is.na(EC_am2[,x])])))))

CCHaPS_gumbelns = as.data.frame(t(sapply(1:nlocs,
                  function(x) loopevd::evd_params(CCHaPS_am2[,x][!is.na(CCHaPS_am2[,x])],
                  evd_mod_str = "fgumbel",
                  nsloc = data.frame(year = years_CCHaPS[!is.na(CCHaPS_am2[,x])])))))

WHACS_gumbelns = as.data.frame(t(sapply(1:nlocs,
                  function(x) loopevd::evd_params(WHACS_am2[,x][!is.na(WHACS_am2[,x])],
                  evd_mod_str = "fgumbel",
                  nsloc = data.frame(year = years_WHACS[!is.na(WHACS_am2[,x])])))))

obs_nsAIC = apply(cbind(obs_gumbel$AIC,obs_gumbelns$AIC),1,which.min)
CCHaPS_nsAIC = apply(cbind(CCHaPS_gumbel$AIC,CCHaPS_gumbelns$AIC),1,which.min)
WHACS_nsAIC = apply(cbind(WHACS_gumbel$AIC,WHACS_gumbelns$AIC),1,which.min)

nsevds = c("","*")

data.frame(location = usenames,
           Buoy = paste(round(obs_gumbelns$locyear/obs_gumbelns$Scaled_year,3)*20, nsevds[obs_nsAIC]),
           CCHaPS = paste(round(CCHaPS_gumbelns$locyear/CCHaPS_gumbelns$Scaled_year,3)*20, nsevds[CCHaPS_nsAIC]),
           WHACS = paste(round(WHACS_gumbelns$locyear/WHACS_gumbelns$Scaled_year,3)*20, nsevds[WHACS_nsAIC]))
##        location    Buoy CCHaPS  WHACS
## 1        cairns    0.1   0.1 *  0.06 
## 2    townsville  -0.02   0.08  0.12 *
## 3        mackay   0.24  0.18 * 0.24 *
## 4      haypoint  -0.02   0.08   0.08 
## 5       emupark  -0.04  0.14 *   0.1 
## 6    mooloolaba   0.06  0.22 *  0.04 
## 7      brisbane  0.38 * 0.46 *  0.06 
## 8      byronbay  -0.38   0.16  -0.04 
## 9  coffsharbour   -0.2   0.04   0.04 
## 10   crowdyhead -0.38 *  0.14  0.24 *
## 11       sydney  0.48 * 0.54 * 0.54 *
## 12   portkembla   0.06  0.44 * 0.54 *
## 13  batemansbay   -0.1   0.4 *  0.3 *
## 14         eden   0.28  0.26 * 0.36 *

Table 2. Linear trend in the Gumbel location parameter for annual maximum Hs over a 20-year period at each site, for buoys, CCHaPS and WHACS. Values (in metres per 20 years) represent the change in the fitted Gumbel location between the start and end of a 20-year window; an asterisk marks sites where a non-stationary Gumbel model (time-varying location) is preferred over a stationary Gumbel based on AIC.

data.frame(location = usenames,
           Buoy = paste(round(exp(obs_gumbelns$locyear/obs_gumbelns$Scaled_year*20/obs_gumbelns$scale),2), nsevds[obs_nsAIC]),
           CCHaPS = paste(round(exp(CCHaPS_gumbelns$locyear/CCHaPS_gumbelns$Scaled_year*20/CCHaPS_gumbelns$scale),2), nsevds[CCHaPS_nsAIC]),
           WHACS = paste(round(exp(WHACS_gumbelns$locyear/WHACS_gumbelns$Scaled_year*20/WHACS_gumbelns$scale),2), nsevds[WHACS_nsAIC]))
##        location   Buoy CCHaPS  WHACS
## 1        cairns  1.37  1.62 *  1.32 
## 2    townsville  0.96   1.36  1.42 *
## 3        mackay  1.93  1.75 * 1.65 *
## 4      haypoint  0.95   1.47   1.42 
## 5       emupark   0.9   1.8 *  1.39 
## 6    mooloolaba  1.09  1.57 *  1.07 
## 7      brisbane 1.92 * 1.95 *  1.09 
## 8      byronbay  0.65    1.3   0.94 
## 9  coffsharbour  0.75   1.09   1.07 
## 10   crowdyhead 0.55 *  1.26  1.57 *
## 11       sydney 1.99 * 2.22 * 2.48 *
## 12   portkembla  1.09  1.86 * 2.48 *
## 13  batemansbay  0.88  1.75 * 1.76 *
## 14         eden  1.44  1.68 * 2.17 *

Table 3. Approximate multiplicative change in 100-year significant wave height over a 20-year period, inferred from the non-stationary Gumbel fits for buoys, CCHaPS and WHACS. Values greater than 1 indicate an increase in 100-year Hs over 20 years, while values less than 1 indicate a decrease; asterisks again denote sites where the non-stationary model is preferred.

i = 1
par(mfrow = c(5,3),las=1,mar = c(2,4.5,2,1)+0.1)
plot(1,1,col = NA,axes=FALSE,xlab = "",ylab = ""); legend("topleft",c("Buoy","CCHaPS","WHACS"),col = c("black","blue","red"),lty = 1)
yrs = 1980:2020
for(i in 1:14){
  x = EC_am2[,i]
  x = x[!is.na(x)]
  plot(years_WHACS,WHACS_am2[,i],type = "l",ylab ="", main = paste0(EC_buoys$shortName[i]," n=",length(x)," depth=",round(EC_buoys$Depth_m[i])," [",round(EC_buoys$CCHaPSdept[i]),"m]"),ylim = c(1,9),col = "red");grid()
  lines(years_CCHaPS,CCHaPS_am2[,i],col="blue")
  lines(years,EC_am2[,i],col = "black",lwd=2)
}

Figure 5. Time series of annual maximum Hs at each East Coast site from buoy observations (black), CCHaPS (blue) and WHACS (red). Each panel shows the annual maxima from 1985–2024 for New South Wales sites and from 1995–2024 for Queensland sites, illustrating how well each hindcast reproduces the year-to-year variability and clustering of extreme wave heights.

i = 1
par(mfrow = c(5,3),las=1,mar = c(2,4.5,2,1)+0.1)
plot(1,1,col = NA,axes=FALSE,xlab = "",ylab = ""); legend("topleft",c("Buoy","CCHaPS","WHACS"),col = c("black","blue","red"),lty = 1)
yrs = 1980:2020
for(i in 1:14){
  x = EC_am2[,i]
  x = x[!is.na(x)]
  plot(years_WHACS,WHACS_am2[,i],type = "l",ylab ="", main = paste0(EC_buoys$shortName[i]," n=",length(x)," depth=",round(EC_buoys$Depth_m[i])," [",round(EC_buoys$CCHaPSdept[i]),"m]"),ylim = c(1,9),col = "red");grid()
  lines(years_CCHaPS,CCHaPS_am2[,i],col="blue")
  lines(years,EC_am2[,i],col = "black",lwd=2)
  lines(yrs,obs_gumbelns$loc[i]+obs_gumbelns$locyear[i]*(yrs - obs_gumbelns$Centred_year[i])/obs_gumbelns$Scaled_year[i],lty=2)
  lines(yrs,CCHaPS_gumbelns$loc[i]+CCHaPS_gumbelns$locyear[i]*(yrs -  CCHaPS_gumbelns$Centred_year[i])/CCHaPS_gumbelns$Scaled_year[i],lty=2,col="blue")
}

Figure 6. As in Figure 5, but with dashed lines showing the fitted non-stationary Gumbel location parameter through time for the buoys (black) and CCHaPS (blue). These trend lines summarise gradual changes in the typical annual maximum Hs at each site and highlight locations where both observations and CCHaPS suggest strengthening extremes.

par(mfrow = c(2,2),las=1,mar = c(5,4,1,1))
plot(obs_gev$loc,CCHaPS_gev$loc,col = "blue",xlab = "Obs GEV location",ylab = "Mod  GEV location",asp=1);abline(a=0,b=1);grid()
points(obs_gev$loc,WHACS_gev$loc,col = "red")
plot(obs_gev$scale,CCHaPS_gev$scale,col = "blue",xlab = "Obs GEV scale",ylab = "Mod GEV scale",asp=1);abline(a=0,b=1);grid()
points(obs_gev$scale,WHACS_gev$scale,col = "red")
plot(obs_gev$shape,CCHaPS_gev$shape,col = "blue",xlab = "Obs GEV shape",ylab = "Mod GEV shape",asp=1);abline(a=0,b=1);grid()
points(obs_gev$shape,WHACS_gev$shape,col = "red")

Figure 7. Comparison of GEV parameter estimates for CCHaPS (blue) and WHACS (red) against buoy-based estimates across all sites: (top left) location, (top right) scale and (bottom left) shape parameters. Points show individual sites and the 1:1 line indicates perfect agreement, allowing visual assessment of how closely each hindcast reproduces the observed GEV extreme value structure.

Contact

If you have questions about this analysis, would like access to additional figures or code, or are interested in collaborating on related extreme wave or coastal hazard work, please get in touch:

Julian O’Grady
CSIRO Environment – Sea Level, Waves & Remote Sensing Team
Email: julian.ogrady@csiro.au