Reliable estimates of sea level return levels are crucial for coastal flooding risk assessments and for coastal flood defence design. We describe a novel method for estimating extreme sea levels that is the first to capture seasonality, interannual variations and longer term changes. We use a joint probabilities method, with skew surge and peak tide as two sea level components. The tidal regime is predictable but skew surges are stochastic. We present a statistical model for skew surges, where the main body of the distribution is modelled empirically whilst a non-stationary generalised Pareto distribution (GPD) is used for the upper tail. We capture within-year seasonality by introducing a daily covariate to the GPD model and allowing the distribution of peak tides to change over months and years. Skew surge-peak tide dependence is accounted for via a tidal covariate in the GPD model and we adjust for skew surge temporal dependence through the subasymptotic extremal index. We incorporate spatial prior information in our GPD model to reduce the uncertainty associated with the highest return level estimates. Our results are an improvement on current return level estimates, with previous methods typically underestimating. We illustrate our method at four UK tide gauges.

Keywords: Extreme sea levels; Generalised Pareto distribution; Joint probabilities method; Non-stationarity; Skew surge

§ INTRODUCTION

Extreme sea levels pose an increasing risk to coastline communities. In the absence of any mitigation, the impacts can be severe: fatality, infrastructure damage and habitat destruction. Estimates of sea level return levels are fundamental for various purposes, including coastal flood defence design and flood risk assessments. A return level is the value we expect the annual maximum sea level to exceed with probability $p$. For a stationary series, this corresponds to a value exceeded once every $1/p$ years, on average. We are particularly interested in rare events, where $p\in[10^{-4},10^{-1}]$ to cover sea levels that are important to a range of industries affected by coastal flooding. For example, we obtain data from the Heysham tide gauge; Heysham is a coastal town in north-west England that is home to two nuclear power stations. Nuclear regulators require accurate return level estimates for $p=10^{-4}$. We also consider a further three gauges at Lowestoft, Newlyn and Sheerness; Sheerness is particularly important since it is located on the river Thames estuary where extreme sea levels can propagate down the river towards London. The annual return level, or equivalently its exceedance probability, changes with year if the series is non-stationary and discussions of appropriate design levels involve the intended lifespan of the defence <cit.>. We focus on a fixed $p$ and recognise that the return level varies across years.

Since the UK is regularly subject to coastal flooding, estimating extreme sea levels is crucial. The 1953 North Sea flood is the worst on record for the 20th century in the UK. Coastal defences were breached in 1,200 places leading to the evacuation of 30,000 people, damage to 24,000 properties and a death toll of 307 in England alone. The damage was estimated at £1.2 billion in 2014 <cit.>. Following this event, coastal flood defences were upgraded around most of the UK. Recently, there has been growing concern regarding anthropogenic sea level rise due to climate change. Rises in the mean sea level, coupled with changes in storm frequency and size, can increase the likelihood of coastal flooding. Therefore, it is increasingly important to accurately estimate sea level return levels so that coastline communities are protected against events such as that in 1953. The latest best estimates for extreme sea levels can be found in <cit.>.

<cit.> give a comprehensive overview of sea level processes. Sea levels are a combination of mean sea level, tide, surge and waves. We consider still water level (with waves filtered out), for simplicity we refer to this as sea level. The mean sea level trend is removed so that tide and surge are the only components. Tides are the regular and predictable changes in sea levels driven astronomically; these are well understood and perfectly forecast <cit.>. Surges define any departure from the predicted tidal regime, often resulting from meteorological forces such as storms and are hence sometimes called storm surges or non-tidal residual (as in Figure <ref>); these are stochastic. However, surges can also include gauge recording errors and tidal prediction errors, but these errors are typically negligible at well maintained gauges, such as those on the UK National Tide Gauge Network (see Section <ref>). Surges can also be influenced by tide-surge interaction, this is the change in the distribution of surge that is dependent on the tidal level, because the surge is essentially a wave which is influenced by the water depth. This interaction is more prominent in shallow water areas, with the largest surges typically occurring mid-tide on the rising tide and the smallest at high tide <cit.>.

Early methods to estimate extreme sea levels modelled the observed sea levels directly, ignoring the known tidal component. More recent approaches focus on the convolution of the surge and tide distributions. Due to the complex dependence between surge and tide, we instead consider skew surge and peak tide; skew surge is the difference between the maximum observed sea level and the maximum predicted tide (peak tide) within a tidal cycle, regardless of their timing (see Figure <ref>). There is a single skew surge value associated with each peak tide value (every 12 hours 26 minutes), as opposed to hourly or 15 minute interval observations of surge and tide. Even though there are fewer skew surge and peak tide observations, these components are preferred because they have much weaker dependence; <cit.> demonstrate that it is reasonable to assume these are independent at most UK sites.

Sea levels during the passage of a surge for a single tidal cycle  <cit.>.

We build on the skew surge joint probabilities method (SSJPM) of <cit.>, which assumes peak tide and skew surge are independent and that they are stationary processes within and across years, with extreme skew surges modelled by a generalised Pareto distribution (GPD). In our approach, we account for seasonality in the skew surges by adding a daily covariate to the rate of threshold exceedance and scale parameter of the GPD model. We also introduce a tidal covariate to capture skew surge-peak dependence. Skew surges also exhibit temporal dependence, we account for this using the extremal index <cit.> but in a different way to previous analyses. Since tides are deterministic, we choose our tidal samples so that they account for monthly and interannual variations. We estimate return levels by deriving distributions for the annual maxima sea levels. We obtain the estimates using our proposed method and for the SSJPM, as well as some intermediate methods, at four sites on the UK coastline, finding that the SSJPM tends to underestimate return levels.

It is fundamental to recognise the uncertainty associated with return level estimates of the highest order to add value to the point estimates. We construct confidence intervals for our estimates using a stationary bootstrap procedure <cit.>. This preserves the realism of the sea level processes. Our approach is the first to use a bootstrap procedure for uncertainty quantification on extreme sea level estimates. We add a prior distribution to the GPD shape parameter based on spatial information <cit.>; this significantly reduces the uncertainty associated with our return level estimates.

We discuss the relevant extreme value theory and existing methodology for extreme sea level estimation in Section <ref>. In Section <ref> we introduce the data and explore the seasonality of each component as well as the dependence between them. Section <ref> describes the methodology for deriving the annual maxima distribution, starting with an idealised solution derived under simplified assumptions that are then relaxed. Our return level estimates are compared with previous methods in Section <ref>. Additional supporting material is presented in <cit.>; when cross referencing, the figure, table or equation number is preceded by `S'.

§ BACKGROUND TO METHOD

In the UK, the tidal regime is semi-diurnal because a full tidal cycle of one high tide and one low tide takes approximately 12 hours 25 minutes. The timing and magnitude of low and high tides varies across the UK coastline. The range of semi-diurnal tides typically follow a 14 day cycle, where the greatest ranges are observed when the Sun and Moon are in line with the Earth at full or new Moon, this is known as spring tide. Periods of smaller tidal ranges are called neap tides, when the sun and moon are perpendicular to the Earth. Interannual variations in tides result from modulations in the semi-diurnal tidal ranges. Longer-term tidal cycles are also driven by astronomically, including the 18.6 year lunar nodal cycle and the 8.85 year cycle of lunar perigee.

Storm surges are short-term sea level changes resulting from the weather. They are caused by low atmospheric pressure and strong winds forcing high waters towards the coastline. Storm surges can only be forecast with the same time horizon as the weather. Due to the complex dependence structure between tide and storm surge, we focus on skew surge as an alternative measure of the meteorological forcing on sea levels (see Figure <ref>). <cit.> demonstrate that is it reasonable to assume skew surge and peak tide are independent.

§.§ Relevant Extreme Value Methods

When deriving a model for extremes, it is natural to first consider the maximum $M_n$ of a sequence of independent and identically distributed (iid) continuous random variables $Z_1,\ldots,Z_n$, i.e., $M_n=\max\{Z_1,\ldots,Z_n\}$. This sequence has marginal distribution function $F$ and upper end point $z^F$. If there exists sequences of constants $\{a_n>0\}$ and $\{b_n\}$, so that the rescaled block maximum $(M_n-b_n)/a_n$ has a nondegenerate limiting distribution as $n\rightarrow\infty$, then this limit must have the form

\begin{equation}
\end{equation}

where $x_+=\max\{x,0\}$ for parameters ($\mu,\sigma,\xi)\in\mathbb{R}\times\mathbb{R}_{+}\times\mathbb{R}$ representing the location, scale and shape, respectively <cit.>. This is the generalised extreme value distribution (GEV) that encompasses three families: for $\xi>0$, this corresponds to the Fréchet distribution, $\xi<0$ the Weibull and $\xi=0$ the Gumbel. Note $\xi=0$ should be interpreted as the limit as $\xi\rightarrow 0$. This result provides asymptotic motivation for using the GEV as a parametric model for observed block maxima and thus a basis to estimate and extrapolate to high return levels. However, this assumes an underlying iid process, which is unrealistic.

Now we relax the independence assumption, so that $Z_1,\ldots,Z_n$ is a stationary sequence with the same marginal distribution function $F$. This corresponds to a series whose variables may be mutually dependent, but whose statistical properties are homogeneous through time. The limiting distribution of the rescaled block maxima of a stationary process satisfying a long-range dependence condition, which ensures events long apart in time are near independent, is $G^{\theta}(z)$ with $G(z)$ as in equation (<ref>) and $\theta\in(0,1]$ the extremal index <cit.>. For an independent series $\theta=1$, but the converse is not true.

When a process exhibits extremal dependence, groups of extreme events form above high thresholds; these groups are called clusters. We define different clusters as those separated by some number of non-extreme values (i.e., below a high threshold $z$), this is called the run length $r$. Within a cluster, extreme events are considered as dependent whilst exceedances in different clusters are assumed to be independent. The extremal index tells us about clusters since it can be estimated empirically as the reciprocal of the mean cluster size of exceedances of $z$. This is known as the runs method <cit.>. Clusters are identified using an arbitrary choice of run length; this is the main pitfall with the approach. <cit.> propose the intervals estimator, based on the limiting distribution of normalised times between exceedances of $z$. This distribution is exponential for independent random variables, otherwise it is a mixture distribution of an exponential with mean $\theta^{-1}$ and a degenerate probability distribution at zero, with probabilities $\theta$ and $1-\theta$, respectively. <cit.> also propose an automatic declustering scheme using the intervals estimate. Both methods only estimate $\theta$ in the observed range of $z$. In practice the runs and intervals estimators of $\theta$ are actually estimators of the subasymptotic extremal index, for threshold level $z$ and run length $r$, defined by <cit.> as
\begin{equation}
\theta(z,r)=\mathbb{P}(\max\{ Z_2,\ldots, Z_r\}<z| Z_1>z).\label{subas_exi}
\end{equation}
Then the extremal index is the limit of expression (<ref>) as $z\rightarrow z^F$ and $r\rightarrow\infty$, with $z$ and $r$ tending to their respective limits in a related fashion.

We can also define extremes as exceedances of a high threshold $u$. If expression $\eqref{GEV_form}$ holds, then for an arbitrary term $Z$ in the sequence $Z_1,\ldots,Z_n$,
\begin{equation}
\mathbb{P}(Z>b_n+a_n z\;|\;Z>a_n+b_n u)\rightarrow H_u(z) \quad \text{ where } \quad
\end{equation}
for $z>u$ as $n\rightarrow\infty$, with $a_n$, $b_n$ as previously and $(\sigma_u,\xi)\in\mathbb{R}_{+}\times\mathbb{R}$ the scale and shape parameters, respectively <cit.>. The shape parameter is the same as that for the GEV, whilst the scale is threshold dependent since $\sigma_u=\sigma+\xi(u-\mu)$ for $\mu$ and $\sigma$ the GEV parameters. This is the generalised Pareto distribution (GPD). If $Z_1,\ldots,Z_n$ are iid, then exceedances of a high threshold $u$ are iid and have limiting GPD tail model
\begin{equation*}
\mathbb{P}(Z>z)=\lambda_u\bigg[1+\xi\bigg(\frac{z-u}{\sigma_u}\bigg)\bigg]_+^{-1/\xi}
\end{equation*}
for $z>u$ where $\lambda_u=\mathbb{P}(Z>u)$. Again, if $Z_1,\ldots,Z_n$ are stationary, a common approach is to identify clusters and decluster them to yield an approximately independent sequence of cluster maxima for which the GPD remains a valid model <cit.>.

For non-stationary processes, it is common to allow the parameters of a stationary statistical model to vary with time or another covariate. In a block maxima framework, observations in a block are assumed to be iid, so covariates in the GEV parameters cannot change within a block. In contrast, the GPD allows the covariates to vary uncontrolled over consecutive observations. A range of methods can be adopted to incorporate covariates in the model parameters $\lambda_u, \sigma_u, \xi$, including via harmonics <cit.>, splines <cit.> or generalised additive models <cit.>.

§.§ Existing Methodology

The earliest methods to estimate sea level return levels fit a GEV to the annual maxima <cit.> or the annual $r$-largest observed sea levels <cit.>. These direct approaches ignore the known tidal component even though it induces non-stationarity into the sea level series, so that the GEV limit may not be a good approximation. <cit.> show that these direct methods underestimate extreme sea levels, especially at longer return levels for tidally dominant sites.

<cit.> were the first to exploit the  decomposition of sea levels into tide and surge for extreme sea level estimation. This method is known as the joint probabilities method (JPM), where the probability distribution of extreme sea levels is derived by convolution of the distributions of these two components. The JPM gives consideration to all surge values in the data regardless of when they occurred relative to high tide. This approach forms the basis of the subsequent methods, but has some restrictive and unrealistic assumptions. Since the empirical surge distribution is used in the JPM, return level estimates are constrained by the sum of the highest predicted tide and the highest observed surge. Hourly surge observations are assumed to be independent, this is unrealistic as surge exhibits strong temporal dependence <cit.>. <cit.> do account for dependence between tide and surge by dividing the tidal range into ordered bands of equal probability and estimating the surge distribution for each band, but this gives results which are sensitive to the choice of the number of bands and their boundary levels.

The revised joint probabilities method (RJPM) of <cit.> attempts to address these limitations. For the upper tail of surges, an extreme value distribution is used to allow extrapolation beyond what has been observed, hence improving return level estimates. To account for temporal dependence, the extremal index is used. <cit.> and <cit.> model the surge-tide dependence by allowing parameters of the GEV to be functions of the tidal level. This is a difficult task because the relationship is complex <cit.>. <cit.> and <cit.> find that the RJPM is the best performing of these methods for observed and simulated data.

<cit.> propose the skew surge joint probabilities method (SSJPM) to overcome the requirement of modelling the dependence of surge on tide. Their approach assumes that skew surge and peak tide are independent. <cit.> demonstrate this is a good approximation both physically and empirically, for many sites in Europe and the USA. The SSJPM fits a GPD to the upper tail of skew surges, whilst the empirical distribution is used for the main body of the distribution. The extremal index is used to measure dependence for sea levels in adjacent tidal peaks. They find $\theta\approx{1}$ for all sites, suggesting no evidence of dependence in the upper tail. <cit.> adapt the SSJPM to account for interannual variations in the tidal regime, considering summer and winter separately and assuming peak tides and skew surge are stationary within each season.

All the methods discussed so far estimate extreme sea levels at a single site, but this gives uncertain return level estimates for sites with shorter observation periods. Spatial pooling can improve estimates at sites with limited or no data. <cit.> use regional frequency analysis to estimate extreme surges. This involves grouping statistically similar sites into homogeneous regions, then fitting an extreme value model with a constant shape parameter over all sites <cit.>. More recent approaches account for uncertainty in the region selection <cit.>. <cit.> use hindcast sea level data to interpolate SSJPM estimates along the UK coast. Other methods allow parameters of the extreme value distribution to vary with spatial covariates. <cit.> do this, allowing the location parameter of the GEV for sea level annual maxima to depend on the harmonic tidal constituents. Whereas, <cit.> spatially smooth parameters of the extreme value model for surges used in the RJPM, using predictable tidal variations along the UK coastline.

§ EXPLORATORY ANALYSIS

§.§ Data

We use data from the UK National Tide Gauge Network obtained from the British Oceanographic Data Centre; these undergo rigorous quality control before release. Sea level elevations are recorded at 44 sites along the UK coastline. We consider data from Heysham, Lowestoft, Newlyn and Sheerness only. Heysham is located on the west coast of England and has records from 1964-2016, with 17% missing. Lowestoft is on the east coast of England, with data available from 1964-2020 (4% missing). Sheerness is at the Thames Estuary, also on the east coast, with data available from 1980-2016 and 9% missing. Newlyn is located on the south coast of England and has records from 1915-2020, with 17% missing. The data are in metres above chart datum.

To compare the relative importance of skew surge and peak tide across sites, we define a surge-tide index as the observed range of skew surges divided by the range of peak tides. Heysham and Newlyn are tidally dominant with indices 0.65 and 0.70, respectively, whilst Lowestoft and Sheerness are surge dominant with indices 3.36 and 1.66, respectively. We choose to study these sites because they are typically affected by different storms and all have a long observational duration. Heysham and Lowestoft are of particular interest because of their widely different surge-tide indices. Newlyn has the longest study period in the tide gauge network. <cit.> use climate model simulations to conclude that skew surge is dependent on peak tide at Sheerness, so we also study this site. The highest astronomical tide (HAT) observed is 10.72m, 2.92m, 6.10m and 6.26m for Heysham, Lowestoft, Newlyn and Sheerness, respectively. The data were preprocessed in an attempt to remove the linear non-stationary effect of sea level rise on skew surges caused by climate change and isostatic rebound, see <cit.> for details. Therefore all results are presented relative to the mean sea level in 2017.

§.§ Seasonality and Temporal Dependence of Skew Surge and Peak Tide

Figures <ref> and <ref> show monthly boxplots of skew surge and peak tide observations at Heysham and the remaining sites, respectively. The median of the monthly skew surges remain relatively constant across months at Heysham, Newlyn and Sheerness but there is some variation at Lowestoft. The skew surge range varies across the year, exhibiting clear seasonality. Extreme skew surges and extreme peak tides typically occur in winter and at the equinoxes, respectively. Tides also exhibit interannual variability, including the 18.6 year lunar nodal cycle and the 8.85 year cycle of lunar perigee <cit.>.

Monthly boxplots of skew surge (left) and peak tide (right) at Heysham.

To further assess seasonality, we estimate the probability that a randomly selected peak tide $X$ is from month $j$ where $j=1-12$, given it is higher than some value $x_q$, the $q$th quantile of the distribution of peak tides, with $q\in[0,1]$. This is defined as
\begin{equation}
\tilde P_{X}(j;x_q)=\hat{\mathbb{P}}(m(X)=j\;|\;X>x_q),\label{prob}
\end{equation}

where $m(X)$ denotes the month of occurrence of $X$ and $\hat{\mathbb{P}}(\cdot)$ is calculated empirically. If peak tides are identically distributed over the year and all months have equal duration, then $\tilde P_{X}(j;x_q)=1/12$ for all $j$ and $q$. Of course, months vary in length, but a significant departure from $\tilde P_{X}(j;x_q)=1/12$ indicates that peak tides are not identically distributed over a year. Months with the largest peak tides will have higher values for $\tilde P_{X}(j;x_q)$ for large $q$.

Figure <ref> shows the estimates $\tilde P_{X}(j;x_q)$ for a range of $q$ at Heysham (see Figure <ref> for the remaining sites). These show that $\tilde P_{X}(j;x_{0.5})\approx1/12$ for all $j=1-12$ at Heysham, Newlyn and Sheerness. At Lowestoft this varies, agreeing with Figure <ref>. For all $q\geq 0.9$, there is clear evidence that $\tilde P_{X}(j;x_q)$ is largest in months nearest the equinoxes at all sites. At Heysham, $\tilde P_{X}(j;x_{0.999})>0$ for $j=2,3,8,9,10$; we find the same results at Newlyn but $\tilde P_{X}(1;x_{0.999})>0$ also. However, at Lowestoft and Sheerness, $\tilde P_{X}(j;x_{0.999})=0$ in months close to the spring equinox. Confidence intervals on these estimates are constructed by exploiting the property that the number of times that an event $X>x_q$ occurs in month $j$ follows a multinomial distribution, with probabilities $\tilde P_{X}(j;x_{q})$. Figures <ref> and <ref> also show 95% confidence intervals for $\tilde P_{X}(j;x_{0.999})$ at each site. These indicate that the differences in $\tilde P_{X}(j;x_{0.999})$ discussed above for each month are statistically significant.

Estimates $\tilde P_{X}(j;x_{q})$ for $j=1-12$ and $q=0.5$ (black), $0.9$ (red), 0.95 (magenta), 0.99 (purple), 0.999 (blue) at Heysham, with 95% confidence intervals for $q=0.99$.

Skew surge and peak tide both exhibit temporal dependence; we are only interested in modelling the former since peak tides are deterministic. Figure <ref> shows the autocorrelation function (acf) plots for skew surge at each site. The correlation is stronger at lower lags, as expected, and tends to zero at higher lags; although this doesn't always reach zero due to seasonality. At Heysham and Newlyn, the correlation is significantly higher for all lags above 1, than at Lowestoft and Sheerness. In Section <ref> of <cit.> we explore temporal dependence in the skew surges further by looking at the two mean measures of extremal dependence $\chi$ and $\bar\chi$ <cit.> at different high thresholds and lags. We find similar patterns over lags, but weaker dependence, for the extreme values.

§.§ Skew Surge and Peak Tide Dependence

<cit.> examine the tide-surge interaction, both empirically and physically, for sites in the southern North Sea, finding that the most extreme surges occur on rising tides. This complex dependence structure motivates the use of skew surge and peak tide in the SSJPM <cit.>. <cit.> conclude that the assumption of skew surge-peak tide independence is broadly well supported empirically and based on valid physical reasoning. However, <cit.> and <cit.> identify a weak correlation at Sheerness. <cit.> use oceanographic numerical model simulations to show that the highest skew surges tend to occur on lower peak tides at Sheerness. These findings raise questions about the validity of the independence assumption in the SSJPM. Ignoring the observed dependence would result in overestimates of return levels. We use a range of exploratory data analyses to investigate these claims in Section <ref> of <cit.>. We find that it is reasonable to assume skew surge-peak tide independence at Heysham, Lowestoft and Newlyn, but not at Sheerness. Our results suggest that this dependence structure is changing throughout the year, with the strongest dependence found in May.

§ NOVEL METHODOLOGY

§.§ Introduction

We begin by deriving the distribution of both the sea level monthly maxima $M_{j}$ for $j=1-12$ and annual maxima $M$ under several simplifying assumptions that we subsequently relax. We use the fact that the peak sea level $Z_i$ in tidal cycle $i$ can be written as the sum of skew surge $Y_i$ and peak tide $X_i$ in that tidal cycle, for all $i=1,\ldots,T$ where $T$ is the total number of tidal cycles. In Section <ref> we develop a simplified model, assuming skew surges are iid and independent of peak tide. We relax these assumptions to develop a novel method that reflects the realism of the sea level process; accounting for skew surge seasonality (Section <ref>), skew surge-peak tide dependence (Section <ref>) and temporal skew surge dependence (Section <ref>). Unless stated otherwise, inference is conducted in a likelihood framework. We provide 95% confidence intervals for parameter estimates based on the Hessian; this assumes extreme skew surges are independent, which is reasonable for model selection but in Section <ref> we use a bootstrap procedure for uncertainty quantification.

Extreme sea levels up to the $\sim20$ year return period can occur with various combinations of skew surge and peak tides, e.g., a typical skew surge value combined with the largest peak tide, through to an extreme skew surge with a typical peak tide. Since peak tides are bounded above by HAT, return periods over 20 years can only be achieved with extreme skew surges, with the very largest return periods requiring skew surges bigger than already observed. We are interested in return levels corresponding to return periods of 1 year and above, so we require an estimate for the distribution of all possible skew surges, though we are particularly interested in modelling the upper tail for inference at high return periods. Thus in each variant of our method, we develop a model for the whole skew surge distribution under the associated assumptions about the skew surge process.

§.§ Idealised Case

F_Y(y)= F̃_Y(y)        if y≤u

1-λ_u[1+ξ(y-u/σ_u)]^-1/ξ     if y> u
where $u$ is a high threshold, $\tilde F_Y$ the empirical distribution, $\lambda_u=1-\tilde{F}_Y(u)$ and $\sigma_u,\xi$ parameters of the GPD as in equation (<ref>).

Assuming that skew surge is independent of peak tide, we have that for all $i$
\begin{equation}
\mathbb{P}(Z_i\leq z)=\mathbb{P}(Y_i\leq z-X_i)=F_{Y}(z-X_i),
\end{equation}

for $-\infty<z<\infty$. To overcome sampling variation of the tides, we average over a number of tidal cycle samples $S$ and denote each sample $s$ for $s=1,\ldots,S$. Assuming all months have an equal number of tidal cycles $T$, we write the distribution of $M_j$ as

\begin{equation}
\mathbb{P}(M_j\leq z)=\frac{1}{S}\sum\limits_{s=1}^{S}\bigg(\prod\limits_{i=1}^{T}\mathbb{P}(Y_i\leq z- X_{i}^{(s)})\bigg)
=\frac{1}{S}\sum\limits_{s=1}^{S}\bigg(\prod\limits_{i=1}^{T}F_{Y}(z- X_{i}^{(s)})\bigg),\label{monmax_noseasonality}
\end{equation}

where $\{X_i^{(s)}\}_{i=1}^{T}$ is a sequence of consecutive peak tides. Since peak tides are not temporally independent of each other, the sample must be contiguous but this sequence can start anywhere in the time series, under the assumption of peak tide stationarity. The annual maxima distribution is given by,

\begin{equation}
\mathbb{P}(M\leq z)  =\prod\limits_{j=1}^{12}\bigg[\frac{1}{S}\sum\limits_{s=1}^{S}\bigg(\prod\limits_{i=1}^{T}F_{Y}(z- X_{i}^{(s)})\bigg)\bigg].\label{annmax_noseasonality}
\end{equation}
We fit this model to our data, with $\lambda_u=0.05$, and derive return level estimates in Section <ref>.

Under the assumption that the skew surges are identically distributed, we estimate the distribution of skew surges below a threshold $u$ using the empirical distribution $\tilde F_Y$. This is adequate because tide gauges on the UK National Tide Gauge Network typically have long observation lengths (>20 years) so the empirical distribution describes the main body of the data well. To enable extrapolation in the tail we use the GPD model (<ref>), with constant parameters. Then, our model is

\begin{align}
F_{Y}(y)=\begin{cases} \tilde{F}_Y(y) \quad\quad &\text{if }y\leq u\\ \\ 1-\lambda_u\big[1+\xi\big(\frac{y-u}{\sigma_u}\big)\big]_+^{-1/\xi}  &\text{if }y> u \end{cases}\label{ss_nonstat_model}
\end{align}

where $\lambda_u=1-\tilde{F}_Y(u)$ and $\sigma_u,\xi$ are the parameters of the GPD. We take $\lambda_u=0.05$. To simplify notation, we subsequently drop the $u$ subscript on the scale $\sigma$ and rate $\lambda$ parameters.

We also assume that skew surge is independent of peak tide, and as peak tides are deterministic, knowing the tidal cycle $i$ determines the peak tide $X_i$, so we have
\begin{equation*}
\mathbb{P}(Z_i\leq z)=\mathbb{P}(X_i+Y_i\leq z)
=\mathbb{P}(Y_i\leq z-X_i)=F_{Y}(z-X_i), \mbox{ for }-\infty<z<\infty.
\end{equation*}

Let $T_j$ denote the number of tidal cycles in month $j$. We use sequential monthly peak tide samples $\{X_{j_i}; j=1-12,i=1,\ldots,T_j\}$ where $j_i$ denotes the $i$th peak tide in month $j$. Then, if skew surges are independent and peak tides repeat exactly on an annual cycle, the distribution of month $j$ maximum sea level, $M_j$, is

\begin{equation}
\mathbb{P}(M_j\leq z)=\prod\limits_{i=1}^{T_j}\mathbb{P}(Y_i\leq z- X_{j_i})
=\prod\limits_{i=1}^{T_j}F_{Y}(z- X_{j_i}).
\label{monmax_noseasonality}
\end{equation}

Then the annual maxima skew surge distribution is given by,

\begin{equation}
\mathbb{P}(M\leq z)  = \prod\limits_{j=1}^{12 } \prod\limits_{i=1}^{T_j}F_{Y}(z- X_{j_i}).
\label{annmax_noseasonality}
\end{equation}

Since peak tides are temporally dependent, the tidal samples $\{X_{j_i}\}$ used in expressions (<ref>) and (<ref>) must be from contiguous peak tides.

Although cycles of periods up to a year dominate, peak tides have longer term periodicities. Previous methods, such as <cit.>, assumed peak tides are stationary within and across years so that the annual maxima distribution is given by
\begin{equation}
P(M\leq z)=\bigg(\prod\limits_{i=1}^{T} F_Y(z-X_{i})\bigg)^{1/K},\label{eqn::RJPM_ss}
\end{equation}
for $T$ the total number of observations, $K$ the number of years of observation and $F_Y$ the stationary skew surge distribution (<ref>). We incorporate interannual peak tide variations by taking the average of the yearly patterns over $K$ years, where $K\geq 19$ to address all nodal cycle variations. We denote peak tide on the $i$th tidal cycle in month $j$ of year $k$ by $X_{j_i}^{(k)}$ for $j=1-12,\;i=1,\ldots,T_j^{(k)}$ and $k=1,\ldots,K$. Here $T_j^{(k)}$ is the number of tidal cycles in month $j$ and year $k$, which varies over $k$ because the timing of the first cycle in the month can vary annually. Then the monthly and annual maxima sea level distributions are

\begin{equation}
\mathbb{P}(M_j\leq z)=\frac{1}{K}\sum_{k=1}^K\prod\limits_{i=1}^{T_j^{(k)}}F_{Y}(z- X^{(k)}_{j_i}),
\mbox{ and }
\mathbb{P}(M\leq z)  =
\frac{1}{K}\sum\limits_{k=1}^{K}\prod\limits_{j=1}^{12}\prod\limits_{i=1}^{T_j^{(k)}}F_{Y}(z- X_{j_i}^{(k)}),
\label{aver_max}
\end{equation}

respectively. This no longer has the property that the distribution of the annual maxima is the product of the monthly maxima distribution, as the monthly maxima are now dependent due to associations in peak tides in different months across years.

§.§ Skew Surge Seasonality

As discussed in Section <ref>, skew surges are not identically distributed across a year. Here we describe how we capture seasonality in three ways: for values above the threshold, below the threshold and in the exceedance probability. We define extreme skew surges using a month-specific threshold $u_j$, defined as a quantile of each monthly distribution. We use the 0.95 monthly empirical quantiles. In Section <ref> we discuss the merits of this approach relative to taking a threshold that varies on a daily scale.

Firstly we look at exceedances of the monthly thresholds and model these using the GPD as the asymptotic justification, given in Section <ref>, is still likely to hold at a monthly level. To account for seasonality we allow the parameters to change over time in a periodic fashion, as discussed in Section <ref>. We consider four models to describe how these parameters change with time, expressed in terms of a daily or monthly covariate. As a first approach to account for seasonality in extremes, <cit.> suggested allowing each month to have a separate distribution. To do so we fit a GPD with a monthly covariate on the scale $\sigma_j$ and shape $\xi_j$ parameters; we refer to this as Model $S0$. We use this basic approach for comparison only. To obtain a more parsimonious model, we fix the shape parameter to be the same value over all months but keep a monthly covariate on the scale parameter; we call this Model $S1$. These models give discontinuities at transitions between months, they assume skew surges are identically distributed within a month and have a large number of parameters. Therefore we introduce a covariate $d=1-365$ that denotes the day in year to capture within-year variations smoothly. For notational simplicity, we assume $d$ is entirely defined using the tidal cycle index $i$ and month $j$, so that there exists a function $h$ where $d=h(j,i)$. We consider two harmonic parameterisations of the scale parameter with a daily covariate. Model $S2$ uses a single harmonic defined by
\begin{equation}
\sigma_d=\alpha_\sigma+\beta_\sigma\sin\bigg(\frac{2\pi}{f}(d-\phi_\sigma)\bigg),\label{sinescale}
\end{equation}
for parameters $\alpha_\sigma>\beta_\sigma>0$, $\phi_\sigma\in[0,365)$ and periodicity $f=365$. Model $S3$ uses two harmonics with periodicities of $f$ and $f/2$. The shape parameter is equal across months for both of these models, so that Models $S2$ and $S3$ have 4 and 6 parameters, respectively.

Scale (left) and shape (right) parameter estimates for Model $S2$ (blue) and Model $S0$ (black) at Heysham. 95% confidence intervals are added to Model $S0$ parameter estimates (shown by black error bars) and to the shape parameter estimate of Model $S2$ (blue dashed lines).

We now fit the models outlined above to skew surge threshold exceedances at each site. Figures <ref> and <ref> show Model $S0$ parameter estimates for Heysham and the remaining sites, respectively. Since the 95% confidence intervals for the monthly shape parameters have considerable overlap for all sites, it is reasonable to have a common value across months as in Models $S1$, $S2$ and $S3$. Table <ref> reports the AIC and BIC scores for each model at all sites. These scores suggest that restricting the shape parameter to be an unknown constant, over months, improves model fit compared to Model $S0$ in all cases. Model $S2$ is selected as our final model for extreme skew surges as the evidence across all sites shows this is reasonable. Table <ref> gives the parameter estimates for Model $S2$ at each site. Notice that the shape parameter estimates have 95% confidence intervals that encompass zero at each site. Figures <ref> and <ref> show the scale and shape parameter estimates graphically; $\hat\sigma_d$ is higher in the winter and lower in the summer, as expected. The estimated scale parameters of Model $S2$, when averaged over each month, do not differ significantly from the estimates of Model $S0$ but do capture a smooth transition within each month.

We also considered using a log link function for the scale parameter instead of the identity link function for Models $S2$ and $S3$ but found that the identity link yields a better model fit. Furthermore the identity link function has the advantage of preserving the threshold stability property when covariates are included <cit.>.

2c     2cHeysham     2cLowestoft     2cNewlyn     2cSheerness

Model     No. Parameters     AIC     BIC     AIC     BIC     AIC     BIC     AIC     BIC

$S1$     13     -66.66     -119.52     -57.40     -118.72     -77.36     -143.83     -67.51     -123.52

$S2$     4      -69.5     -176.61     -60.47     -171.93     -88.89     -209.74     -76.18     -178.03

$S3$     6      -65.65     -161.98     -56.53     -156.85     -86.06     -194.82     -75.06     -166.73

$S4$     5     2.99     8.42     1.82     7.39     5.77     11.82     -181.57     3.53

AIC and BIC scores at each site for each skew surge model relative to Model $S0$ scores, except Model $S4$ which is measured relative to Model $S2$.

Another fundamental part of capturing within year seasonality is through the rate parameter $\lambda$. So far, we have assumed this is constant as the threshold has been set at the same quantile across months. Here, we add a daily covariate $d$ to capture smooth changes in $\lambda$ throughout each month, whist recognising that the average exceedance rate across a month is equal for all months. Let $V_d$ be a binary random variable representing whether a skew surge value exceeds its month-specific threshold $u_j$ or not, so that $V_d\sim\text{Bernoulli}(\lambda_d)$. Therefore, we use a logit link function $g(\cdot)$ to model $\lambda_d$ and capture daily changes using a generalised linear model (GLM). To account for within month variations, we relate $g(\lambda_d)$ to the day in month $d_j\in[1,31]$, standardised by the monthly mean day $\bar d_j$, so that $(d_j-\bar d_j)\in[-15,15]$ approximately. We parametrise the gradient using the day in year $d$ to account for different gradients across the year using a harmonic with periodicity $f=365.$ We refer to this model as Model $R0$ and this is given is given by

\begin{equation}
g(\lambda_d)=g(\lambda)+(d_j-\bar d_j)\beta_\lambda\sin\bigg(\frac{2\pi}{f}(d-\phi_\lambda)\bigg),\label{rate_day}
\end{equation}
for $\beta_\lambda>0$, $\phi_\lambda\in[0,365)$ which are parameters to be estimated and $\lambda$ the exceedance rate in a month (here $\lambda=0.05$). Fitting this model to our data demonstrates there is more variation in exceedance probabilities during spring and autumn compared to summer and winter. At all sites, the greatest range in $\lambda_d$ was $\sim 0.03$ in April and October. This agrees with our scale parameter model where the steepest gradient is found in spring and autumn (see Figures <ref> and <ref>). At all sites, the fitted model for $\lambda_d$ has negative gradient in months at the beginning of the year, so that the exceedance probability is higher earlier in the month (closer to winter), whereas the slope is positive later in the year.

Finally, we use a month-specific empirical distribution $\tilde F_Y^{(j)}$ for skew surges below the threshold. Bringing this together with the parameterisations of $\sigma_d$ and $\lambda_d$ (equations (<ref>) and (<ref>), respectively) for the upper tail model, the final full skew surge distribution $F_Y$ is dependent on month $j$ and day of the year $d$, and is given by
\begin{align}
F^{(d,j)}_{Y}(y)=\begin{cases} \tilde{F}^{(j)}_Y(y) &\text{if }y\leq u_j\\ \\ 1-\lambda_d\big[1+\xi\big(\frac{y-u_j}{\sigma_d}\big)\big]_+^{-1/\xi} &\text{if }y> u_j. \end{cases}\label{final_ss_model}
\end{align}
As $F^{(d,j)}_{Y}(y)$ changes with day $d$ and $d=h(j,i)$, it also changes for every tidal cycle $i$. Consequently the estimated monthly and annual maxima distributions of sea levels are
\begin{align}
&\mathbb{P}(M_j\leq{z})=\frac{1}{K}\sum\limits_{k=1}^{K}\prod\limits_{{i}=1}^{T_j^{(k)}}F_Y^{(d,j)}({z- X_{j_i}^{(k)}}),\label{sldist_ss_stat} \\
&\mathbb{P}(M\leq{z})=\frac{1}{K}\sum\limits_{k=1}^{K}\prod\limits_{{j}=1}^{12}\prod\limits_{{i}=1}^{T_j^{(k)}}F_Y^{(d,j)}({z- X_{j_i}^{(k)}}).\label{sldist_ss_stat2}
\end{align}

Heysham     Lowestoft     Newlyn     Sheerness
5lModel $S2$
$\alpha_\sigma$     0.14 (0.13, 0.15)     0.15 (0.14, 0.16)     0.076 (0.073, 0.080)     0.11 (0.10, 0.12)

$\beta_\sigma$     0.060 (0.050, 0.070)     0.080 (0.070, 0.090)     0.024 (0.020, 0.028)     0.052 (0.043, 0.061)

$\phi_\sigma$     271.51 (262.77, 280.23)     266.32 (260.01, 272.63)     273.58 (265.06, 282.10)     272.11 (262.66, 281.56)

$\xi$     0.002 (-0.042, 0.051)     0.024 (-0.023, 0.071)     -0.040 (-0.074, 0.006)     0.037 (-0.029, 0.10)
5lModel $S4$
$\gamma_\sigma^{(x)}$     0.002 (-0.005, 0.009)     -0.0051 (-0.040, 0.030)     0.0048 (-0.00048, 0.010)     -0.012 (-0.026, 0.0011)
5lModel $R1$
$\beta_\lambda$     0.0087 (0.0004, 0.017)     0.022 (0.015, 0.030)     0.024 (0.018, 0.030)     0.022 (0.014,0.032)

$\phi_\lambda$     155.66 (100.74, 210.59)     175.16 (155.86, 194.46)     209.50 (195.48, 223.52)     184.31 (160.94,207.69)

$\alpha_\lambda^{(x)}$     -0.13 (-0.18, -0.079)     -0.055 (-0.101, 0.009)     -0.063 (-0.099, 0.108)     -0.32 (-0.37,-0.26)

$\beta_\lambda^{(x)}$     0.14 (0.068, 0.21)     -0.016 (-0.084, 0.051)     0.061 (0.014, 0.108)     0.23 (0.14, 0.31)

$\phi_\lambda^{(x)}$     311.86 (281.78, 341.94)     359.95 (265.77, 454.15)     352.38 (299.32, 405.44)     278.54 (260.44, 293.63)

Parameter estimates for the scale parameter for Models $S2$ and $S4$, and the rate parameter for Model $R1$ with 95% confidence intervals, at each site.

So far, we have assumed within-year stationarity of tides and therefore selected sequences $\{X_{j_i}^{(s)}\}$ of consecutive peak tides starting from any point within the time series. We have also assumed the tidal regime is the same across years, but interannual variations exist. In this section, we describe how to account for these features.

To account for within-year variations, we simply choose our contiguous peak tide samples from the month of interest, denoted $\{X_{i}^{(j)}\}_{i=1}^{T_j}$. To account for interannual variations, we begin by deriving the distribution of the monthly maxima for a fixed year. We adapt our notation so that we know the year, as well as the month, of observation. We denote the year by $k=1,\ldots,K$ for $K$ total years of observation. Then the tidal sample is taken from the month and year of interest, $X_{i}^{(j,k)}$ for $i=1,\ldots,T_{j}^{(k)}$ where $T_j^{(k)}$ reflects that the number of tidal cycles changes for a fixed month across years because the timing of the first cycle in the month can vary annually. Assuming skew surges are independent and incorporating the seasonal model from equation (<ref>), we write the distribution of the monthly maxima in a given year $k$, $M_j^{(k)}$, as

\begin{equation}
\mathbb{P}(M_j^{(k)}\leq{z})=\prod\limits_{{i}=1}^{T^{(k)}_j}\mathbb{P}(Y_{i}^{(j)}\leq{z- X_{i}^{(j,k)}})=\prod\limits_{{i}=1}^{T^{(k)}_j}F_Y^{(d,j)}(z- X_{i}^{(j,k)}).
\end{equation}

Similarly, to find the distribution of the annual maxima, we first find the annual maxima distribution for a fixed year $k$, denoted $M^{(k)}$. This is given by
\begin{equation}
\mathbb{P}(M^{(k)}\leq{z})=\prod\limits_{j=1}^{12}\mathbb{P}(M_j^{(k)}\leq{z})=\prod\limits_{j=1}^{12}\Bigg(\prod\limits_{{i}=1}^{T^{(k)}_j}F_Y^{(d,j)}(z- X_{i}^{(j,k)})\Bigg)\label{annmax_seasonality_givenyr}.
\end{equation}

Then we assume each observed annual tidal regime has an equal probability of occurrence and average over all possibilities, we write the distribution an arbitrarily selected monthly maxima $M_j$ as
\begin{equation}
\mathbb{P}(M_j\leq{z})=\frac{1}{K}\sum_{k=1}^{K}\bigg(\prod\limits_{{i}=1}^{T^{(k)}_j}F_Y^{(d,j)}(z- X_{i}^{(j,k)})\bigg),\label{monmax_seasonality_overyears}
\end{equation}
and the distribution of the annual maxima $M$ as
\begin{equation}
\mathbb{P}(M\leq z)=\frac{1}{K}\sum\limits_{k=1}^{K}\Bigg[\prod\limits_{j=1}^{12}\Bigg(\prod\limits_{{i}=1}^{T^{(k)}_j}F_Y^{(d,j)}(z- X_{i}^{(j,k)})\Bigg)\Bigg]\label{annmax_seasonality_allyr}.
\end{equation}

Note that $\mathbb{P}(M\leq z)\neq \prod_{j=1}^{12}\mathbb{P}(M^{(k)}\leq{z})$. The right hand side of this inequality does not impose a link in the tides to be from the same year across months.

§.§ Skew Surge Dependence on Peak Tide

In Section <ref> we conclude that skew surge-peak tide independence is a reasonable assumption at Heysham, Lowestoft and Newlyn, but not at Sheerness. Here, we describe how we account for this dependence in our skew surge model. We do this for the upper tail by adding a tidal covariate to the GPD scale and rate parameters, but not to the shape parameter $\xi$ to avoid additional uncertainty. For values below the threshold, we use different empirical distributions of skew surges depending on their associated peak tide band, a similar approach to <cit.>.

We first consider how the threshold exceedance probability varies with peak tide. Here, we extend the GLM parametrisation of expression (<ref>) by adding a peak tidal covariate $x$ to the rate parameter $\lambda_d$. We linearly standardise peak tide via $(x-\bar x)/s_x$, where $\bar x$  is the mean and $s_x$ the standard deviation of all peak tide observations at each site. We parametrise the gradient in the same was as expression (<ref>), using the day in year $d$ and a harmonic with periodicity $f=365$ to capture smooth daily changes within a year. This captures the time varying dependence structure between skew surge and peak tide that we found in Section <ref> of <cit.>. We use the notation $\lambda_{d,x}$ and the following model, denoted $R1$,
\begin{equation}
g(\lambda_{d,x})=g(\lambda_d)+\bigg(\frac{x-\bar x}{s_x}\bigg)\bigg[\alpha_\lambda^{(x)}+\beta_\lambda^{(x)}\sin\bigg(\frac{2\pi}{f}(d-\phi_\lambda^{(x)})\bigg)\bigg], \label{eqn::ss_rate_tide}
\end{equation}

for $g(\cdot)$ the logit link function, $g(\lambda_d)$ defined by expression (<ref>) and $\alpha_\lambda^{(x)}\in\mathbb{R}$, $\beta_\lambda^{(x)}>0$ and $\phi_\lambda^{(x)}\in[0,365)$ being parameters to be estimated. The full model for $\lambda_{d,x}$ has 5 parameters.

We fit this GLM to model the exceedance probabilities at each site and give parameter estimates in Table <ref>. Figure <ref> shows these results graphically at Sheerness for March, June, September and December (see Figure <ref> for the other sites). As expected the estimated exceedance probability is lower at higher tides, this result is more significant in months where we found skew surge-peak tide dependence to be stronger (April - September). As in Section <ref>, $\lambda_{d,x}$ changes most with day in spring and autumn corresponding to the greatest range of skew surge within a month. Figure <ref> shows the gradient of the tidal covariate for each month, when averaged over days. At Heysham, Newlyn and Sheerness, the gradients are greatest near the equinoxes; at Lowestoft there is little variation across months.

We compare the Model $R1$ (<ref>) with Model $R0$ (<ref>) using AIC and BIC scores. Model $R1$ minimises the AIC score at Heysham, Newlyn and Sheerness; this is most notable at Sheerness with a reduction of 139. We obtain almost identical AIC scores at Lowestoft. BIC is also minimised by Model $R1$, at Heysham and Sheerness, but not at Lowestoft and Newlyn. We also compare Model $R1$ with the model in expression (<ref>) but with $\beta_\lambda^{(x)}=0$, so that the gradient term for tide does not vary with day; we found that AIC favoured Model $R1$ at all sites indicating that the dependence structure varies seasonally.

Estimated exceedance probability $\lambda_{d,x}$ under Model $R1$ in March, June, September and December with respect to $x$ being peak tide (in metres) and $d_j$ being day in month at Sheerness. Darker (lighter) points represent days later (earlier) in the month.

We also investigate adding a tidal covariate $x$ to the scale parameter of the GPD for extreme skew surges. We build on the existing parametrisation of Model $S2$ in equation (<ref>) to give Model $S4$,
\begin{equation}
\sigma_{d,x}=\alpha_\sigma+\beta_\sigma\sin\bigg(\frac{2\pi}{f}(d-\phi_\sigma)\bigg)+\gamma_\sigma^{(x)} x,\label{eqn::ss_scale_tide1}
\end{equation}
for $\alpha_\sigma>\beta_\sigma>0$ , $\phi_\sigma\in[0,365)$, $\gamma_\sigma^{(x)}\in\mathbb{R}$ parameters to be estimated, $f=365$ the periodicity and $d$ the day in year. Since we found evidence that the dependence of skew surge on peak tide changes within a year, we also considered parameterisations that allow the tidal effect to vary with day or month but found no improvement in fit relative to the loss of parsimony.

Table <ref> gives estimates of $\hat \gamma_\sigma^{(x)}$ at each site; this tells us how $\sigma_{d,x}$ changes with tide (see Table <ref> for all Model $S4$ parameter estimates). At Lowestoft and Sheerness $\hat\gamma_\sigma^{(x)}<0$ which agrees with our results of Section <ref> that extreme skew surges tend to occur on lower tides, although their 95% confidence intervals contain 0. To formally compare the fit of Model $S4$ with $S2$ we use AIC and BIC scores (see Table <ref>). Model $S4$ is only favoured by the AIC at Sheerness, this suggests the simpler model without a tidal covariate is preferred elsewhere. Given that our findings of skew surge-peak tide dependence were different at Sheerness in comparison to the other sites, we examined this further to check it wasn't an artifact of the data measurement or tide extraction processes. We explored the fits of Models $R1$ and $S4$ on a 483-year data set from a hydrodynamical model driven by a regional climate model <cit.>. We find that adding peak tide covariates to the scale and rate parameters improves fit, with similar estimates to the observations. See Section <ref> of <cit.> for details. Therefore, we proceed with Model $S4$ for the scale parameter at all sites.

Below the threshold, we found it sufficient to capture the dependence by splitting the empirical distribution of skew surges into three associated peak tide bands, i.e.,
\begin{equation}
\tilde F_{j,x}(y)=\tilde F_{j}^{(1)}(y) \; \text{ if } x\leq x_{0.33}^{(j)}, \text{ or } \tilde F_{j}^{(2)}(y) \; \text{ if } x_{0.33}^{(j)}<x\leq x_{0.67}^{(j)}, \text{ or }
\tilde F_{j}^{(3)}(y) \; \text{ if } x>x_{0.67}^{(j)},\label{eqn::ss_emp_tide}
\end{equation}
where $x_q^{(j)}$ is the $q$ quantile of the peak tide distribution in month $j$ and $\tilde F_{j}^{(l)}$ for $l=1,2,3$ is the empirical distribution of skew surges in month $j$ which are associated with the lowest ($l=1$), medium ($l=2$) and highest ($l=3$) band of peak tides. The choice of 3 tidal bands is somewhat arbitrary, but appears sufficient given the weak dependence on peak tide. We could have used a kernel smoother to ensure continuity across bands, as in <cit.>. Since our interest lies with the extreme values, we did not explore this further as this would not have made a practical difference to our return level estimates.

So our skew surge model, that is dependent on peak tide, is given by
\begin{align}
\tilde F_{j,x}(y) \quad &\text{if } y\leq u_j\\
\\
1-\lambda_{d,x}\big[1+\xi\big(\frac{y-u_j}{\sigma_{d,x}}\big)\big]_+^{-1/\xi} \quad &\text{if } y>u_j,
\end{cases}\label{ss_tidedep_model}
\end{align}
with $\tilde F_{j,x}(\cdot)$, $\sigma_{d,x}$ and $\lambda_{d,x}$ defined in expressions (<ref>), (<ref>) and (<ref>), respectively.

§.§ Skew Surge Temporal Dependence

So far, we have assumed that skew surges are independent. We now describe how we account for their temporal dependence across tidal cycles. As discussed in Section <ref>, temporal dependence causes clusters of events above high thresholds; if this is ignored, sea level return levels of annual maxima will be overestimated. The best measure of dependence for the extreme values of a stationary sequence is the extremal index $\theta$. <cit.> and <cit.> use the extremal index of the highly non-stationary sea level series. We model the sub-asymptotic extremal index $\theta(y,r)$ (defined by equation (<ref>)) of skew surges using varying thresholds (levels) $y$ with a fixed run length $r$. Then the distribution of the monthly and annual maxima sea levels are given by

\begin{align}
&\mathbb{P}(M_j\leq{z})=\frac{1}{K}\sum_{k=1}^{K}\prod\limits_{{i}=1}^{T^{(k)}_j}\big[F_Y^{(d,j,x)}(z- X_{j_i}^{(k)})\big]^{\hat\theta(z-X_{j_i}^{(k)},r)}\nonumber,\\
&\mathbb{P}(M\leq z)=\frac{1}{K}\sum\limits_{k=1}^{K}\prod\limits_{j=1}^{12}\prod\limits_{{i}=1}^{T^{(k)}_j}\big[F_Y^{(d,j,x)}(z- X_{j_i}^{(k)})\big]^{\hat\theta(z-X_{j_i}^{(k)},r)}\label{eqn::sldist_exi}.
\end{align}

The empirical estimates $\tilde\theta(y,r)$ of the sub-asymptotic extremal index are shown in Figures <ref> and <ref> for Heysham and the remaining sites, respectively. These demonstrate substantial variation with $y$. We are interested in $\theta(y,r)$ across all ranges of skew surge, but particularly in values greater than those observed to allow extrapolation for return level estimation.

Estimates of the subasymptotic extremal index $\theta(y,r)$ for different skew surge levels $y$ and run length $r=2$ at Heysham using the runs estimate (black) and our model estimate (red) of expression (<ref>). The dashed black line is the 0.99 skew surge quantile.

We estimate $\theta(y,r)$ in two stages. For $y\leq v$ where $v$ is a high threshold, we use the empirical runs estimate $\tilde\theta(y,r)$ as this is smooth over $y$ in this range. For computational efficiency purposes, we evaluate $\tilde\theta(y,r)$ on a regular grid of $y$ values in the range $[y_F, v]$ where $y_F$ is the minimum observed skew surge. Then we linearly interpolate these for any $y<v$ of interest. For $y>v$, the empirical estimate is increasingly variable so we adopt a parametric model for $\theta(y,r)$. We propose
\begin{equation}
\hat\theta(y,r)=\begin{cases}
\tilde\theta(y,r) \quad\quad &\text{if }y\leq v \\
\theta-[\theta-\tilde\theta(v,r)]\exp\Big[\scalebox{0.75}[1.0]{\( - \)}\frac{(y-v)}{\psi}\Big] \quad\quad &\text{if }y>v,
\end{cases}\label{eqn::exi_model}
\end{equation}
where $\psi>0$ and $\tilde\theta(v,r)\leq\theta\leq 1$ are parameters to be estimated. This parametric form ensures the estimate asymptotes to the extremal index $\theta$ and is continuous at $v$.

The parameters $\psi$ and $\theta$ are estimated using a weighted least squares approach with weight $w(y)=\sqrt{c(y,r)-1}$ where $c(y,r)$ is the number of clusters above $y$ separated by run length $r$. This gives a greater weight when there are more clusters, and a weight of zero when there is a single cluster. We choose $v$ to be the 0.99 quantile of all skew surges but this choice is subjective. The run length $r$ represents the approximate duration of a storm across a single site, we take $r=2,10,1$ and $10$ for Heysham, Lowestoft, Newlyn and Sheerness, respectively. These choices are supported by estimating the run length using the intervals estimator for each season, where we expect the stationary assumption to be reasonably justified <cit.>. At Heysham, $\hat\theta=1$ and $\hat\psi=0.40$; the estimates for the other sites are given in Section <ref> of <cit.>. Figures <ref> and <ref> show this model fit for Heysham and the remaining sites, respectively, compared with the entirely empirical estimates.

So far we have assumed that the only non-stationary features of extreme skew surges are in their mean value (addressed by removal of a linear trend in Section <ref>) and within year variations (described in Section <ref>). In this section, we investigate if there are any longer-term trends in the extreme skew surges that we should incorporate into our model. We can capture these trends through the parameters of the GPD model for extreme skew surges. We do not consider modelling the shape parameter as we have fixed this across months, and this parameter is associated with the most uncertainty. We adapt the parameterisations of $\sigma$ and $\lambda$ in our model to account for longer-term trends and check for any improvement in model fit.

Here, we adapt the skew surge model of Section <ref>, with scale parameter defined by expression (<ref>). We introduce a linear trend for each season using an indicator function, so that the scale parameter is given by
\begin{equation}
\sigma_{(d, k_n)} =a+b\sin\bigg(\frac{2\pi}{365}(d-\phi)\bigg)+\sum\limits_{s=1}^{4}\gamma_s k_n \mathbb{I}_{\{d\in S_s\}}\label{scale_lineartrend}
\end{equation}
where $S_s$ for $s=1,\ldots,4$ indicates the seasons Winter, Spring, Summer and Autumn respectively. Then each season has it's own linear trend given by $\gamma_s k_n$ where $k_n\in[-1,1]$ denotes the year on a normalised scale common across sites and $\gamma_s$ the trend parameter.

We fit the skew surge model with the updated scale parameter at each site using likelihood techniques. The estimates of $a,b,\phi$ and $\xi$ do not change significantly from those given in Table <ref>. Table <ref> gives the linear trend parameter estimates $\hat\gamma_s$. The results across sites vary, based on AIC and BIC scores reported in Table <ref>. Since the 95% confidence intervals of some parameter estimates contain zero, we also considered cases where we set $\gamma_s=0$ for some $s$, so that not every season has a linear trend. We omit the details.

Heysham     Lowestoft     Newlyn     All
$\hat\gamma_1$     0.031 (-0.022, 0.084)      -0.0048 ( -0.058, 0.048)     0.0049 (-0.009, 0.018)     0.007 (-0.005, 0.019)

$\hat\gamma_2$     0.031 (-0.0054, 0.068)     -0.0136 (-0.044, 0.017)      0.011 (0.002,0.019)     0.010 (0.002, 0.018)

$\hat\gamma_3$     -0.029 (-0.053, -0.0043)     0.018 (-0.0051, 0.042)     -0.002 (-0.005, 0.004)     -0.002 (-0.007, 0.004)

$\hat\gamma_4$     -0.026 (-0.069,0.016)     0.012 (-0.029, 0.052)     0.0089 (0.001, 0.018)     0.0075 (-0.001, 0.016)

Estimates of the seasonal linear trend parameters when fitting the skew surge model of equation $\eqref{final_ss_model}$ with scale parameter in equation $\eqref{scale_lineartrend}$ at Heysham, Lowestoft and Newlyn. The final column shows seasonal linear trend parameter estimates when fitting the same model but with the scale parameter in equation (<ref>), so that the trend is the same at all sites. 95% confidence intervals are given in parentheses.

Notice in Table <ref>, there is considerable overlap between $\hat\gamma_s$ for each $s$ across sites. Therefore, we investigate fitting a common seasonal linear trend parameter across sites. We add another covariate on the scale parameter to indicate the site $L=1,2,3$ corresponding to Heysham, Lowestoft and Newlyn respectively. We still use a separate harmonic at each site because the seasonality of the skew surges are not identical. The scale parameter is now given by
\begin{equation}
\sigma_{(d,k_n,L)} =\sum\limits_{l=1}^{3}\bigg[a_l+b_l\sin\bigg(\frac{2\pi}{365}(d-\phi_l)\bigg)\bigg]\mathbb{I}_{\{L=l\}}+\sum\limits_{s=1}^{4}\gamma_s k_n \mathbb{I}_{\{d\in S_s\}}.\label{scale_lineartrend_allsites}
\end{equation}

A common shape parameter for the GPD is used across sites, this is justified by the overlap of the confidence intervals for the separate site analysis (see Table <ref>). Comparing AIC and BIC in Table <ref>, we do not observe an improvement on model fit when this linear trend is incorporated, although the AIC scores are close.

2c|Heysham     2c|Lowestoft     2c|Newlyn     3cAll

None     All     None     All     None     All     None     All     Final
AIC     0     -4.17     0     3.50     0     -2.69     0      1.06     -4.39

BIC     0     -64.15     0     25.23     0     21.20     0      24.87     -0.69

AIC and BIC scores for the fit of the skew surge model of equation $\eqref{final_ss_model}$ with scale parameter in equation $\eqref{sinescale}$ with a linear trend for each season (All) relative to no linear trend in $\eqref{scale_lineartrend}$ (None) at Heysham, Lowestoft and Newlyn. The final column shows the scores when fitting the same model but with a linear trend common to all sites for each season (All) and a common linear trend for autumn, winter and spring, but no trend in Summer (Final).

We then investigate pooling various combinations of seasons together whilst setting other seasons to zero; this is motivated by the confidence intervals of parameter estimates overlapping or containing zero. AIC and BIC scores both suggest that the best model fit is one with a common linear trend across sites and across autumn, winter and spring; so that the summer linear trend parameter is set to zero Differences in AIC/BIC are small though. For this model, the scale parameter is defined by

\begin{equation}
\sigma_{(d,k_n,L)} =\sum\limits_{l=1}^{3}\bigg[a_l+b_l\sin\bigg(\frac{2\pi}{365}(d-\phi_l)\bigg)\bigg]\mathbb{I}_{\{L=l\}}+\gamma_s k_n \mathbb{I}_{\{d\in S_1,S_2,S_4\}}.\label{scale_lineartrend_final}
\end{equation}

So that if the day observation $d$ lies in seasons $S_1,S_2$ or $S_4$ (winter, spring or autumn), then the linear trend parameter $\gamma_s$ is used; otherwise if the day lies in summer, there is no seasonal trend parameter. Using likelihood inference we estimate this parameter as $\hat\gamma_s=0.00838$ (95% confidence intervals is [0.00299, 0.0138]), this corresponds to an annual increase of 0.171mm in extreme skew surges.

§ RESULTS

§.§ Introduction

Using the distribution of the monthly and annual maxima sea levels derived in Section <ref>, we estimate return levels by solving
\begin{equation*}
\mathbb{P}(M_j\leq z_{j,p})=1-p \quad \text{and} \quad \mathbb{P}(M\leq z_p)=1-p,
\end{equation*}
respectively for $p\in[0,1]$. We are interested in return levels up to the 10,000 year level, corresponding to annual exceedance probability $p=10^{-4}$. In the monthly case, this is the level we expect the monthly maxima $M_j$ to exceed, on average, every $1/p$ of that particular month (for example, every 10,000 Januarys). Whilst in the annual case, this is the level we expect the annual maxima to exceed every $1/p$ years, on average.

To assess the importance of each of the novel modelling steps in Section <ref>, we derive return level estimates from each stage; accounting for non-stationarity in each component, the dependence between them and skew surge temporal dependence. Each stage is detailed below in a nested list numbered (i)-(vii), so that each model below builds on the previous method, except (iv) from (iii). We compare these to the current approach used in practice and a baseline approach, both derived under simplifying and false assumptions. We subsequently refer to each model by name, given below, using italics. The notation follows from Section <ref>.

* Current: The methodology currently used in practice, where skew surges are assumed to be iid and peak tides are stationary, with annual maxima distribution (<ref>).

* Baseline: As in (i) but we recognise interannual variations in the tide by averaging over yearly tidal samples. The annual maxima distribution is given by
\begin{equation}
P(M\leq z)=\frac{1}{K}\sum_{k=1}^{K}\prod\limits_{i=1}^{T^{(k)}} F_Y(z-X_{i}^{(k)})\label{eqn::baselineF}
\end{equation}
where $X_i^{(k)}$, $i=1,\ldots,T^{(k)}$ represents an annual tidal sample for year $k=1,\ldots,K$, where we choose $K=100$ arbitrary but contiguous samples.

* Seasonal surge: As in (ii) but we account for the within-year seasonality of skew surges using the skew surge distribution (<ref>).

* Seasonal tide: As in (ii), but conversely to (iii) we account for within-year seasonality in peak tide and not in skew surge, so the skew surge model of expression (<ref>) is used in the monthly and annual maxima distributions (<ref>).

* Full seasonal: As in (ii) but we account for within-year seasonality of both components, with monthly and annual maxima distributions (<ref>) and (<ref>), respectively.

* (Skew surge-peak tide) interaction: As in (v) but with skew surge-peak tide dependence captured, so that the skew surge distribution (<ref>) is used.

* (Skew surge) temporal dependence: As in (vi) but accounting for temporal dependence in the skew surge series, with annual maxima distribution (<ref>).

We compare our model estimates to empirical estimates. These are restricted to return periods within the range of data, but are useful for checking whether the model is capturing the distributional properties of observed sea levels. The empirical estimates act as a guide to the truth and we do not expect our model to fit these perfectly since each empirical estimate is specific to a particular annual tidal regime, whereas in our approaches (ii)-(vii) we account for tidal variations by averaging over different samples. The empirical estimates are also sensitive to missing data and so can be biased. In Sections <ref> and <ref>, we make comparisons using point estimates to assess sensitivity to the model choice, in each case results in Section <ref> identify statistically significant differences between these models. Once the model choice is finalised we present measures of uncertainty and assess model fit in Section <ref>. In Section <ref> we use our model to find the probability that an extreme sea level occurs in a specific month, given it is higher than some level with a given return period of interest. See <cit.> for a comprehensive overview of our return level estimates.

§.§ Return Levels: Accounting for Non-stationarity

We investigate how accounting for skew surge and peak tide seasonality influences sea level return level estimates by comparing monthly return level estimates from the baseline, seasonal surge, seasonal tide and full seasonal models. We are mainly interested in differences between the baseline and full seasonal estimates. We do not expect the full seasonal model to match the empirical estimates since we have not accounted for skew surge-peak tide dependence and temporal dependence yet. The intermediate solutions (seasonal surge and seasonal tide) allow us to understand which components' non-stationarity is influencing the return levels most in different months.

Monthly maxima sea level return level estimates for Lowestoft in March, June, September and December (from left to right), estimated using the baseline (green), seasonal surge (purple dashed), seasonal tide (blue dashed) and full seasonal (cyan) models. Empirical estimates are shown by black points.

Annual maxima sea level return level estimates for Heysham (top left), Lowestoft (top right), Newlyn (bottom left) and Sheerness (bottom right), up to the 100 year level, estimated using the current (black dashed), baseline (green), full seasonal (cyan), interaction (red dot-dashed) and temporal dependence (orange dashed) models. Empirical estimates are shown by the black points.

Figures <ref> and <ref> show monthly return level estimates at Lowestoft and the other sites, respectively. At all sites and in all months, the empirical estimates lie closer to full seasonal model than the baseline. The most noticeable feature is the difference between the baseline and full seasonal model in June, this reaches 1.54m at the 10,000 year level at Lowestoft. In June skew surges are particularly low relative to the rest of the year, so ignoring their seasonality leads to significant overestimation. This overestimation increases with return period. On the other hand, in December when skew surges are relatively high compared to the other months, ignoring seasonality results in underestimates of return levels. Since the tidal range at Lowestoft is narrow compared to the range of skew surges, differences between seasonal tide and full seasonal are small relative to the differences between seasonal surge and full seasonal. For example, when the tidal range is largest at the autumn equinox in September, we observe a difference of 7cm at the 10,000 year level. Whereas at Heysham, where the tidal range is large relative to other sites, we observe significant differences between the two model estimates of 26cm at the 10,000 year level in September.

Figure <ref> shows annual return level estimates from the current, baseline and full sesaonal models, compared to empirical estimates at all sites, up to the 100 year level (Figure <ref> shows these up to the 10,000 year level). The current method gives similar results to the baseline approach for all return periods at all sites. This suggests assuming each year is identically distributed is not unreasonable. Even at Heysham, where year-to-year variations are large, the difference between the baseline and current method is small, with the current method giving a 1 year return level 4cm higher than the baseline. Since the full seasonal model lies closer to the empirical estimates than models (i)-(iv) at all sites, this highlights the importance of accounting for both forms of seasonality. At the 10,000 year level, the baseline gives a return level 6cm, 23cm and 58cm lower than the full seasonal method at Heysham, Lowestoft, and Sheerness, respectively, whereas it is 1cm higher at Newlyn.

§.§ Return Levels: Accounting for Dependence

Here, we build on the full seasonal model to capture skew surge-peak tide dependence (interaction) and skew surge temporal dependence (temporal dependence). We compare these with the empirical and full seasonal estimates. We present monthly maxima return levels of the interaction model to demonstrate the changing skew surge-peak tide dependence structure throughout the year, as discussed in Section <ref>. Then we compare annual maxima return level estimates from both models.

Monthly return level estimates for Sheerness in March, June, September and December, estimated using the full seasonal (cyan) and interaction (red dot-dashed) models when we account for and ignore skew surge-peak tide dependence, respectively. Empirical estimates are shown by the black points.

Figure <ref> shows monthly return level estimates in March, June, September and December for the full seasonal and interaction models, as well as the empirical estimates at Sheerness. We illustrate the results for Sheerness because this dependence is stronger compared to the other sites, as discussed in Section <ref>. Ignoring skew surge-peak tide dependence overestimates return levels compared to the full seasonal model for return periods greater than 10 years. At the 10,000 year level, these are overestimated by 15, 14, 17 and 10cm in March, June, September and December, respectively, with similar values down to the 1 year level. This is slightly surprising as in Section <ref> we found evidence that it is reasonable to assume skew surge and peak tide are independent in March, September and December. For Sheerness the results from the interaction model tend to lie closer to the empirical estimates than those from the full seasonal model (see Figure <ref>); this suggests that accounting for this dependence is important and better reflects the process properties. At the other sites, the return level estimates from the full seasonal and interaction model are almost indistinguishable, which echos that adding peak tide covariates into the skew surge model  when independence is a reasonable assumption does not alter results significantly.

As skew surge temporal dependence is not yet accounted for in the interaction model, its return level estimates are anticipated to slightly overestimate the empiricals. Figure <ref> compares annual return level curves for the temporal dependence and interaction models. These are very close for high return periods (>100 years). For lower return periods (<10 years) the temporal dependence estimates lie closer to the empiricals at all sites. This is a natural consequence of the extremal index model given by expression (<ref>) for temporal dependence, since the estimated extremal index is closer to zero for lower skew surges, corresponding to shorter return periods. But for the highest skew surges, and for those greater than observed, $\hat\theta(y,r)\rightarrow \hat\theta\approx 1$ and $y\rightarrow y^F$, so the largest return level estimates remain unchanged. Accounting for temporal skew surge dependence has the greatest influence at Heysham, where the acf estimates are higher compared to the other sites (see Figure <ref>). Here, the temporal dependence model gives a return values 10cm, 4cm and 2cm lower than the interaction model at the 1, 5 and 10 year return periods to better match the empirical estimates.

Return level estimates from the final model (temporal dependence) (orange), with the maximum and minimum year-specific return level estimates (dashed orange) and empirical estimates (black) at Heysham (left) and Lowestoft (right).

§.§ Assessment of Fit and Uncertainty

Figure <ref> shows annual maxima return level estimates from the temporal dependence model, judged our best model, at Heysham and Lowestoft (see Figure <ref> for the other sites). We also add the estimated maximum and minimum year-specific return levels, showing the effect of interannual peak tide variation, as given by (<ref>) without averaging over $K$ years. Figure <ref> shows yearly estimates at each site. The range of yearly estimates for Lowestoft is very narrow since the tidal range is small, compared to Heysham which has the largest tidal range of the four sites. Empirical estimates may deviate away from the return level point estimates if that specific value occurs in a year with particularly low or high tides compared to the average. Figures <ref> and <ref> show that the empirical estimates lie within the bounds for each site, at most return periods. <cit.> demonstrate that the 4.4 year perigean cycle has a greater effect on modulations in return levels than the nodal cycle in the UK. Our results are less clear (see Figure <ref>), with Sheerness year-specific return levels having an $\sim{20}$ year cycle.

We also assess goodness-of-fit by transforming the observed annual maxima to a uniform distribution using the probability integral transform with the distribution function of the annual maxima for their respective year. We do this for three cases; the baseline approach (<ref>), our final model (<ref>) and the year-specific final model. If the model is a fits well we expect the transformed annual maxima to be $\text{Uniform}(0,1)$ distributed. Figure <ref> shows PP plots for Sheerness, and the remaining sites are shown in Figure <ref>; a good fit is indicated by the empirical and model probabilities being equal so that the line of equality lies between the tolerance bounds. Figure <ref> demonstrates the improvement in fit across the three models. We formally test this using the Kolmogorov-Smirnov test at each site for the final model, with $p$ values 0.0066, 0.51, 0.0044 and 0.10 at Heysham, Lowestoft, Newlyn and Sheerness, respectively. At Heysham and Newlyn the $p$ value is significant at the 5% level so we cannot conclude a good fit, but it is a much better than for the baseline fits.

PP plots for the transformed annual maxima sea levels at Sheerness through the baseline (left), final (temporal dependence) (centre) and year specific final (right) distribution functions. The black line shows the line of equality, $y=x$. 95% tolerance bounds are shown by blue dashed lines.

We also assess the model fit uncertainty by estimating 95% confidence intervals on the 1, 10, 100, 1000 and 10,000 year return levels using a stationary bootstrap procedure. We first transform the skew surges using our final model with their corresponding seasonal and tidal covariates, call this series $\{U^{Y}_i\}$ for $i=1\ldots,T$ for $T$ the total number of observations, where $U^{Y}_i=F_{Y}^{(d,j,x)}(Y_i)$. This gives the basis for another approach to assess model fit. If the model is `correct', these transformed observations will be Uniform(0,1) and we check this formally using the Kolmogorov-Smirnov test over all years, and for each year. We discuss this further in Section <ref> of <cit.>. To account for temporal dependence, we sample $\{U^{Y}_i\}$ using a stationary bootstrap with block length $L\sim\text{Geometric}(1/10)$, so the mean block length corresponds to 5 days (the maximum duration of a storm), and the total length of each sample is equal to that of the original data <cit.>. Call this bootstrapped series of uniform variables $\{U^{B}_i\}.$ Then we transform the bootstrap sample back to the original scale using our final model, but preserving the original covariate information of $\{Y_i\}$, i.e., $\{X_{j_i},d,j\}$, call this series $\{Y^{B}_i\}$. Then we fit our tail model to the $\{Y^B\}$ series to re-estimate all of the parameters $\alpha_\sigma, \beta_\sigma, \phi_\sigma, \gamma_\sigma^{(x)}, \xi, \beta_{\lambda}, \phi_\lambda, \alpha_\lambda^{(x)}, \beta_\lambda^{(x)}, \phi_\lambda^{(x)}, \theta$ and $\psi$, as well as the monthly empirical distributions and thresholds. These estimates are then used to derive the annual maxima distribution of sea levels and estimate return levels.

We take 200 bootstrap samples and present 95% confidence intervals in Figures <ref> and <ref> at Sheerness and the other sites, respectively. The uncertainty associated with each return level increases with return period, as expected, since it becomes uncertain as we extrapolate. Uncertainty at the largest return periods can be attributed to the difficulty in accurately and precisely estimating the shape parameter of the GPD. This can be reduced by specifying a prior distribution on the shape parameter. <cit.> use a $\text{Normal}(0.0119, 0.0343^2)$ prior, based on site specific parameter estimates. We incorporate this prior information into our model, using a penalised likelihood framework for parameter estimation. The scale and rate parameter estimates remain unchanged, but our shape estimates now lie closer to zero with narrower confidence intervals at -0.019 (-0.021, 0.059), 0.014 (-0.025, 0.052), -0.027 (-0.058, 0.004) and 0.008 (-0.039, 0.054) for Heysham, Lowestoft, Newlyn and Sheerness, respectively, with 95% confidence intervals given in parentheses. Figure <ref> demonstrates how adding this prior information reduces uncertainty, since the distribution of $\hat\xi$ over 200 bootstrap samples is much narrower. Figure <ref> also shows the return level estimates and 95% bootstrap confidence intervals for this updated model; the return level curve is only affected at large return levels, as expected, where the estimates are less bounded. There is a dramatic reduction in the uncertainty associated with high return levels; at the 10,000 year level, adding prior information has reduced the range of confidence interval by 2.5m. Clearly, adding prior information on the shape parameter is important. However, we didn't do this sooner because it is important to allow the data to speak for itself when trying to the reflect the realism of the sea level process in other aspects of our methodology.

95% bootstrap confidence intervals on the final (temporal dependence) return level estimates at Sheerness (left) before (orange) and after (blue) adding a prior distribution to the shape parameter. Empirical estimates are shown by black points. Histograms of the shape parameter estimates and their densities (right) for these two models in their corresponding colours.

§.§ Sea Level Seasonality

It is helpful to understand when in the year the most extreme sea levels might occur, such as for coastal defence maintenance planning. We use our final model (temporal dependence) to investigate the seasonality of extreme sea levels. We are interested in the probability that a randomly selected sea level annual maxima $M$ is from month $j$ given it is equal to some level $z$. We consider $z_p$, a sea level with an associated annual exceedance probability $p\in[0,1]$, so that $z_{p}$ is derived from expression (<ref>). Then the probability of interest is defined by,
\begin{equation}
\hat P_M(j;z)=\hat{\mathbb{P}}(m(M) = j|M=z),\label{eqn::annmaxprob_model}
\end{equation}

where $m(M)$ denotes the month of occurrence of $M$ and $\hat{\mathbb{P}}(\cdot)$ is under our final model. This is a similar metric to that used for assessing peak tide seasonality (Section <ref>), however, since we have a model for the sea levels we condition on $M=z$ because we can obtain the form of the probability density function. Increasing $z$ to rare return levels uncovers the seasonal variations in extreme sea levels. For months where the most extreme sea levels occur, $\hat P_M(j;z)$ takes its largest values, with $\hat P_M(j;z)\rightarrow 0$ for months with the least extreme sea levels. We derive the form of $\hat P_M(j;z)$ in Section <ref> of <cit.>.

Estimates of $\hat P_M(j;z)$ for months $j=1-12$ at Sheerness, for $p=1$ (dark green), 0.1 (green), 0.01 (cyan), 0.001 (blue), and 0.0001 (dark blue). The dashed dark green line is the empirical estimate $\tilde P_M(j;z_1)$, 95% confidence intervals are for $\hat P_M(j;z_1)$.

Figure <ref> shows the estimates $\hat P_M(j;z)$ to demonstrate sea level seasonality of extreme levels at Sheerness, as well as the empirical estimate for sea levels equal to the 1 year return level (see Figure <ref> for the remaining sites). These empirical estimates tend to lie within the 95% confidence intervals for the corresponding model estimates, at the 1 year return level. This demonstrates that our model for extreme sea levels is capturing the seasonality well. Our model estimates allow us to extrapolate beyond the observed data. At Sheerness, the  empirical and model estimates at the 1 year return level are influenced by the tidal seasonality since they peak near the equinoxes. However, once we extrapolate to higher return levels, the sea level seasonality is almost entirely influenced by skew surge seasonality. This is the case for all sites at high return levels, since tides are bounded above by HAT but this is less than the 10 year sea level return level for all four sites. Therefore, extremely large skew surges are required to exceed the highest return levels, which are most likely to occur in winter.

§ DISCUSSION

By accounting for seasonality in skew surge and peak tide, the dependence between them, and temporal dependence in skew surges, we have developed a novel methodology for estimating extreme sea levels. Our results show a significant improvement on current methods, which ignored these features of the sea level processes and instead made several simplifying assumptions. Our model also allows us to study the seasonality of sea levels exceeding levels previously unobserved; this can be useful for coastal defence maintenance planning. The return levels estimated from our model present a more accurate representation of future extreme levels, and should be used for future coastal defence upgrades. Our methodology can be applied to all UK National Tide Gauge Network sites with a sufficient data record length, and is not limited to UK locations. We conclude by addressing some of the limitations with our model and suggest avenues for further exploration.

Our method assumes a steady state climate, since the existing mean sea level rise trend was removed from before analysis, and we have not included any longer-term trends in our model. Note the mean sea level trend needs to be added back onto the return level estimates presented in this paper when used in practice. Since there is growing concern with anthropogenic induced climate change and sea level rise, this should be accounted for. Whilst most of this is captured by mean sea level rise, it has been suggested there is a small, but different, trend in the extreme sea levels <cit.>. There are two approaches to do this; the preprocessing model of <cit.> where trends are modelled and removed before analysis, then covariates can be added to the extreme value model to explain trends in the tail of the distribution that are different from that in the body. Alternatively, covariates can be added to model parameters straight away without an attempt to remove any trends from the data <cit.>.

Spatial pooling provides a promising framework to capture longer-term trends due to climate change; single site trends are subtle but sharing information across sites could give more significant results. Spatial pooling also enables inference at locations with limited or no data, where our current singe site model would not perform well. In Section <ref> we add a prior penalty for the shape parameter in our model, based on spatial information; this is a first approach to spatial pooling and drastically reduced the uncertainty associated with longer term return levels. Table <ref> shows that the 95% confidence intervals of $\hat\xi$ for Model $S4$ (no penalty) overlap at all sites, therefore it may be reasonable to fix the shape parameter across sites, as an alternative approach for borrowing information. Fixing the shape parameter in homogeneous regions is a crucial step in regional frequency analysis, originally introduced by <cit.>. This is also an assumption in the Bayesian hierarchical modelling framework presented by <cit.> and a common assumption for spatial extreme value theory (see <cit.>). We refer the reader to <cit.>, <cit.> and <cit.> for different approaches to spatial pooling for extreme sea level estimation. Whilst spatial pooling is clearly an important aspect of extreme sea level estimation, it is also fundamental that the marginal site estimates are accurate, therefore our approach for capturing the realism of the sea level process should be adopted within a spatial framework.

We have shown that skew surges clearly exhibit within-year seasonality. We used a monthly threshold, defined as a quantile of the monthly skew surge distribution, to define extreme values and develop our non-stationary model. Using a quantile ensures there are a similar number of exceedances to model per month. This approach is similar to that of <cit.>, where we first assume stationarity within months and then build in the seasonal variation on a shorter temporal resolution through covariates in the GPD parameters. This meant that we were able to capture most of the non-stationarity, as well as skew surge-peak tide dependence, at the same stage of the modelling process. However, we could have considered a smoother threshold choice by using quantile regression <cit.> but we did not try this as our monthly threshold appeared sufficient.

§ ACKNOWLEDGMENTS

This paper is based on work completed while Eleanor D'Arcy was part of the EPSRC funded STOR-i centre for doctoral training (EP/S022252/1). Thanks to Jenny Sansom of the Environment Agency and Joanne Williams of National Oceanography Centre for providing the data.