# Accounting for climate change in extreme sea level estimation

Eleanor D’Arcy∗ STOR-i Centre for Doctoral Training, Department of Mathematics
and Statistics, Lancaster University, LA1 4YR, UK
∗ Correspondence to<EMAIL_ADDRESS>Jonathan A. Tawn  STOR-i Centre
for Doctoral Training, Department of Mathematics and Statistics, Lancaster
University, LA1 4YR, UK
∗ Correspondence to<EMAIL_ADDRESS>
###### Abstract

Extreme sea level estimates are fundamental for mitigating against coastal
flooding as they provide insight for defence engineering. As the global
climate changes, rising sea levels combined with increases in storm intensity
and frequency pose an increasing risk to coastline communities. We present a
new method for estimating extreme sea levels that accounts for the effects of
climate change on extreme events that are not accounted for by mean sea level
trends. We follow a joint probabilities methodology, considering skew surge
and peak tides as the only components of sea levels. We model extreme skew
surges using a non-stationary generalised Pareto distribution (GPD) with
covariates accounting for climate change, seasonality and skew surge-peak tide
interaction. We develop methods to efficiently test for extreme skew surge
trends across different coastlines and seasons. We illustrate our methods
using data from four UK tide gauges.

Keywords: Extreme sea levels; Skew surge; Generalised Pareto distribution;
Non-stationarity; Climate change

##  1 Introduction

The UK coastline is one of the largest in Europe at approximately 8000km for
mainland Britain and is regularly subject to coastal flooding (Zsamboky et
al., , 2011). Coastal flooding is defined as a natural phenomenon where
coastal land is inundated by sea levels above the normal tidal conditions.
This has the potential to devastate coastal towns, damage infrastructure and
destroy habitats. In extreme cases, coastal flooding has led to the loss of
human life. The likelihood of coastal flooding is increasing with
anthropogenic induced climate change (see Figure 1) so it is increasingly
important to protect coastline communities or at least have a well-founded
scientific basis for the proposal for a managed retreat. Coastal flood
defences, such as a sea wall, protect against these consequences if they are
built to withstand the most extreme sea levels. However, resources are wasted
in building defences that are too conservative. Estimates of sea level return
levels provide crucial information for this design process; a return level is
the value we expect the annual maximum sea level to exceed with probability
$p$, i.e., once every $1/p$ years, on average, for a stationary series. We
estimate these levels for $p\in[10^{-4},10^{-1}]$ to cover a range of industry
interest, from agricultural preservation to nuclear power plant protection.

Coastal flooding is influenced by a combination of tide, surge and waves. We
are interested in the still water level, i.e., the sea level with waves
filtered out, but for simplicity we refer to this as sea level. Therefore,
tide and surge are the only components of sea levels that we consider. Tides
are the regular and predictable changes in sea levels driven astronomically;
these changes are well understood and perfectly forecast (Egbert and Ray, ,
2017). High tides generally occur once every 12 hours and 25 minutes, although
variations are possible. We refer to the maximum tide in this cycle as the
peak tide. Surges are stochastic, transient changes in sea levels often caused
by strong winds and low atmospheric pressure due to a storm, hence are often
referred to as storm surges. Surges are sometimes called the non-tidal
residual as they define any departure from the predicted tidal regime so can
also include gauge recording errors, tidal prediction errors and effects of
the tide-surge interaction. These are often available at hourly or 15 minute
intervals on the UK National Tide Gauge Network. We refer the reader to Pugh
and Woodworth,  (2014) for a complete overview of sea level processes.

An alternative decomposition of sea levels is to consider the maximum level in
a tidal cycle that can be written as the sum of skew surge and peak tide. Skew
surge is the difference between the maximum observed sea level and the peak
tide in a tidal cycle, regardless of their timing. In this case, we have less
data since observations are available once every tidal cycle. However, skew
surge and peak tide exhibit a much weaker dependence than surge and tide
(which has a complex dependence structure), so they are often preferred.
Williams et al.,  (2016) show that it is reasonable to assume skew surge and
peak tide are independent at most sites on the UK National Tide Gauge Network;
though there is physical evidence that this is not always true (Howard and
Williams, , 2021).

Long term changes in mean sea level have been widely studied (Woodworth and
Player, , 2003; Wahl et al., , 2013) via empirical assessments and using
hydrodynamic models linked to climate models. Typically linear models are fit
to estimate these trends. Similar statistical methods have been used for
extreme sea levels using regression of annual or monthly maximum data on
either sea levels or skew surges. Interestingly these simple exploratory
methods find no significant evidence for the trend in extreme sea levels to
differ from that for mean sea levels (Wong et al., , 2022; Woodworth et al., ,
2011; Weisse et al., , 2014). Complications to these methods are the large
interannual variability, the presence of seasonality and the inefficient usage
of extreme event data (through the use of maxima rather than all large
values). The difference between extreme and mean sea level trends is likely to
be of a smaller order than for mean sea level trend, hence they are more
difficult to estimate. Furthermore, only trends in average extreme values are
looked for, not changes in their variability over time. As a consequence,
inference for these properties at a single site is likely to be overloaded by
uncertainty, resulting in the hypothesis of identical trends in extreme and
mean sea levels not being rejected.

Figure 1: Global mean temperature anomalies from 1915 to 2020, relative to the
period 1961-1990, with associated uncertainties in red.

We propose a different approach which is integrated into the estimation of
return levels for extreme sea levels; this accounts for short term variations
in skew surges such as seasonality, uses all extreme skew surges above a high
month specific quantile, allows for the distribution of the extreme skew
surges, and enables pooling of information about the trend across sites.
Critically, we separately assess changes in the rate of which extreme skew
surges events occur and changes in the distribution (e.g., the mean) of these
extreme events once they have occurred.

The earliest methods for estimation of sea level return levels modelled the
sea level data directly, whilst the next approaches used a joint probabilities
method to consider surge and tide components. More recent approaches use skew
surge and peak tides. Section 2.3 gives an overview of the history of
methodology developments. We extend the most recent method, that of D’Arcy et
al.,  (2022), to account for the effects of climate change on extreme sea
level estimation. They model skew surges and combine this with the known peak
tide regime. They particularly focus on the tail of the skew surge
distribution, using a generalised Pareto distribution (GPD) to model
exceedances of a high threshold (Coles, , 2001). Covariates are added to this
model representing day of the year and peak tide, to account for seasonality
and skew surge-peak tide dependence, as well as capturing the temporal
dependence of extreme skew surges. Their results demonstrate a considerable
improvement on previous approaches since the realism of the sea level
processes is captured, and significant improvements in goodness-of-fit are
achieved. However, their model assumed that skew surges were identically
distributed across years, after a linear mean sea level trend was removed. If
climate change impacted the within year skew surge variance, or even its
distribution in a more subtle way than simply changing its mean value, then
the extreme sea levels relative to the mean sea level will also change.
Therefore we need a methodology that can incorporate such changes through to
the return level estimation. Here we will develop methods to account for non-
stationarity in this mean adjusted skew surge data to help quantify any
remaining non-stationarity in extreme skew surges.

In Sections 2.1 and 2.2 we introduce the data and relevant extreme value
theory, respectively. Section 2.3 reviews the exisiting methods used for
extreme sea level estimation, with a particular focus on that from D’Arcy et
al.,  (2022). In Section 2.4 we propose methods for investigating longer-term
trends in extreme skew surges, with respect to time and global mean
temperature anomaly (GMT), at a single site. We consider pooling information
about the trends in extreme value model across sites in Section 2.5, and
suggest methods for exploring pairwise extremal dependence in skew surges
across sites. We present the results for the single site and pooled methods in
Sections 3.2 and 3.3, respectively. Section 4 concludes this paper with a
summary of our findings and suggestions for future work.

##  2 Materials and Methods

###  2.1 Data

Sea level observations are taken from the UK National Tide Gauge Network
maintained by the Environment Agency. The data undergo rigorous quality
control and can be obtained from the British Oceanographic Data Centre (BODC).
This network is part of the National Tidal and Sea Level Facility where tidal
elevations are recorded at 44 sites along the UK coastline. We consider data
from Heysham, Lowestoft, Newlyn and Sheerness, located on the west, east,
south and east (at the Thames Estuary) coast of England, respectively. Newlyn
has the longest observational period from 1915-2020 (17% missing), whereas
Heysham, Lowestoft and Sheerness have data available for 1964-2016 (17%
missing), 1964-2020 (4% missing) and 1980-2016 (9% missing), respectively.
Observations are given as a height in metres above chart datum.

We chose to study these sites because they have different characteristics, are
typically affected by different storms and all have a long observational
periods. Heysham has the second largest tidal range on the network and is a
tidally dominant site, whereas Lowestoft is surge dominant. Sheerness is the
only site we study where it is unreasonable to assume skew surge and peak tide
are independent (D’Arcy et al., , 2022). The highest astronomical tide (HAT)
observed is 10.72m, 2.92m, 6.10m and 6.26m for Heysham, Lowestoft, Newlyn and
Sheerness, respectively. A linear mean sea level trend was removed from the
data at each site therefore all of our results are presented relative to the
mean sea level in the year 2017. Environment Agency,  (2018) details this
preprocessing stage; they remove estimated linear trends of 1.52, 2.27, 1.73
and 1.81mm per year from Heysham, Lowestoft, Newlyn and Sheerness,
respectively. Of course, these trends incorporate land level changes as well
as climate caused sea level changes, and also are based on different time
periods as they correspond to the sample record at each site.

###  2.2 Extreme Value Inference

Within extreme value inference, it is natural to first consider modelling the
maximum of a sequence $M_{n}=\max\\{Z_{1},\ldots,Z_{n}\\}$. We first assume
this sequence is independent and identically distributed (iid) with marginal
distribution $F$ and upper end point $x^{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 the distribution function $G$ of the maximum must
be of the form

$G(z)=\exp\bigg{\\{}-\bigg{[}1+\xi\bigg{(}\frac{z-\mu}{\sigma}\bigg{)}\bigg{]}^{-1/\xi}_{+}\bigg{\\}},$
(2.1)

where $x_{+}=\max\\{x,0\\}$ whatever the distribution function $F$. This
distributional model $G$ has three parameters $\mu\in\mathbb{R}$,
$\sigma\in\mathbb{R}_{+}$ and $\xi\in\mathbb{R}$ representing the location,
scale and shape, respectively (Coles, , 2001). This is known as the
generalised extreme value distribution (GEV). For $\xi>0$ this corresponds to
the Fréchet distribution, $\xi<0$ the Weibull and $\xi=0$ the Gumbel (although
$\xi=0$ should be interpreted as the limit as $\xi\rightarrow 0$). This
result, often referred to as the extremal types theorem, gives an asymptotic
justification to use the GEV as a model for block maxima, often taken to be
annual maxima in environmental applications. However, in these settings, an
iid assumption is usually unrealistic. A more commonly accepted assumption is
stationarity, where the series can exhibit mutual dependence, but the
statistical properties are homogeneous through time. If we now assume that
$Z_{1},\ldots,Z_{n}$ are from a stationary series that satisfies a long-range
dependence condition, so that events far enough apart in time are near
independent. Under these conditions this limiting distribution must be of the
form $G^{\theta}(z)$ with $G(z)$ in expression (2.1) and $\theta\in(0,1]$ the
extremal index (Leadbetter et al., , 1983).

If a process exhibits dependence, values above a high threshold $z$ form
clusters, for example during a storm that spans multiple days we might observe
several extreme skew surge values consecutively. We identify clusters as those
separated by a run length $r$ defined as the number of consecutive
observations below the high threshold $z$, i.e., non-extreme values. Choosing
this run length can be subjective, though Ferro and Segers,  (2003) propose an
automated selection procedure based on the distribution of all times between
consecutive exceedances of $z$. We can reasonably assume that observations in
different clusters are independent, but this is not the case for observations
in the same cluster. The extremal index $\theta$ provides information about
clusters because it can be empirically estimated (known as the runs estimate)
as the reciprocal of the mean cluster size (Smith and Weissman, , 1994). These
are both actually estimates of the subasymptotic extremal index

$\theta(z,r)=\mathbb{P}(\max\\{Z_{2},\ldots,Z_{r}\\}<z|Z_{1}>z).$ (2.2)

Then the extremal index is defined as the limit of expression (2.2) as
$z\rightarrow z^{F}$ and $r\rightarrow\infty$ in a related fashion (Ledford
and Tawn, , 2003).

An alternative, and more popular, approach to defining extreme values is as
exceedances of a high threshold $u$. If the extremal types theorem holds, then
for an arbitrary term $Z$ in the series $Z_{1},\ldots,Z_{n}$,

$\mathbb{P}(Z>b_{n}+a_{n}z\;|\;Z>a_{n}+b_{n}u)\rightarrow H_{u}(z)\quad\text{
where }\quad
H_{u}(z)=\bigg{[}1+\xi\bigg{(}\frac{z-u}{\sigma_{u}}\bigg{)}\bigg{]}_{+}^{-1/\xi}$
(2.3)

for $z>u$ as $n\rightarrow\infty$, with $\\{a_{n}>0\\}$ and $\\{b_{n}\\}$
sequences of constants and $H_{u}$ is non-degenerate. This is known as the
generalised Pareto distribution (GPD) and has two parameters
$\sigma_{u}\in\mathbb{R}_{+}$ and $\xi\in\mathbb{R}$ representing the scale
and shape, respectively (Coles, , 2001). Notice the scale parameter is
threshold dependent since $\sigma_{u}=\sigma+\xi(u-\mu)$ for $\mu$ and
$\sigma$ the GEV parameters; the shape parameter is the same as that for the
GEV. Assuming $Z_{1},\ldots,Z_{n}$ are iid, exceedances of a high threshold
$u$ are also iid and have limiting GPD tail model

$\mathbb{P}(Z>z)=\lambda_{u}\bigg{[}1+\xi\bigg{(}\frac{z-u}{\sigma_{u}}\bigg{)}\bigg{]}_{+}^{-1/\xi}$
(2.4)

for $z>u$ where $\lambda_{u}=\mathbb{P}(Z>u)$. We can write the mean of
excesses of the threshold $u$ as

$\mathbb{E}(Z-u|Z>u)=\frac{\sigma_{u}}{1-\xi}.$ (2.5)

However, if $Z_{1},\ldots,Z_{n}$ are a dependent stationary series, a common
approach is to identify clusters and decluster them (e.g., by considering
cluster maxima only) to yield an approximately independent sequence so that
the asymptotic justification for the GPD remains valid (Smith et al.,  (1997),
Fawcett and Walshaw,  (2007)). We subsequently drop the $u$ subscript on the
scale $\sigma$ and rate $\lambda$ parameters.

###  2.3 Existing Methodology

The earliest methods directly modelled sea levels, but this ignores the known
tidal component (Graff, , 1978; Coles et al., , 1999; Tawn, , 1988). Dixon and
Tawn,  (1999) demonstrate that these approaches underestimate return levels.
Pugh and Vassie,  (1978) were the first to exploit the components of sea
levels in their joint probabilities method (JPM) using surge and tide. Tawn,
(1992) presents the revised joint probabilities method (RJPM) to address
limitations associated with the JPM; mainly, they use an extreme value
distribution to model the upper tail of surges to allow extrapolation beyond
the range of observed values and, through a parameteric model, attempt to
account for tide-surge dependence. The main pitfall with both of these
approaches is that surge and tide have a complex joint distribution which is
difficult to model effectively. Batstone et al.,  (2013) proposed the skew
surge joint probabilities method (SSJPM) to avoid this complexity. This uses
skew surge and peak tide as two components of sea levels, since they have a
much weaker dependence and can be reasonably assumed to be independent at most
sites on the UK National Tide Gauge Network (Williams et al., , 2016). Baranes
et al.,  (2020) build on this by accounting for interannual tidal variations
and considering separate distributions for summer and winter skew surges; this
is the quasi non-stationary skew surge joint probabilities method (qn-SSJPM).

We build on the sea level model presented by D’Arcy et al.,  (2022) that uses
skew surge and peak tide as two components of sea levels in a joint
probabilities framework. This was first approach to capture the within year
seasonality of each component and the dependence between them by adding
covariates to the model parameters. It also accounts for skew surge temporal
dependence which addresses previous issues of overestimation at short return
periods. We describe their model for the annual maxima sea levels $M$. For a
tidal cycle $i$, the peak sea level $Z_{i}$ can be written as the sum of the
deterministic peak tide $X_{i}$ and stochastic skew surge $Y_{i}$. We first
present their skew surge model, then describe how this is combined with the
known peak tides to derive a sea level distribution. Lastly, we detail their
model for the extremal index of skew surges used to derive the annual maxima
distribution.

Since extreme sea levels can occur with various combinations of skew surge and
peak tide, it is important to have a model for the entire skew surges
distribution. To split the distribution into the body and tail, they use a
monthly threshold $u_{j}$ for $j=1,\ldots,12$ to account for seasonality, with
$u_{j}$ being a quantile, for a fixed percentile, of month $j$’s skew surge
distribution. This choice ensures a similar number of exceedances each month.
They use the 0.95 percentile, this is chosen based on monthly parameter
stability plots (Coles, , 2001). Skew surges below these thresholds are
modelled using the monthly empirical distribution $\tilde{F}_{j,x}$ to capture
within year non-stationarity, that is also dependent on tides $x$ to account
for skew-surge peak tide dependence. This empirical distribution is split into
three associated peak tide bands:

$\tilde{F}_{j,x}(y)=\begin{cases}\tilde{F}_{j}^{(1)}(y)\quad\text{if }x\leq
x_{0.33}^{(j)}\\\ \\\ \tilde{F}_{j}^{(2)}(y)\quad\text{if
}x_{0.33}^{(j)}<x\leq x_{0.67}^{(j)}\\\ \\\
\tilde{F}_{j}^{(3)}(y)\quad\text{if }x>x_{0.67}^{(j)},\end{cases}$ (2.6)

where $x_{q}^{(j)}$ denotes the $q$ quantile of the peak tide distribution for
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. Since tide
gauges on the UK National Tide Gauge Network usually have long observational
periods, this can reliably model the main body of the data. For exceedances of
the monthly threshold, they use a non-stationary GPD dependent on day in year
$d=1,\ldots,365$, month $j$ and peak tide $x$. Therefore, the full skew surge
model is given by

$\displaystyle
F_{Y}^{(d,j,x)}(y)=\begin{cases}\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}$ (2.7)

where $\lambda_{d,x}$, $\sigma_{d,x}$ and $\xi$ are parametric functions to be
estimated. Notice that the shape parameter $\xi$ does not vary with any
covariate; this is kept fixed to avoid introducing additional uncertainty
associated with estimating this parameter. The rate and scale parameters both
depend on day and peak tide. They model the scale parameter using a harmonic
for seasonal variations and a linear trend in terms of tide,

$\sigma_{d,x}=\alpha_{\sigma}+\beta_{\sigma}\sin\bigg{(}\frac{2\pi}{f}(d-\phi_{\sigma})\bigg{)}+\gamma_{\sigma}x,$
(2.8)

for $\alpha_{\sigma}>\beta_{\sigma}>0$ , $\phi_{\sigma}\in[0,365)$,
$\gamma_{\sigma}\in\mathbb{R}$ parameters to be estimated and $f=365$ the
periodicity. The rate parameter is modelled similarly, using a generalised
linear model with logit link function and a harmonic to capture seasonal
variations. They also use a harmonic to capture how skew surge-peak tide
dependence changes with time; D’Arcy et al.,  (2022) show that this
relationship varies throughout the year at Sheerness, with the strongest
dependence occurring in May. This parameterisation is given by

$g(\lambda_{d,x})=g(\lambda)+(d_{j}-\bar{d}_{j})\beta_{\lambda}^{(d)}\sin\bigg{(}\frac{2\pi}{f}(d-\phi_{\lambda}^{(d)})\bigg{)}+\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{]},$
(2.9)

for $g(\cdot)$ the logit link function, $\lambda$ the constant exceedance
probability in a month, $d_{j}\in[1,31]$ day in month (standardised by the
monthly mean day $\bar{d}_{j}$), $\bar{x}$ is the mean and $s_{x}$ the
standard deviation of peak tides, and $\alpha_{\lambda}^{(x)}\in\mathbb{R}$,
$\beta_{\lambda}^{(d)},\beta_{\lambda}^{(x)}>0$,
$\phi_{\lambda}^{(d)},\phi_{\lambda}^{(x)}\in[0,365)$ are parameters to be
estimated.

To derive a distribution for the sea levels, D’Arcy et al.,  (2022) use a
joint probabilities method and the fact that peak tides are deterministic. So
that

$\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.$ (2.10)

Let $T_{j}^{(k)}$ denote the number of tidal cycles in month $j$ and year $k$.
They capture within and across year peak tide non-stationarity by using
sequential monthly and yearly peak tide samples $X_{j_{i}}^{(k)}$, so that
$j_{i}$ denotes the $i$th peak tide in month $j$, and $k=1,\ldots,K$
represents the year. Since peak tides are temporally dependent, the samples
$\\{X_{j_{i}}^{(k)}\\}$ are from contiguous peak tides. Then the distribution
of the annual maxima sea level $M$ is

$\mathbb{P}(M\leq
z)=\frac{1}{K}\sum_{k=1}^{K}\prod\limits_{j=1}^{12}\prod\limits_{i=1}^{T_{j}^{(k)}}F_{Y}^{(d,j,x)}(z-X^{(k)}_{j_{i}})^{\theta(z-X_{j_{i}}^{(k)},r)}$
(2.11)

where $F_{Y}^{(d,j,x)}$ is the skew surge model (2.7) and
$\theta(z-X_{j_{i}}^{(k)},r)$ is a model for the extremal index, dependent on
skew surge level $y=z-X_{j_{i}}^{(k)}$ and run length $r$, to capture temporal
dependence of skew surges. This model is given by

$\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{(}-\frac{y-v}{\psi}\big{)}\quad\quad&\text{if
}y>v,\end{cases}$ (2.12)

where $v$ is a high threshold (they take the 0.99 percentile), $\psi>0$ and
$\tilde{\theta}(v,r)\leq\theta\leq 1$ are parameters to be estimated and
$\tilde{\theta}(y,r)$ is the empirical runs estimate.

###  2.4 Incorporating Interannual Variations to Skew Surge Distribution

We provide a framework to explore longer-term trends in extreme skew surges
that may result from an increase in storm frequency and intensity. After
removing the mean sea level trend, we follow the approach of Eastoe and Tawn,
(2009) by adding yearly and global mean temperature anomaly (GMT) covariates
into the scale and rate parameters to the GPD model for extreme skew surges
presented by D’Arcy et al.,  (2022). We do not consider adding covariates to
the shape parameter or to the empirical distribution used for non-extreme skew
surges. Another option would be to add covariates to the threshold choice, but
it is difficult to account for the uncertainty in the threshold selection in
extreme value inference (Wadsworth and Tawn, , 2012; Northrop et al., , 2017).
Since we are interested in the temporal change in extreme events, it seems
problematic to allow the threshold to also vary with time.

The model of D’Arcy et al.,  (2022) already accounts for short term variations
in the threshold exceedance rate and the GPD scale parameter. So, here we are
focusing on the additional long term changes in these two features, knowing
that estimates of these are not contaminated by short term features. Trends in
the two features tell us about different aspects of the occurrence of extreme
skew surge events. An increase in the threshold exceedance rate tells us
simply that more extreme events are occurring over time or with GMT increases.
In contrast, increases in the scale parameter inform us that the nature of the
extreme events are changing, in that their average size is getting larger. So
it is of interest to explore both aspects.

First we consider how the threshold exceedance probability varies with year
$k$ and GMT in year $k$ measured in ${}^{\circ}C$, denoted $m_{k}$. This tells
us how the frequency of extreme skew surges is changing in response to climate
change. We refer to the model for $\lambda_{d,x}$ introduced by D’Arcy et al.,
(2022) as $R0$, given by (2.9). We propose four model extensions of R0 to
account for how the threshold exceedance rate also changes with $k$ (Models
$R1$ and $R2$) and with $m_{k}$ (Models $R3$ and $R4$); with the odd numbered
models having a single trend across the year and the even numbered models
having a different linear trend per season, with seasons
$\\{\mathcal{S}_{s},s=1,2,3,4\\}$ denoting winter (December, January,
February), spring (March, April, May), summer (June, July, August) and autumn
(October, November December), respectively. These models are parametrised as
follows,

$\displaystyle\text{Model~{}}R1\text{: }\quad
g(\lambda_{d,x,\tilde{k}})=g(\lambda_{d,x})+\delta_{\lambda}^{(\tilde{k})}\tilde{k},$
(2.13) $\displaystyle\text{Model~{}}R2\text{: }\quad
g(\lambda_{d,x,\tilde{k}})=g(\lambda_{d,x})+\sum\limits_{s=1}^{4}\delta_{\lambda,s}^{(\tilde{k})}\tilde{k}\bm{1}_{\\{d\in\mathcal{S}_{s}\\}},$
(2.14) $\displaystyle\text{Model~{}}R3\text{: }\quad
g(\lambda_{d,x,m_{k}})=g(\lambda_{d,x})+\delta_{\lambda}^{(m)}m_{k},$ (2.15)
$\displaystyle\text{Model~{}}R4\text{: }\quad
g(\lambda_{d,x,m_{k}})=g(\lambda_{d,x})+\sum\limits_{s=1}^{4}\delta_{\lambda,s}^{(m)}m_{k}\bm{1}_{\\{d\in\mathcal{S}_{s}\\}},$
(2.16)

where $g(\cdot)$ is the logit link function, which is selected to help our
modelling of probabilities with linear models,
$\delta_{\lambda}^{(\tilde{k})},\delta_{\lambda,s}^{(\tilde{k})},\delta_{\lambda}^{(m)},\delta_{\lambda,s}^{(m)}\in\mathbb{R}$
are parameters to be estimated ($s=1,2,3,4$), and $\tilde{k}\in\mathbb{R}$
denotes the standardised year defined as $\tilde{k}=(k-1968)/53$ where $k$ is
the year of observation. Here 1968 is the midpoint of all years of observation
and 53 is half of the range. We standardise in this way so that the results
are comparable across sites with different observation lengths. For our study
period, the covariates take values $\tilde{k}\in[-1,1]$ and
$m_{k}\in(-0.56,0.94)$; recall GMT is an anomaly centred at the temperature in
the period 1961-1990, so it has been standardised somewhat already.

We consider these four models to investigate whether time or GMT is the best
linear predictor of extreme skew surge non-stationarity over our observation
period, and to explore if the longer-term trends are non-stationary within a
year, for example, if extreme skew surges are becoming more frequent in the
winter but less so in summer. For Model $R1$, we are particularly interested
in the change in threshold exceedance probability over the period 1920-2020
(100 years), this is given by
$\Delta_{\lambda}^{(\tilde{k})}=\lambda_{d,x,b}-\lambda_{d,x,a}$, for
$a=-0.91$ (1920) and $b=1$ (2020). Similarly for Model $R3$, we define the
change in exceedance probability with an increase in GMT of $1^{\circ}$C as
$\Delta_{\lambda}^{(m)}=\lambda_{d,x,1}-\lambda_{d,x,0}=\lambda_{d,x,1}-\lambda_{d,x}$.

Next, we investigate how the GPD scale parameter changes with year and GMT to
understand if the magnitude of extreme events is changing due to climate
change. We extend the $\sigma_{d,x}$ parameterisation (2.8) of D’Arcy et al.,
(2022) (call this Model $S0$) and consider four models that capture changes
with year, GMT and season as we did for the threshold exceedance rate,

$\displaystyle\text{Model~{}}S1\text{:
}\quad\sigma_{d,x,k}=\sigma_{d,x}+\delta_{\sigma}^{(\tilde{k})}\tilde{k},$
(2.17) $\displaystyle\text{Model~{}}S2\text{:
}\quad\sigma_{d,x,k}=\sigma_{d,x}+\sum\limits_{s=1}^{4}\delta_{\sigma,s}^{(\tilde{k})}\tilde{k}\bm{1}_{\\{d\in\mathcal{S}_{s}\\}},$
(2.18) $\displaystyle\text{Model~{}}S3\text{:
}\quad\sigma_{d,x,m_{k}}=\sigma_{d,x}+\delta_{\sigma}^{(m)}m_{k},$ (2.19)
$\displaystyle\text{Model~{}}S4\text{:
}\quad\sigma_{d,x,m_{k}}=\sigma_{d,x}+\sum\limits_{s=1}^{4}\delta_{\sigma,s}^{(m)}m_{k}\bm{1}_{\\{d\in\mathcal{S}_{s}\\}},$
(2.20)

with parameters
$\delta_{\sigma}^{(\tilde{k})},\delta_{\sigma,s}^{(\tilde{k})},\delta_{\sigma}^{(m)},\delta_{\sigma,s}^{(m)}\in\mathbb{R}$
to be estimated and $\tilde{k}$, $m_{k}$, $\mathcal{S}_{s}$ as in
(2.13)-(2.16).

###  2.5 Spatial Pooling

####  2.5.1 Improved Inference by Pooling

So far we have described the modelling of extreme skew surges at a single
site. However, this approach can be very inefficient, particularly for sites
with short records or where the physical processes exhibit similarly over the
sites, e.g., the same storm events effect all of the different sites. In such
cases we would anticipate certain parameters of the extreme surge skew
distribution to be similar, or even identical, in value across sites. By
imposing this feature into the inference and carrying out joint inferences
across sites, known as pooling, this can lead to large improvements in
parameter estimation, by effectively sharing information about extreme events
across sites, which in turn reduces estimation uncertainty resulting in
narrower confidence intervals.

In the model of D’Arcy et al.,  (2022) the benefits of pooling was illustrated
for the shape parameter. This parameter is known to very difficult to estimate
with much precision, with the variability in its estimator being the primary
source of uncertainty in return level estimation. This parameter has been
recognised across a wide spectrum of problems as being very similar for a
given process over large spatial regions, e.g., for rainfall (Davison et al.,
, 2012), sea levels (Dixon et al., , 1998) and air temperature (Huser and
Genton, , 2016) with different values for the shape parameter for plains and
mountains. D’Arcy et al.,  (2022) use information from Environment Agency,
(2018) that the shape parameter estimates, estimated separately from each site
over the UK, followed a normal distribution with mean $0.0119$ and variance
$0.0343$. They account for this in the likelihood inference, using this
distribution as a prior penalty function. D’Arcy et al.,  (2022) obtained
shape parameter estimates which were more similar over sites with much reduced
uncertainty, thus resulting in uncertainty reduction of high return level
estimates. For example, for the 10,000 year return level at Sheerness, the 95%
confidence interval was reduced by 2.5m, corresponding to a factor of 6.

In our context the difficult parameters to estimate are those of the longer-
term trends in expressions (2.13)-(2.16) for the threshold exceedance rate and
(2.17)-(2.20) for the GPD scale parameter. Here we also want to share spatial
information through pooling. Given that we do not know if these trend
parameters are identical over sites, and we only are illustrating the method
for four sites, we undertake a formal likelihood testing method to assess the
evidence to see if we can treat these trend parameters as constant over sites,
without reducing the quality of the fit relative to the improved parsimony.

The pooled inference procedure involves a combined likelihood function
$L({\bm{\theta}})$ which combines the likelihood functions
$L_{i}({\bm{\theta}}_{i})$ from each of the $i=1,\ldots,4$ sites through

$L({\bm{\theta}})=\prod_{i=1}^{4}L_{i}({\bm{\theta}}_{i}),$

where ${\bm{\theta}}_{i}$ are the parameters for site $i$ and
${\bm{\theta}}=({\bm{\theta}}_{1},\ldots,{\bm{\theta}}_{4})$. This likelihood
enables hypothesis testing to be carried out, to assess the evidence for
whether certain parameters are the same at all, or a subset set, of the sites,
i.e., is the time trend gradient parameter the same at all sites. The joint
likelihood function then enables the sharing of information about this common
parameter across sites whilst allowing the other parameters to vary over
sites. The choice of this joint likelihood function has potential
restrictions; since it is a product over sites, this implicitly implies that
extreme skew surges are being assumed to be independent across the sites. In
cases where this assumption is unreasonable, the point estimates of the
parameters will still be good (asymptotically consistent) but the variance of
the estimates and the confidence intervals for the parameters will be
underestimated, with the degree of underestimation dependent on the level of
ignored true dependence there is between skew surges at the different sites.
Therefore before exploiting the pooling strategy it is important to check that
the independence assumption, for the extreme values of skew surge, is a
reasonable approximation.

####  2.5.2 Spatial Independence Diagnostics

We discuss how to check for pairwise dependence between skew surges at
different sites. Kendall’s $\tau$ correlation coefficient can be used to check
for dependence in all skew surge observations; this is a measure of rank
correlation so is robust to outliers but it is a measure across all values of
the variables. However, since our interest lies with the dependence of the
extreme values, it is natural to also study the two main measures of extremal
dependence $\chi$ and $\bar{\chi}$ (Coles et al., , 1999), as described next.

Let $Y^{A}$ and $Y^{B}$ denote skew surge random variables at two different
sites $A$ and $B$, in the same tidal cycle with marginal distribution
functions $F_{A}$ and $F_{B}$ respectively. The simplest measure of dependence
is to see how the joint probability of $Y^{A}$ and $Y^{B}$ both being above
their respective $(1-p)$th marginal quantiles, compares to $p$ (the value of
this probability under perfect dependence of $Y^{A}$ and $Y^{B}$) and relative
to $p^{2}$ (the value under independence of $Y^{A}$ and $Y^{B}$). Under
positive dependence we would expect that

$p^{2}<\mathbb{P}\\{Y^{A}>F_{A}^{-1}(1-p),Y^{B}>F_{B}^{-1}(1-p)\\}<p.$ (2.21)

Coles et al.,  (1999) formalise this intuition to define the measure of
extremal dependence as $p\rightarrow 0$, i.e., as we look above increasing
quantiles. Specifically they take

$\chi=\lim_{p\rightarrow
0}\mathbb{P}\\{Y^{B}>F_{B}^{-1}(1-p)~{}\mid~{}Y^{A}>F_{A}^{-1}(1-p)\\}$

where $\chi\in[0,1]$. Increasing values of $\chi$ corresponding to stronger
extremal dependence, and $\chi=1$ corresponds to perfect dependence between
$Y^{A}$ and $Y^{B}$. Thus $\chi$ is the limiting probability of one variable
being extreme given that the other is equally extreme. If $\chi\in(0,1]$, we
say $Y^{A}$ and $Y^{B}$ are asymptotically dependent, with there being a non-
zero probability of $Y^{B}$ being large when $Y^{A}$ is large at extreme
levels. Though the class of extremal dependence where $\chi>0$ is widely
studied, this only corresponds to cases where the joint probability in (2.21)
is of $O(p)$, i.e., decays as a multiple of $p$ as $p\rightarrow 0$. We find
that $\chi=0$ in all other dependence cases as well as when $Y^{A}$ and
$Y^{B}$ are actually independent, this class of variables is known as being
asymptotically independent, and $\chi$ doesn’t give us any information on the
level of asymptotic independence. We need a more refined measure of extremal
dependence than $\chi$ to enable us to separate between when there is some
dependence of large values and independence of $Y^{A}$ and $Y^{B}$. Coles et
al.,  (1999) also define

$\bar{\chi}=\lim_{p\rightarrow
0}\frac{2\log\mathbb{P}\\{Y^{A}>F_{A}^{-1}(1-p)\\}}{\log\mathbb{P}\\{Y^{B}>F_{B}^{-1}(1-p),Y^{A}>F_{A}^{-1}(1-p)\\}}-1.$
(2.22)

where $\bar{\chi}\in(-1,1]$. Here $\bar{\chi}=1$ and $-1<\bar{\chi}<1$
correspond to asymptotic dependence and asymptotic independence, respectively.
When $\bar{\chi}=0$ this shows there is no dependence in the tails of
$(Y^{A},Y^{B})$ as it arises when $Y^{A}$ and $Y^{B}$ are independent, with
$0<\bar{\chi}<1$ and $-1<\bar{\chi}<0$ indicating positive and negative
dependence in the joint tails of $Y^{A}$ and $Y^{B}$ respectively.

To assess inter-site dependence in extreme skew surges we evaluate these
dependence measures using empirical estimates of the associated probabilities
using texmex R package (Southworth et al., , 2020). Specifically, we use skew
surge daily maxima for each pairwise combination of the four study sites,
using data on the same day and with lags of $\pm 1$ day to account for time
lags between the peak of surge reaching each site, when events last multiple
days. Here we have lags $t=1$ and $t=-1$ so that site $A$ is one day ahead or
behind site $B$, respectively. Since the variables are not identically
distributed, due to seasonality for example, this can effect the evaluation of
$\chi$ and $\bar{\chi}$. We address this potential concern by also using the
marginal distributional model of D’Arcy et al.,  (2022) $F_{Y}^{(d,j,k)}$,
given by expression (2.7), to account for this through a transform the
variables to identical uniform margins and then re-evaluate these measures.
These results are discussed in Section 3.3.

##  3 Results

###  3.1 Introduction

We now present the results from applying the extreme skew surge models
discussed in Sections 2.4 and 2.5, in Sections 3.2 and 3.3, respectively to
data from our four study sites. Here we define extreme skew surges as being
exceedances of the monthly 0.95 quantile, as in D’Arcy et al.,  (2022). All
models are fit in a likelihood framework, with 95% confidence intervals
provided for parameter estimates based on the hessian, i.e., using asymptotic
normality of maximum likelihood estimators. The likelihood is constructed
under the assumption that extreme skew surges are temporally independent for
single site inference, but also that observations at different sites are
independent for spatial pooling. These are not unreasonable assumptions for
model selection, the former being found as a reasonable approximation in
D’Arcy et al.,  (2022) as the extremal index is near one for large levels and
the validity of the latter being assessed before we apply any spatial pooling.
We compare models using AIC and BIC scores; these are commonly used measures
of the quality of a statistical model for a particular data set relative to
the parsimony of the model. The chosen best fitting model should minimise
these scores. Recall that all estimates presented here are after the mean sea
level trends have been removed. An estimated change here means that the change
is the same as in the mean sea level, so negative trend estimates correspond
to the extreme sea levels not rising as fast as the mean the level at the
site.

###  3.2 Single-site Analysis

We fit the models of Section 2.4 to the GPD rate and scale parameters for
extreme skew surges individually at each site. We start with the threshold
exceedance probability parameter, $\lambda$, fitting Models $R0-4$; AIC/BIC
scores and estimates of $\delta_{\lambda}^{(\tilde{k})}$,
$\delta_{\lambda,s}^{(\tilde{k})}$, $\delta_{\lambda,s}^{(m)}$ and
$\delta_{\lambda}^{(m)}$ are given in Table 1. Since the parameter estimates
are not intuitive, we consider the change in exceedance probability with the
particular covariate of interest for the annual trends in Models $R1$ and
$R3$.

| Heysham | Lowestoft | Newlyn | Sheerness
---|---|---|---|---
Model $R0$
AIC | 12234.21 | 15312.08 | 24498.77 | 9286.58
BIC | 12275.89 | 15354.88 | 24543.93 | 9326.94
Model $R1$
$\delta_{\lambda}^{(\tilde{k})}$ | 0.091 (-0.008, 0.191) | -0.061 (-0.150, 0.028) | 0.215 (0.154,0.276) | -0.114 (-0.013, -0.010)
AIC | 12232.89 | 15312.26 | 24453.12 | 9283.96
BIC | 12282.91 | 15363.61 | 24507.31 | 9332.40
Model $R2$
$\delta_{\lambda,1}^{(\tilde{k})}$ | 0.161 (-0.033, 0.335) | 0.063 (-0.106, 0.232) | 0.114 (-0.011, 0.238) | -0.032 (-0.228, 0.164)
$\delta_{\lambda,2}^{(\tilde{k})}$ | 0.034 (-0.161, 0.230) | -0.141 (-0.322, 0.040) | 0.197 (0.077, 0.316) | -0.250 (-0.468, -0.032)
$\delta_{\lambda,3}^{(\tilde{k})}$ | 0.207 (0.013, 0.400) | -0.094 (-0.266, 0.078) | 0.209 (0.089, 0.328) | -0.189 (-0.405, 0.026)
$\delta_{\lambda,4}^{(\tilde{k})}$ | -0.047 (-0.261, 0.167) | -0.081 (-0.264, 0.102) | 0.338 (0.217, 0.460) | -0.021 (-0.221, 0.178)
AIC | 12235.30 | 15315.29 | 24452.52 | 9286.48
BIC | 12310.32 | 15392.33 | 24533.80 | 9359.14
Model $R3$
$\delta_{\lambda}^{(m)}$ | 0.204 (0.074, 0.334) | -0.012 (-0.12, 0.427) | 0.336 (0.245, 0.427) | -0.164 (-0.304, -0.024)
AIC | 12227.14 | 15314.04 | 24451.34 | 9283.20
BIC | 12277.16 | 15365.39 | 24505.53 | 9331.64
Model $R4$
$\delta_{\lambda,1}^{(m)}$ | 0.256 (-0.002, 0.514) | 0.103 (-0.107, 0.312) | 0.135 (-0.058, 0.329) | -0.079 (-0.340, 0.181)
$\delta_{\lambda,2}^{(m)}$ | 0.111 (-0.143, 0.365) | -0.067 (-0.282, 0.149) | 0.322 (0.144, 0.501) | -0.374 (-0.672, -0.076)
$\delta_{\lambda,3}^{(m)}$ | 0.416 (0.167, 0.665) | -0.048 (-0.259, 0.163) | 0.393 (0.212, 0.574) | -0.273 (-0.568, 0.022)
$\delta_{\lambda,4}^{(m)}$ | 0.010 (-0.274, 0.293) | -0.040 (-0.262, 0.182) | 0.478 (0.300, 0.655) | -0.008 (-0.253, 0.269)
AIC | 12227.92 | 15318.49 | 24450.27 | 9284.69
BIC | 12302.94 | 15395.53 | 24531.56 | 9357.34

Table 1: Parameter estimates for the Models $R1-R4$ with AIC and BIC scores
for each model fit at each site (including Model $R0$). The minimum AIC/BIC
scores are highlighted in bold for each site.

We find that Model $R3$ minimises AIC at Heysham, Lowestoft and Sheerness,
whilst at Newlyn Model $R4$ is preferable. The BIC is minimised by Model $R3$
at Newlyn, but elsewhere Model $R0$ is favourable. This suggests that if any
longer-term trends are included in the model to capture changes in the rate of
extreme events (relative to mean sea level), GMT should be used as a covariate
as opposed to the year.

We look at Models $R1$ and $R3$ in more detail, these have a fixed trend
parameter within the year with respect to year and GMT, respectively. At
Newlyn we find a significant increasing trend for both models, since the
confidence intervals do not contain zero. We also find positive trends at
Heysham, but only the GMT trend in Model $R3$ is significant. Neither trends
are significant at Lowestoft, but we find a significant decreasing trend for
both models at Sheerness. Figure 2 shows histograms of the estimates of
$\Delta_{\lambda}^{(\tilde{k})}$ and $\Delta_{\lambda}^{(m)}$ (defined in
Section 2.4), based on all combinations of day $d$ and peak tide $x$, so these
do not account for uncertainty in $\delta^{(\tilde{k})}_{\lambda}$ or
$\delta_{\lambda}^{(m)}$ estimates but are simply a reflection that the rate
of threshold exceedance varies over the short term. For Model $R1$, we find an
increase in $\lambda_{d,x,\tilde{k}}$ over 100 years at Newlyn, with
$\max\Delta^{(\tilde{k})}_{\lambda}=3\%$, so that the exceedance probability
almost doubles from 3.5% to 6.5% in 1920-2020. However, we observe decreases
in exceedance probability at Sheerness. For Model $R3$, we also find an
increase in exceedance probability with a $1^{\circ}$C increase in GMT at
Newlyn, where $\max\Delta_{\lambda}^{(m)}=3\%$, but a negative trend at
Sheerness. If the trends were statistically significant at Heysham and
Lowestoft, the exceedance probability would be increasing and decreasing with
both trend parameters, respectively.

Figure 2: Histograms of $\Delta_{\lambda}^{(\tilde{k})}$ over 100 years (left)
and $\Delta_{\lambda}^{(m)}$ with a $1^{\circ}$C increase in GMT, as
percentages, for all day $d$ and peak tide $x$ combinations at Heysham (cyan),
Lowestoft (grey), Newlyn (blue) and Sheerness (green).

Now we look at Models $R2$ and $R4$ with season-specific trend parameters for
year and GMT, respectively. The trends at Newlyn are significant in both
models, except for winter, whilst at Heysham only the increasing trends in
summer are significant. None of the seasonal trends are significant at
Lowestoft but we find a significant decreasing trend for spring in both models
at Sheerness. As for Models $R1$ and $R3$, we obtain a variety of results
across sites; Newlyn has an increasing exceedance probability with year and
GMT in all seasons, whilst Sheerness has a decreasing trend parameter for all
seasons, with the greatest decline in spring. However, for Heysham and
Lowestoft, we have a mixture of positive and negative parameters throughout
the year. The confidence intervals for the four parameter estimates in Models
$R2$ and $R4$ at Heysham, Lowestoft and Sheerness overlap, suggesting that
there isn’t significant within-year variation of the longer-term trend
parameters so that the simpler Models $R1$ and $R3$ are sufficient. Whilst at
Newlyn, this overlap is small (see Figure 3). Here, we find the greatest trend
in autumn, which is not concerning for extreme sea level estimation since the
most extreme sea levels occur in winter (D’Arcy et al., , 2022), but using
Models $R1$ and $R3$ with common trend parameters across the year could result
in overestimation of the trends in winter, hence influencing sea level return
level estimation.

Figure 3: Confidence intervals for parameter estimates
$\hat{\delta}_{\lambda,s}^{(\tilde{k})}$ (left) and
$\hat{\delta}_{\lambda,s}^{(m)}$ (right) at Newlyn for $s=1,2,3,4$ denoting
winter, spring, summer and autumn, respectively.

Next, we consider models for the scale parameter at each site individually
(Models $S0-4$, introduced in Section 2.4). Table 2 shows the parameter
estimates for each model, along with AIC and BIC scores. Models $S1$ and $S3$
have a single parameter denoting a common long-term trend across the year, but
neither of these are an improvement on Model $S0$ (without a long-term trend)
at any site. All of the 95% confidence intervals for
$\delta^{(\tilde{k})}_{\lambda}$ or $\delta_{\lambda}^{(m)}$ estimates contain
zero, suggesting these trends are not significant. If we ignore this
uncertainty, the point estimates suggest small changes in the scale parameter;
at Heysham and Lowestoft our results show a decrease with both year and GMT,
suggesting that the magnitude of extreme skew surge events are getting smaller
with anthropogenic climate change effects. In the 100 year period 1920-2020 at
Newlyn, the point estimate $\delta_{\sigma}^{(\tilde{k})}$ corresponds to an
increase in mean excesses (see expression (2.5)) of 2mm (relative to a mean of
94mm in 1920), whilst at Sheerness in the years of observation 1980-2016 this
corresponds to an increase of 10mm relative to a mean of 125mm in 1980. Notice
there is overlap in the parameter estimates for
$\delta_{\sigma}^{(\tilde{k})}$ and $\delta_{\sigma}^{(m)}$ across sites; in
Section 3.3 we fit similar model with these trend parameters common across
sites (see Figure 4).

| Heysham | Lowestoft | Newlyn | Sheerness
---|---|---|---|---
Model $S0$
AIC | -3091.53 | -3672.07 | -10152.63 | -2974.317
BIC | -3064.77 | -3644.20 | -10122.42 | -2948.854
Model $S1$
$\delta_{\sigma}^{(\tilde{k})}$ | -0.009 (-0.032, 0.013) | -0.006 (-0.024, 0.011) | 0.001 (-0.003, 0.005) | 0.016 (-0.013, 0.044)
AIC | -3088.558 | -3670.55 | -10150.80 | -2973.05
BIC | -3056.448 | -3637.11 | -10114.54 | -2942.50
Model $S2$
$\delta_{1}^{(\tilde{k})}$ | 0.022 (-0.034, 0.078) | -0.041 (-0.090, 0.007) | 0.004 (-0.009, 0.016) | 0.023 (-0.032, 0.078)
$\delta_{2}^{(\tilde{k})}$ | 0.022 (-0.014, 0.059) | -0.030 (-0.055, -0.006) | 0.006 (-0.002, 0.014) | -0.001 (-0.036, 0.035)
$\delta_{3}^{(\tilde{k})}$ | -0.025 (-0.051, 0.001) | 0.012 (-0.010, 0.034) | -0.003 (-0.009, 0.003) | 0.023 (-0.010, 0.055)
$\delta_{4}^{(\tilde{k})}$ | -0.035 (-0.079, 0.008) | -0.015 (-0.053, 0.023) | 0.001 (-0.008, 0.011) | 0.008 (-0.039, 0.054)
AIC | -3095.28 | -3674.12 | -10146.27 | -2971.19
BIC | -3047.12 | -3623.96 | -10091.89 | -2925.36
Model $S3$
$\delta_{\sigma}^{(m)}$ | -0.011 (-0.033, 0.011) | -0.008 (-0.026, 0.009) | -0.001 (-0.008, 0.006) | 0.006 (-0.020, 0.032)
AIC | -3088.42 | -3670.90 | -10149.07 | -2972.43
BIC | -3056.32 | -3637.46 | -10112.82 | -2941.87
Model $S4$
$\delta_{1}^{(m)}$ | 0.036 (-0.027, 0.099) | -0.050 (-0.105, 0.004) | -0.0003 (-0.021, 0.020) | 0.025 (-0.042, 0.091)
$\delta_{2}^{(m)}$ | 0.029 (-0.012, 0.070) | -0.037 (-0.061, -0.012) | 0.005 (-0.008, 0.018) | -0.015 (-0.054, 0.023)
$\delta_{3}^{(m)}$ | -0.027 (-0.054, -0.00009) | 0.017 (-0.006, 0.039) | -0.003 (-0.013, 0.006) | 0.013 (-0.018, 0.045)
$\delta_{4}^{(m)}$ | -0.030 (-0.081, 0.021) | -0.024 (-0.066, 0.017) | -0.005 (-0.021, 0.010) | -0.006 (-0.053, 0.040)
AIC | -3093.73 | -3677.95 | -10145.39 | -2970.79
BIC | -3045.56 | -3627.79 | -10091.01 | -2924.96

Table 2: Parameter estimates for the Models $S1-4$ with AIC and BIC scores for
each model fit at each site (including Model $S0$). The minimum AIC/BIC scores
are highlighted in bold for each site.

Models $S2$ and $S4$ have four additional parameters relative to Model $S0$,
these denote a separate trend for each season with respect to year and GMT.
AIC and BIC are still minimised by Model $S0$, except AIC scores for Heysham
and Lowestoft, which favour Models $S2$ and $S4$, respectively. However, the
four confidence intervals overlap at each site, suggesting a fixed trend
within a year is sufficient. At Heysham, the overlap across all seasons is
small but there is considerable overlap between winter and spring, with
positive trend parameters, and likewise for summer and autumn with negative
trend parameters for both models. If these trends were statistically
significant it would suggest that the magnitude of extreme skew surges is
increasing with increases in GMT in December-April, but decreasing for the
rest of the year. Given the timing of the most extreme events, this could be
important if statistically significant.

Figure 4: Confidence intervals for parameter estimates
$\hat{\delta}_{\lambda}^{(\tilde{k})}$ (left) and
$\hat{\delta}_{\sigma}^{(\tilde{k})}$ (right) for all sites.

###  3.3 Spatial Pooling

We present the results from pooling information across sites, for the long-
term trend parameters, when refitting the models of Section 2.4. Before
pooling information, we use the dependence measures discussed in Section 2.5
to check if it is reasonable to assume each pair of sites are independent in
their extreme skew surge values. We estimate the dependence measures for the
observed skew surges and a standardised transformation of them to remove
sources of within-year non-stationarity via mapping to uniform margins through
the distribution function (2.7). The results are shown in Table 3. For most
combinations of sites at lags $t=-1,0,1$ the dependence is weak, with
Kendall’s $\tau$ of approximately $0.3$ for most pairs, except for Newlyn with
the east coast site where Kendall’s $\tau$ is near 0, whilst for Lowestoft and
Sheerness this value is approximately 0.5. The effect of de-seasonalising the
data (by transforming to uniform margins) has typically decreased dependence
by at most $5\%$. Examining the extremal dependence measure $\chi$, we find
that the de-seasonalised estimate is down from 0.3 to approximately 0.15 for
most pairs, though unchanged for Lowestoft and Sheerness. With the exception
of Lowestoft and Sheerness, it is not unreasonable to make an independence in
extremes approximation for the data. Finally, looking at $\bar{\chi}$ we find
$\bar{\chi}<1$ for all pairs, giving evidence of asymptotic independence with
weak dependence in the observed tails of the variables. The strongest
dependence is found between Lowestoft and Sheerness at lag $t=0$. This is not
surprising since these sites are close in proximity, with extreme skew surges
progressing south down the east coast through Lowestoft onto Sheerness.
Therefore they are highly likely to be affected by the same storms. Despite
this pair of sites giving clear evidence of dependence, we continue under the
belief that it is reasonable to assume skew surge daily maxima at all pairs of
sites are sufficiently close to being independent for the purposes of spatial
pooling.

| HEY-LOW | HEY-NEW | HEY-SHE | LOW-NEW | LOW-SHE | NEW-SHE
---|---|---|---|---|---|---
lag | -1 | 0 | 1 | -1 | 0 | 1 | -1 | 0 | 1 | -1 | 0 | 1 | -1 | 0 | 1 | -1 | 0 | 1
Observations
$\tau$ | 0.133 | 0.160 | 0.309 | 0.287 | 0.322 | 0.259 | 0.153 | 0.149 | 0.298 | 0.089 | 0.040 | 0.034 | 0.155 | 0.510 | 0.238 | 0.137 | 0.168 | 0.196
$\chi$ | 0.095 | 0.129 | 0.270 | 0.127 | 0.145 | 0.076 | 0.092 | 0.111 | 0.259 | 0.017 | 0 | 0 | 0.145 | 0.509 | 0.200 | 0.054 | 0.077 | 0.121
$\bar{\chi}$ | 0.200 | 0.251 | 0.424 | 0.249 | 0.276 | 0.160 | 0.195 | 0.224 | 0.412 | 0.040 | -0.018 | -0.055 | 0.274 | 0.645 | 0.344 | 0.120 | 0.158 | 0.237
Transform to Uniform(0,1)
$\tau$ | 0.103 | 0.130 | 0.289 | 0.285 | 0.318 | 0.244 | 0.108 | 0.102 | 0.262 | 0.086 | 0.036 | 0.028 | 0.143 | 0.523 | 0.228 | 0.139 | 0.173 | 0.200
$\chi$ | 0.026 | 0.040 | 0.180 | 0.103 | 0.122 | 0.053 | 0.056 | 0.036 | 0.173 | 0 | 0 | 0 | 0.095 | 0.494 | 0.174 | 0.003 | 0.016 | 0.050
$\bar{\chi}$ | 0.069 | 0.100 | 0.321 | 0.215 | 0.236 | 0.114 | 0.123 | 0.084 | 0.313 | -0.012 | -0.107 | -0.134 | 0.198 | 0.634 | 0.316 | 0.003 | 0.035 | 0.114

Table 3: Kendall’s $\tau$, $\chi$ and $\bar{\chi}$ measures of dependence for
daily maximum skew surge observations at pairs of sites. We show the
dependence over lags -1 (LHS site is 1 day behind RHS), 0 and 1 (LHS site is 1
day ahead of RHS); in bold we show the largest dependence over these lags.
$\chi$ and $\bar{\chi}$ are measures of extremal dependence for exceedances of
the 0.95 quantile.

Firstly, we focus on pooling information across sites regarding the long-term
trend parameters with respect to year $k$ and GMT $m_{k}$. Figure 4 shows
there is considerable overlap in the confidence intervals for
$\hat{\delta}_{\lambda}^{(\tilde{k})}$ at Lowestoft and Sheerness; similarly
there is some overlap for Heysham and Newlyn. Although pooling information
across randomly selected subsets of sites should be discouraged, here we note
that the pairs of sites with similarities are on different coastlines, so we
explore pooling over sites on the east coast (Lowestoft and Sheerness) and
separately on the west coast (Heysham and Newlyn). Here, we consider refitting
Models $R1$ and $R3$ (i.e., a fixed trend parameter within a year) with common
trend parameters $\delta_{\lambda}^{(\tilde{k})}$ and $\delta_{\lambda}^{(m)}$
between the pairs of sites. We obtain negative trends parameters
$\hat{\delta}_{\lambda}^{(\tilde{k})}=-0.084\;(-0.151,-0.016)$ and
$\hat{\delta}_{\lambda}^{(m)}=-0.070\;(-0.156,0.015)$ for Sheerness and
Lowestoft, whilst at Newlyn and Heysham we obtain statistically significant
positive trend parameters
$\hat{\delta}_{\lambda}^{(\tilde{k})}=0.180\;(0.128,0.231)$ and
$\hat{\delta}_{\lambda}^{(m)}=0.285\;(0.208,0.359)$. Both of these models are
an improvement on the previous results, where a separate long-term trend
parameter is used for each site; the model with a yearly trend parameter
reduces AIC by 48 and the BIC by 0.5, whilst the model with a GMT parameter
reduces AIC and BIC by 54 and 6, respectively. This highlights the importance
of sharing information spatially.

There is also information to be gained from sharing spatial information about
long-term trends in the scale parameter since there is considerable overlap in
the confidence intervals for $\hat{\delta}_{\sigma}^{(\tilde{k})}$ (see Figure
4) and $\hat{\delta}_{\sigma}^{(m)}$ (see Table 2). We refit the models of
Section 2.4 for the scale parameter with common longer-term trend parameters
across sites, but neither parameter estimates are significant. We find that
$\hat{\delta}_{\sigma}^{(\tilde{k})}=4.8\times 10^{-4}$, corresponding to an
increase in scale parameter of 0.48mm over 1915-2020. For GMT
$\hat{\delta}_{\sigma}^{(\tilde{k})}=-0.0024$, a 24mm decrease in scale
parameter with a $1^{\circ}$C increase in temperature. Neither of these models
improve the fit relative to having no long-term trends (in addition to those
in mean sea level), although the AIC scores are close. We also fit a model
similar to that of Models $S2$ and $S4$ so there is a common seasonal trend
across sites, with respect to year and GMT but find that neither of these
improve model fit. This agrees with our single-site results of Section 3.2
where we found no evidence of changes in the magnitude of extreme skew surge
events with respect to year or temperature.

##  4 Discussion

We have presented a framework to investigate the effects of anthropogenic
climate change on extreme skew surges as any increases in the magnitude or
frequency in these events can have catastrophic consequences if not included
in extreme sea level estimation for coastal flood defence design. These trends
can be different to those observed in the main body of the data, such as mean
sea level rise. We use year and GMT as covariates in our statistical model for
extreme event occurrence, building on a model developed by D’Arcy et al.,
(2022) that accounts for seasonality and skew surge-peak tide dependence.
Recall that our results are relative to the mean sea level trend in 2017 so
this would need to be added onto any sea level return level estimates when
used in practice. We show that there is evidence of an increase in the
probability of an extreme skew surge event with GMT increases at Heysham and
Newlyn, but evidence of a decrease in the likelihood of these events at
Lowestoft and Sheerness. We do not find any significant changes in the
magnitude of extreme skews surges, i.e., in the scale parameter, and hence in
the mean of the skew surge excesses of the threshold. We also looked at any
seasonal variations in these longer-term trends and found evidence of a
greater change in exceedance probability with year and GMT in autumn compared
to the other seasons, at Newlyn. The ideas presented in this paper could be
applied to more locations, but also to other environmental variables to
investigate trends in extreme values.

We demonstrate the advantages of pooling information across sites, although
this is only primarily illustrative since we consider just four sites here.
There are 44 sites on the UK National Tide Gauge Network where this
methodology could be extended. It would be interesting to apply our
methodology within a spatial framework, for example in regional frequency
analysis where sites in a homogeneous region not only have a common shape
parameter, but also common longer-term trends due to anthropogenic climate
change.

Skew surges are also believed to change over decadal time scales with climate
indices. The North Atlantic Oscillation index (NAO) describes such time scale
changes in regional weather systems, so is believed to impact storm surges,
and thus skew surge. Araújo and Pugh,  (2008) find a negative correlation
between storm surge and air pressure patterns, using NAO. It would be
interesting to explore how adding an NAO covariate into the GPD for extreme
skew surges would change model fit.

##  5 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.

## References

  * Araújo and Pugh, (2008) Araújo, I. B. and Pugh, D. T. (2008).  Sea levels at Newlyn 1915–2005: analysis of trends for future flooding risks.  Journal of Coastal Research, (24):203–212.
  * Baranes et al., (2020) Baranes, H., Woodruff, J., Talke, S., Kopp, R., Ray, R., and DeConto, R. (2020).  Tidally driven interannual variation in extreme sea level frequencies in the Gulf of Maine.  Journal of Geophysical Research: Oceans, 125(10):e2020JC016291.
  * Batstone et al., (2013) Batstone, C., Lawless, M., Tawn, J. A., Horsburgh, K., Blackman, D., McMillan, A., Worth, D., Laeger, S., and Hunt, T. (2013).  A UK best-practice approach for extreme sea-level analysis along complex topographic coastlines.  Ocean Engineering, 71:28–39.
  * Coles, (2001) Coles, S. G. (2001).  An Introduction to Statistical Modeling of Extreme Values.  Springer, London.
  * Coles et al., (1999) Coles, S. G., Heffernan, J., and Tawn, J. A. (1999).  Dependence measures for extreme value analyses.  Extremes, 2(4):339–365.
  * D’Arcy et al., (2022) D’Arcy, E., Tawn, J. A., Joly, A., and Sifnioti, D. E. (2022).  Accounting for seasonality in extreme sea level estimation.  Under review.  arXiv:2207.09870.
  * Davison et al., (2012) Davison, A. C., Padoan, S. A., and Ribatet, M. (2012).  Statistical modeling of spatial extremes.  Statistical science, 27(2):161–186.
  * Dixon and Tawn, (1999) Dixon, M. J. and Tawn, J. A. (1999).  The effect of non-stationarity on extreme sea-level estimation.  Journal of the Royal Statistical Society: Series C, 48(2):135–151.
  * Dixon et al., (1998) Dixon, M. J., Tawn, J. A., and Vassie, J. M. (1998).  Spatial modelling of extreme sea-levels.  Environmetrics, 9(3):283–301.
  * Eastoe and Tawn, (2009) Eastoe, E. F. and Tawn, J. A. (2009).  Modelling non-stationary extremes with application to surface level ozone.  Journal of the Royal Statistical Society: Series C, 58(1):25–45.
  * Egbert and Ray, (2017) Egbert, G. D. and Ray, R. D. (2017).  Tidal prediction.  Journal of Marine Research, 75(3):189–237.
  * Environment Agency, (2018) Environment Agency (2018).  Coastal Flood Boundary Conditions for the UK: update 2018. Technical summary report.  https://environment.data.gov.uk/dataset/6e856bda-0ca9-404f-93d7-566a2378a7a8.  Accessed 01/10/21.
  * Fawcett and Walshaw, (2007) Fawcett, L. and Walshaw, D. (2007).  Improved estimation for temporally clustered extremes.  Environmetrics, 18(2):173–188.
  * Ferro and Segers, (2003) Ferro, C. A. T. and Segers, J. (2003).  Inference for clusters of extreme values.  Journal of the Royal Statistical Society: Series B, 65(2):545–556.
  * Graff, (1978) Graff, J. (1978).  Concerning the recurrence of abnormal sea levels.  Coastal Engineering, 2:177–187.
  * Howard and Williams, (2021) Howard, T. and Williams, S. D. P. (2021).  Towards using state-of-the-art climate models to help constrain estimates of unprecedented UK storm surges.  Natural Hazards and Earth System Sciences, 21(12):3693–3712.
  * Huser and Genton, (2016) Huser, R. and Genton, M. G. (2016).  Non-stationary dependence structures for spatial extremes.  Journal of agricultural, biological, and environmental statistics, 21(3):470–491.
  * Leadbetter et al., (1983) Leadbetter, M., Lindgren, G., and Rootzén, H. (1983).  Extremes and Related Properties of Random Sequences and Processes.  Springer-Verlag, New York.
  * Ledford and Tawn, (2003) Ledford, A. W. and Tawn, J. A. (2003).  Diagnostics for dependence within time series extremes.  Journal of the Royal Statistical Society: Series B, 65(2):521–543.
  * Northrop et al., (2017) Northrop, P. J., Attalides, N., and Jonathan, P. (2017).  Cross-validatory extreme value threshold selection and uncertainty with application to ocean storm severity.  Journal of the Royal Statistical Society: Series C (Applied Statistics), 66(1):93–120.
  * Pugh and Vassie, (1978) Pugh, D. and Vassie, J. (1978).  Extreme sea levels from tide and surge probability.  Coastal Engineering, 16:911–930.
  * Pugh and Woodworth, (2014) Pugh, D. and Woodworth, P. (2014).  Sea-Level Science: Understanding Tides, Surges, Tsunamis and Mean Sea-Level Changes.  Cambridge University Press.
  * Smith et al., (1997) Smith, R. L., Tawn, J. A., and Coles, S. G. (1997).  Markov chain models for threshold exceedances.  Biometrika, 84(2):249–268.
  * Smith and Weissman, (1994) Smith, R. L. and Weissman, I. (1994).  Estimating the extremal index.  Journal of the Royal Statistical Society: Series B, 56(3):515–528.
  * Southworth et al., (2020) Southworth, H., Heffernan, J. E., and Metcalfe, P. D. (2020).  texmex: Statistical modelling of extreme values.  R package version 2.4.8.
  * Tawn, (1988) Tawn, J. A. (1988).  An extreme-value theory model for dependent observations.  Journal of Hydrology, 101(1-4):227–250.
  * Tawn, (1992) Tawn, J. A. (1992).  Estimating probabilities of extreme sea-levels.  Journal of the Royal Statistical Society: Series C, 41(1):77–93.
  * Wadsworth and Tawn, (2012) Wadsworth, J. and Tawn, J. (2012).  Likelihood-based procedures for threshold diagnostics and uncertainty in extreme value modelling.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 74(3):543–567.
  * Wahl et al., (2013) Wahl, T., Haigh, I. D., Woodworth, P. L., Albrecht, F., Dillingh, D., Jensen, J., Nicholls, R. J., Weisse, R., and Wöppelmann, G. (2013).  Observed mean sea level changes around the North Sea coastline from 1800 to present.  Earth-Science Reviews, 124:51–67.
  * Weisse et al., (2014) Weisse, R., Bellafiore, D., Menéndez, M., Méndez, F., Nicholls, R. J., Umgiesser, G., and Willems, P. (2014).  Changing extreme sea levels along European coasts.  Coastal engineering, 87:4–14.
  * Williams et al., (2016) Williams, J., Horsburgh, K. J., Williams, J. A., and Proctor, R. N. (2016).  Tide and skew surge independence: new insights for flood risk.  Geophysical Research Letters, 43(12):6410–6417.
  * Wong et al., (2022) Wong, T. E., Sheets, H., Torline, T., and Zhang, M. (2022).  Evidence for increasing frequency of extreme coastal sea levels.  Frontiers in Climate, 4.
  * Woodworth et al., (2011) Woodworth, P. L., Menéndez, M., and Roland Gehrels, W. (2011).  Evidence for century-timescale acceleration in mean sea levels and for recent changes in extreme sea levels.  Surveys in geophysics, 32(4):603–618.
  * Woodworth and Player, (2003) Woodworth, P. L. and Player, R. (2003).  The permanent service for mean sea level: An update to the 21st century.  Journal of Coastal Research, 19(2):287–295.
  * Zsamboky et al., (2011) Zsamboky, M., Fernández-Bilbao, A., Smith, D., Knight, J., and Allan, J. (2011).  Impacts of climate change on disadvantaged UK coastal communities.  Joseph Rowntree Foundation, pages 1–63.


# Crossover of spin Hall to quantum anomalous Hall effect in PbC/MnSe
heterostructures††thanks: Project supported by the National Natural Science
Foundation of China (Grant No. 61604106), and Shandong Provincial Natural
Science Foundation (Grant No. ZR2014FL025).

L. Chen

(Received July 24, 2019, in final form September 26, 2019)

###### Abstract

In this paper, we study topological properties and Hall conductivities in
PbC/MnSe heterostructure under the illumination of a circularly polarized
light. At high frequency regime, energy gap, Chern numbers, and Hall
conductivities are studied based on the Floquet theory and Green’s function
formalism, respectively. The interplay between spin orbit coupling and light
leads to topological phase transition between anomalous Hall states and spin
Hall states, which is related to the emission and absorption of two virtual
photons. The anomalous Hall conductivities are dependent on polarization of
light, while the spin Hall conductivity is independent.

Key words: Hall, light, Flouqet theory, topology

###### Abstract

Ó ñòàòò äîñëäæóþòüñÿ òîïîëîãí âëàñòèâîñò òà õîëëâñüê ïðîâäíîñò ó
ãåòåðîñòðóêòóð PbC/MnSe ïðè îñâòëåíí öèðêóëÿðíî ïîëÿðèçîâàíèì ñâòëîì. Ó
âèñîêîàñòîòíîìó ðåæèì äîñëäæóþòüñÿ çàáîðîíåíà çîíà, èñëà åðíà òà õîëëâñüê
ïðîâäíîñò íà îñíîâ òåîð¿ Ôëîêå ôîðìàëçìó ôóíêöé Ãðíà, âäïîâäíî. Âçàìîçâ’ÿçîê
ìæ ñïí-îðáòàëüíîþ âçàìîäþ ñâòëîì ïðèâîäèòü äî òîïîëîãíîãî ôàçîâîãî ïåðåõîäó ìæ
àíîìàëüíèìè ñòàíàìè Õîëëà òà ñïí-õîëëâñüêèìè ñòàíàìè, ùî, â ñâîþ åðãó,
ïîâ’ÿçàíî ç âèïðîìíþâàííÿì ïîãëèíàííÿì äâîõ âðòóàëüíèõ ôîòîíâ. Àíîìàëüí
ïðîâäíîñò Õîëëà çàëåæàòü âä ïîëÿðèçàö¿ ñâòëà, â òîé àñ ÿê ñïí-õîëëâñüêà
ïðîâäíñòü íåçàëåæíîþ.

Ключов слова: Õîëë, ñâòëî, òåîðÿ Ôëîêå, òîïîëîãÿ

##  1 Introduction

In the recent years, topological phases and topological phase transitions have
been one of the most important topics in modern condensed matter physics [1,
2]. The intrinsic anomalous Hall and spin Hall effect have attracted enormous
attention in the past decade [3, 4]. The quantum anomalous Hall effect has
been theoretically proposed [5, 6] and experimentally realized [7, 8] in
magnetically doped topological insulator films. The quantum spin Hall effect
was theoretically predicted [9] and experimentally observed [10] in HgTe
quantum wells. Li et al. [11] investigated the role of disorder in the
topological insulator, which possesses a pair of helical edge states with
opposing spins moving in opposite directions and exhibits the phenomenon of
quantum spin Hall effect. The Kane and Mele model defined on a hexagonal
lattice can exhibit transition between quantum spin Hall effect phase and a
simple insulator due to Rashba spin orbit coupling [12, 13]. The quantum spin
Hall insulator phase has been investigated in the graphene system with a
coexistence of Coulomb interaction, staggered potential, and intrinsic spin-
orbit coupling [14]. In two-dimensional photonic crystal, a conventional
insulator phase, a quantum spin Hall phase, or a quantum anomalous Hall phase
can be realized by simply adjusting the geometric parameters and magnetic
field [15, 16]. Kort-Kamp [17] unveiled topological phase transitions in the
photonic spin Hall effect in the graphene family materials. The spin Hall
effect can be realized in a Bose-Einstein condensate of neutral atoms
interacting via the magnetic dipole-dipole interactions [18]. A $1$D dynamical
version of the quantum spin Hall effect was implemented in an optical
superlattice with ultracold bosonic atoms [19]. Omnidirectional spin Hall
effect was studied in a Weyl spin-orbit-coupled atomic gas [20]. The edge
state of quantum spin Hall insulators in topological Dirac and Weyl semimetals
has been studied [21]. Recently, interaction-driven anomalous quantum Hall
state was theoretical predicted [22, 23].

Another interesting direction of topological states and and topological phase
transitions studied in the recent years arises from the nonequilibrium
engineering Hall effect under the influence of a periodic drive. An off
resonant circularly polarized light has been found to induce a topological
term in periodically driven quantum systems to break the time-reversal
symmetry. Based on Floquet theory, Chen et al. proposed a topological quantum
phase transition to a quantum anomalous Hall phase induced by off-resonant
circularly polarized light in a two-dimensional system that is initially in a
quantum spin Hall phase or in a trivial insulator phase [24]. It has been
demonstrated that in monolayer graphene [25] and silicene [26], respectively,
irradiated by an off-resonance circularly polarized light a topological phase
transition between different Hall states takes place. Laser-induced quantum
anomalous Hall states in honeycomb lattices have also been investigated [27].
Weyl semimetals showed a large anomalous Hall effect controllable by
illuminating an off-resonant circularly polarized light [28]. Previous studies
were mostly performed on steady states in systems with isotropic low-energy
dispersions while Hall effects in periodically driven heterostructures are
still lacking.

In this paper, we use the Floquet theory to study the Hall effects in PbC/MnSe
heterostructure [29] with spin orbit coupling under the illumination of an
off-resonant circularly polarized light. PbC monolayer grown on an MnSe (111)
surface has a very large energy gap (about 244 meV), and the structure of the
system is mechanically stable. The half-metallic property owned by the
heterostructure is protected by the topology of the system and is robust
against the variation of the chemical potential and strain. Due to these
prominent merits, PbC/MnSe heterostructure provides a novel platform to study
the Hall effects.

The structure of this paper is organized as follows. In section 2 we present a
brief account of the model and a research method of PbC/MnSe heterostructures
with a circularly polarized light, and discuss the evolution of the energy
gap. Based on the Floquet theory, we study the Chern numbers and Hall phases
in our model at high frequency in section 3. In section 4, we study the
anomalous Hall and spin Hall conductivities. In section 5, we summarize our
results.

##  2 Model and method

We begin with the tight-binding Hamiltonian of PbC/MnSe heterostructure in
momentum space on the basis of ${\mid p_{+,\uparrow}\rangle,\mid
p_{-,\uparrow}\rangle,\mid p_{+,\downarrow}\rangle,\mid
p_{-,\downarrow}\rangle}$ near the Fermi energy around the $\Gamma$ point,
which is given by [29]

$\displaystyle H_{0}$ $\displaystyle=$
$\displaystyle\left({\begin{array}[]{cccc}\varepsilon_{1}(k)&f(k)&0&0\\\
f^{*}(k)&\varepsilon_{1}(k)&0&0\\\ 0&0&\varepsilon_{1}(k)&f(k)\\\
0&0&f^{*}(k)&\varepsilon_{2}(k)\end{array}}\right),$ (2.5)

where $\mid p_{\pm}\rangle=(\mid p_{x}\rangle\pm\mathrm{i}\mid
p_{y}\rangle)/\sqrt{2}$. $\mid p_{x,y}\rangle$ denote the orbital states in
Pb, and $\uparrow\downarrow$ represent the spin-up (-down) state. Taking time-
reversal symmetry and the $C_{3v}$ symmetry into consideration, one can obtain
$f(k_{\pm})=\mathrm{e}^{\mathrm{i}\frac{4\piup}{3}}f(\pm\mathrm{i}\frac{2\piup}{3}k_{\pm})$
with $k_{\pm}=k_{x}+\mathrm{i}k_{y}$. Therefore, $f(k)$ will take the form
$f(k)=\beta k^{2}_{-}$, and $\varepsilon_{1}(k)=\varepsilon_{2}(k)$ is
required. Finally, considering the atomic spin orbit coupling (SOC)
interaction, the $k\cdot p$ Hamiltonian for low-energy physics of the PbC/MnSe
heterostructure is

$\displaystyle\hat{H}$ $\displaystyle=$ $\displaystyle\alpha
k^{2}+\beta\left[(k_{x}^{2}-k_{y}^{2})\hat{\sigma}_{x}+2k_{x}k_{y}\hat{\sigma}_{y}\right]+s\lambda\hat{\sigma}_{z}\,,$
(2.6)

where $\alpha$ and $\beta$ are constant taking the symmetry into
consideration. $\hat{\sigma}_{x,y,z}$ denotes the Pauli matrices for the basis
$\mid p_{\pm}\rangle=(\mid p_{x}\rangle\pm\mathrm{i}\mid
p_{y}\rangle)/\sqrt{2}$. $\lambda$ denotes SOC interaction for the spin up
($s=1$) and spin down ($s=-1$) subspace, respectively [30]. For simplicity, we
take $\beta$ as energy unit and set $\beta=1$ throughout the paper.

Figure 1:  (Colour online) Bulk band gap as a function of the light induced
term $\lambda_{\omega}=\Delta_{\omega}|_{k=0}$ when (a) $\kappa=1$ and (b)
$\kappa=-1$. The blue dashed and red solid lines denote spin-up and spin-down
bands, respectively. $\lambda$ is set to be $0.1$.

Considering the circularly polarized light irradiated onto the PbC/MnSe
heterostructure, the vector potential introduced into the Hamiltonian by
substitution $\vec{k}\rightarrow\vec{k}+e\vec{A}(t)$ is as follows:

$\displaystyle\vec{A}(t)=A(\kappa\sin{\omega t},\cos{\omega t}),$ (2.7)

where $\omega$ is the frequency of light, $\kappa=\pm 1$ denote right-hand and
left-hand polarization of light, respectively. $A$ is the amplitude of the
vector potential. Based on Floquet theory [31, 32, 33], the photon dressed
effective Hamiltonian with weak light intensity can be written as follows:

$\displaystyle\hat{H}^{\text{F}}=\hat{H}+\frac{[\hat{H}_{-1},\hat{H}_{1}]}{\hbar\omega}+\frac{[\hat{H}_{-2},\hat{H}_{2}]}{2\hbar\omega}\,,$
(2.8)

where

$\displaystyle\hat{H}_{\pm 1}$ $\displaystyle=$
$\displaystyle\frac{\omega}{2\piup}\int_{0}^{\frac{2\piup}{\omega}}\hat{H}(t)\mathrm{e}^{\pm\mathrm{i}\omega
t},$ (2.9) $\displaystyle\hat{H}_{\pm 2}$ $\displaystyle=$
$\displaystyle\frac{\omega}{2\piup}\int_{0}^{\frac{2\piup}{\omega}}\hat{H}(t)\mathrm{e}^{\pm
2\mathrm{i}\omega t}.$ (2.10)

$\hat{H}(t)$ is the time-dependent Hamiltonian using the substitution
$\vec{k}\rightarrow\vec{k}+e\vec{A}(t)$ in equation (2.6). Straightforward
calculations show

$\displaystyle\hat{H}^{\text{F}}$ $\displaystyle=$ $\displaystyle\alpha
k^{2}+\big{(}k_{x}^{2}-k_{y}^{2}\big{)}\hat{\sigma}_{x}+2k_{x}k_{y}\hat{\sigma}_{y}+[s\lambda+\kappa\Delta_{\omega}]\hat{\sigma}_{z}\,,$
(2.11)

where

$\displaystyle\Delta_{\omega}=\frac{1}{\hbar\omega}\left(\frac{eA}{\hbar}\right)^{2}\left[4k^{2}+\frac{1}{2}\left(\frac{eA}{\hbar}\right)^{2}\right].$
(2.12)

We note that an extra light induced term ($\Delta_{\omega}$) is introduced,
and its sign depends on $\kappa$ rather than on spin. This new term is clearly
different from the intrinsic SOC term ($s\lambda$) and also different from the
conventional Zeeman term. For $\lambda>0$ and $\Delta_{\omega}<0$, the overall
value of the effective SOC ($s\lambda+\kappa\Delta_{\omega}$) will decrease
(increase) if $s$ and $\eta$ have the same (different) sign. The critical
point $\lambda_{\omega}$ is defined by $s\lambda+\kappa\Delta_{\omega}=0$ at
$k=0$, where the $\lambda_{\omega}$ vanishes and band gap closes.
Consequently, the spin-up and -down bands will have an opposite response to
the circularly polarized light of the opposite $\kappa$. Figure 1 (a) shows
the gap as a function of the light induced term
$\lambda_{\omega}=\Delta_{\omega}|_{k=0}$ for $\kappa=1$ and $\lambda=0.1$.
When increasing $\lambda_{\omega}$, the gap of spin-down component decreases
while that of the spin up component increases. Then, the spin-down gap closes
and reopens. Figure 1 (b) shows that similar changes are found for
$\kappa=-1$. It indicates that a topological phase transition takes place. The
topologically different responses arise from the interplay between the SOC and
light.

##  3 Chern numbers and Hall phase

For our system, the Hamiltonian can always be expressed in terms of the Pauli
matrices as follows:

$\displaystyle\hat{H}^{\text{F}}({\bf k})={\bf d}({\bf
k})\cdot\hat{\sigma}+\epsilon({\bf k})\mathbf{I}\,,$ (3.1)

where ${\bf d}({\bf k})=(d_{x}({\bf k}),d_{y}({\bf k}),d_{z}({\bf k}))$,
$\hat{\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z})$ , and $\mathbf{I}$ is the
rank-2 unit matrix. From the above effective Hamiltonian, the Chern number is
calculated as [34]

$\displaystyle
C(s)=\frac{1}{2\piup}\int_{BZ}\mathrm{d}^{2}\textbf{\emph{k}}\Omega(s)\,,$
(3.2)

where $\Omega(s)$ is the Berry curvature in the momentum space over all
occupied states of electrons with spin $s$ component, and it can be written as
follows:

$\displaystyle\Omega(s)=\frac{1}{2}\hat{d}\cdot\left(\frac{\partial\hat{d}}{\partial
k_{x}}\times\frac{\partial\hat{d}}{\partial k_{y}}\right),$ (3.3)

where $\hat{d}=\frac{{\bf d}}{|{\bf d}|}$. Straightforward calculation shows

$\displaystyle\Omega(s)=\frac{2k^{2}\left[s\lambda+\kappa\frac{1}{2\hbar\omega}\left(\frac{eA}{\hbar}\right)^{4}\right]}{\left[k^{4}+(s\lambda+\kappa\Delta_{\omega})^{2}\right]^{3/2}},$
(3.4)

and

$\displaystyle C(s)=\operatorname{sgn}(s\lambda+\kappa\lambda_{\omega})\,,$
(3.5)

where $\operatorname{sgn}$ is the sign function. One can see that the Chern
number will change its sign with the varying of
$\operatorname{sgn}(s\lambda+\kappa\lambda_{\omega})$.

For $\alpha=\lambda=0$, the Hamiltonian [equation (2.11)] is similar to the
bilayer graphene irradiated with circularly polarized light [35, 36, 37, 38].
In the bilayer graphene irradiated with circularly polarized light, the light
induces a term which is valley and polarization dependent. Therefore, the
Berry curvatures are the valley and polarization dependent, and the
integration gives an integer nonzero Chern number per valley, hence realizing
a quantum valley-Hall state. In our system, a quantum spin Hall state probably
exists due to SOC.

Figure 2:  (Colour online)  Chern numbers as a function of the light induced
term $\lambda_{\omega}=\Delta_{\omega}|_{k=0}$ ${C}_{s}$ (dashed blue line)
and $C$ (solid red line) for $\lambda=0.1$ and $\kappa=1$.

In figure 2, we plot Chern numbers as a function of
$\lambda_{\omega}=\Delta_{\omega}|_{k=0}$ for $\lambda=0.1$ and $\kappa=1$.
When the Fermi level lies inside the bulk energy gap, it can be found that the
quantum state is a quantum spin Hall (QSH) insulator when there is no
circularly polarized light is absent, which is characterized by a nonzero spin
Chern number ${C}_{s}=\frac{1}{2}[C(1)-C(-1)]$. We have considered that the
spin $s_{z}$ is a good quantum number. With an increasing circularly polarized
light intensity, the spin Chern number ${C}_{s}$ changes from
$\operatorname{sgn}(\lambda)$ to $0$, while the Chern number $C=C(1)+C(-1)$
changes from $0$ to $2\kappa$ with closure and reopening of the gap. The
nonzero Chern number $C$ means that the insulating state is a quantum
anomalous Hall (QAH) insulator. Therefore, a QSH to QAH phase transition takes
place, which is consistent with our numerical results that the gap closes and
reopens. It should be noted that the topological phase transition results from
emission and absorption of two virtual photons (see equation 2.12) in this
paper. In the next section we will study the anomalous Hall and spin Hall
conductivities to confirm our results. Thus, there is a direct correspondence
between the Chern number and the Hall conductance for two-dimensional
insulators ($\sigma_{xy}=Ce^{2}/h$) when the Fermi level lies inside the bulk
energy gap.

It should be noted that the first term in equation (2.11) in the revised
manuscript does not modify the eigenstates of our system. Therefore, we will
set $\alpha=0$ to explore the Hall effects in the remainder of this paper.
Furthermore, the critical point is
${e}A/{\hbar}=(2s\kappa\lambda\omega)^{\frac{1}{4}}$ when the topological
phase transition take place based on the Hamiltonian (2.11). The atomic spin
orbit is small which is responsible for the opening of the gap. At high
frequency regime, we have set $\omega=6$. Thus the light intensity
${e}A/{\hbar}$ can be tuned to be small, which can be realized experimentally.
On the other hand, the number of the edge state is characterized by Chern
number, and topological phase transition is associated with the closure and
reopening of the energy of the gap. Thus, we can observe the closure and
reopening of the energy, the number of each edge, and propagation direction of
the each edge to confirm the topological phase transition experimentally.

##  4 Anomalous Hall and spin Hall conductivities

To better understand the topological properties, we studied the anomalous Hall
and Spin Hall conductivities in this section. Based on Green’s function theory
[39, 40], anomalous Hall and spin Hall can be expressed as follows:

$\displaystyle\sigma_{xy}$ $\displaystyle=$
$\displaystyle\sigma_{xy}^{c\text{I}}+\sigma_{xy}^{c\text{II}},$ (4.1)
$\displaystyle\sigma_{xy}^{s_{z}}$ $\displaystyle=$
$\displaystyle\sigma_{xy}^{s\text{I}}+\sigma_{xy}^{s\text{II}}.$ (4.2)

The first part is the contribution from the Fermi surface, the other part is
the contribution from the Fermi sea, and

$\displaystyle\sigma_{xy}^{\text{c(s)\text{I}}}$ $\displaystyle=$
$\displaystyle-\frac{e\hbar}{2}\int\frac{\mathrm{d}\varepsilon}{2\piup}\partial_{\varepsilon}f(\varepsilon)\int\frac{\mathrm{d}^{2}\textbf{k}}{(2\piup\hbar)^{2}}{\rm
Tr}\left[\hat{j}_{\text{c(s)}}(\textbf{k})\hat{G}_{0}^{\text{R}}(\varepsilon,\textbf{k})\hat{v}_{y}(\textbf{k})\big{(}\hat{G}_{0}^{\text{A}}(\varepsilon,\textbf{k})-\hat{G}_{0}^{\text{R}}(\varepsilon,\textbf{k})\big{)}\right.$
(4.3)
$\displaystyle\left.-\hat{j}_{\text{c(s)}}(\textbf{k})\big{(}\hat{G}_{0}^{\text{A}}(\varepsilon,\textbf{k})-\hat{G}_{0}^{\text{R}}(\varepsilon,\textbf{k})\big{)}\hat{v}_{y}(\textbf{k})\hat{G}_{0}^{\text{R}}(\varepsilon,\textbf{k})\right],$
$\displaystyle\sigma_{xy}^{\text{c(s)\text{II}}}$ $\displaystyle=$
$\displaystyle-\frac{e\hbar}{2}\int\frac{\mathrm{d}\varepsilon}{2\piup}f(\varepsilon)\int\frac{\mathrm{d}^{2}\textbf{k}}{(2\piup\hbar)^{2}}{\rm
Tr}\left[\hat{j}_{\text{c(s)}}(\textbf{k})\hat{G}_{0}^{\text{A}}(\varepsilon,\textbf{k})\hat{G}_{0}^{\text{A}}(\varepsilon,\textbf{k})\hat{v}_{y}(\textbf{k})\hat{G}_{0}^{\text{A}}(\varepsilon,\textbf{k})\right.$
(4.4)
$\displaystyle\left.-\hat{j}_{\text{c(s)}}(\textbf{k})\hat{G}_{0}^{\text{A}}(\varepsilon,\textbf{k})\hat{v}_{y}(\textbf{k})\hat{G}_{0}^{\text{A}}(\varepsilon,\textbf{k})\hat{G}_{0}^{\text{A}}(\varepsilon,\textbf{k})-\big{(}\hat{G}_{0}^{\text{A}}\rightarrow\hat{G}_{0}^{\text{R}}\big{)}\right],$

where $f(\varepsilon)$ is the Fermi distribution.
$\hat{j_{\text{c}}}=\hat{v}_{x}$ and
$\hat{j_{\text{s}}}=\frac{1}{2}\\{\hat{v}_{x},\hat{s}_{z}\\}$ are the charge
and spin current density operator, respectively.
$\hat{v}_{x,y}=(1/\hbar)(\partial H/\partial p_{x,y})$ is the velocity
operator.

$\hat{G}^{\text{R,A}}_{0}(\varepsilon,\textbf{k})=[\varepsilon-\hat{H}^{\text{F}}(\textbf{k})\pm\mathrm{i}0^{+}]^{-1}$
(4.5)

is the retarded (advanced) Green’s function in the clean limit.

Figure 3:  (Colour online) Spin Hall and Hall conductivity in units of
$e/2\piup$ and $e^{2}/2\piup\hbar$, respectively, under the irradiation of
circularly polarized light. The red and blue solid lines denote the right-hand
and left-hand circularly polarized light, respectively. The typical parameters
$\lambda_{\omega}$ are set to $0.05$ in (a) (c) and $0.15$ in (b) (d).

Figure 3 shows numerical results of the spin Hall ($\sigma^{s_{z}}_{xy}$) and
anomalous Hall conductivities ($\sigma_{xy}$) for both the left-hand and
right-hand circularly polarized lights. As shown in figures 3 (a) and (c), the
system is QSH insulator with quantized spin Hall conductivity of $e/2\piup$
when the Fermi level lies inside the energy gap for both circularly polarized
lights before the closure of the energy gap. However, after reopening of the
energy gap, the spin Hall conductivity is quenched to zero. This is because
the Chern number for both two-spin components is determined by the sign of
$s\lambda+\kappa\lambda_{\omega}$ (see equation 3.5). We can also find that
the Chern number $C=2$ and $C=-2$ for the QAH state under the left-hand and
right-hand circularly polarized lights, respectively, after reopening the
energy gap. This is consistent with the quantized anomalous Hall conductivity
of $e^{2}/2h$ and $-e^{2}/2h$ for the two circularly polarized lights,
respectively, as shown in figures 3 (b) and (d). Consequently, a QSH to QAH
topological phase transition is realized, and it depends on the polarization
of light.

On the other hand, $\sigma^{s_{z}}_{xy}$ [figures 3 (a) and (c)] drops
promptly from the quantized plateau. This behavior is independent of
polarization of light. When the Fermi energy ($|E_{\text{F}}|$) is outside the
energy gap, the $\sigma^{s_{z}}_{xy}$ decreases with increasing
$|E_{\text{F}}|$. The $\sigma_{xy}$ [figures 3 (b) and (d)] is drops/raises
promptly from the quantized plateau for the right-hand/left-hand circularly
polarized light. When $|E_{\text{F}}|$ is outside the energy gap, the
$|\sigma_{xy}|$ decreases with increasing $|E_{\text{F}}|$. There is a link in
$\sigma_{xy}$. The nonzero Berry curvature of the band mainly distributes near
its band edge, leading to this behavior of $\sigma^{s_{z}}_{xy}$ and
$\sigma_{xy}$. Based on equation (3.4), $\sigma^{s_{z}}_{xy}$ and
$\sigma_{xy}$ will always vanish at a higher energy.

##  5 Conclusion and summary

To conclude, we have analyzed a topological phase transition and Hall
conductivities in PbC/MnSe heterostructure under the illumination of a
circularly polarized light. We have employed Floquet theory and Green’s
function formalism to compute Chern numbers and Hall conductivities at high-
frequency regime, respectively. The Chern numbers are determined by the
interplay between SOC and light induced SOC-like term. A photo-induced
transition between quantum spin Hall state and anomalous Hall state was
observed. When the Fermi energy lies in the energy gap, the anomalous Hall
conductivities are dependent on polarization of light, while spin Hall
conductivity is independent. Our findings provide a way to control the
topological phase transition and transport properties of the heterostructure.
This is also important for studying the effect of spin-orbit coupling in a
heterostructure.

## References

  * [1] Ren Y., Qiao Z., Niu Q., Rep. Prog. Phys., 2016, 79, 066501, doi:10.1088/0034-4885/79/6/066501.
  * [2] Zhai H., Rechtsman M., Lu Y.M., Yang K., New J. Phys., 2016, 18, 080201, doi:10.1088/1367-2630/18/8/080201.
  * [3] Nagaosa N., Sinova J., Onoda S., MacDonald A.H., Ong N.P., Rev. Mod. Phys., 2010, 82, 1539,
doi:10.1103/RevModPhys.82.1539.

  * [4] Sinova J., Valenzuela S.O., Wunderlich J., Back C.H., Jungwirth T., Rev. Mod. Phys., 2015,  87 1213,
doi:10.1103/RevModPhys.87.1213.

  * [5] Liu C.X., Qi X.L., Dai X., Fang Z., Zhang S.C., Phys. Rev. Lett., 2008, 101, 146802,
doi:10.1103/PhysRevLett.101.146802.

  * [6] Yu R., Zhang W., Zhang H.-J., Zhang S.-C., Dai X., Fang Z., Science, 2010, 329, 61,
doi:10.1126/science.1187485.

  * [7] Chang C.-Z., Zhang J.S., Feng X., Shen J., Zhang Z.C., Guo M.H., Li K., Ou Y.B., Wei P., Wang L.L. et al., Science, 2013, 340, 167, doi:10.1126/science.1234414.
  * [8] Kou X.F., Guo S.-T., Fan Y.B., Pan L., Lang M.R., Jiang Y., Shao Q.M., Nie T.X., Murata K., Tang J.S., Wang Y., He L., Lee T.-K., Lee W.-L., Wang K.L., Phys. Rev. Lett., 2014, 113, 137201,
doi:10.1103/PhysRevLett.113.137201.

  * [9] Bernevig B.A., Hughes T.L., Zhang S.-C., Science, 2006, 314, 1757, doi:10.1126/science.1133734.
  * [10] König M., Wiedmann S., Brüne C., Roth A., Buhmann H., Molenkamp L., Qi X.L., Zhang S.C., Science, 2007, 318, 766, doi:10.1126/science.1148047.
  * [11] Li J., Chu R.-L., Jain J.K., Shen S.-Q., Phys. Rev. Lett., 2009, 102, 136806,
doi:10.1103/PhysRevLett.102.136806.

  * [12] Kane C.L., Mele E.J., Phys. Rev. Lett., 2005, 95, 146802, doi:10.1103/PhysRevLett.95.146802.
  * [13] Kane C.L., Mele E.J., Phys. Rev. Lett., 2005, 95, 226801, doi:10.1103/PhysRevLett.95.226801.
  * [14] Cao J., Xiong S.-J., Phys. Rev. B, 2013, 88, 085409, doi:10.1103/PhysRevB.88.085409.
  * [15] Chen Z.-G., Mei J., Sun X.-C., Zhang X., Zhao J., Wu Y., Phys. Rev. A, 2017, 95, 043827,
doi:10.1103/PhysRevA.95.043827.

  * [16] Liu F., Deng H.-Y., Wakabayashi K., Phys. Rev. B, 2018, 97, 035442, doi:10.1103/PhysRevB.97.035442.
  * [17] Kort-Kamp W.J.M., Phys. Rev. Lett., 2017, 119, 147401, doi:10.1103/PhysRevLett.119.147401.
  * [18] Oshima T., Kawaguchi Y., Phys. Rev. A, 2016, 93, 053605, doi:10.1103/PhysRevA.93.053605.
  * [19] Schweizer C., Lohse M., Citro R., Bloch I., Phys. Rev. Lett., 2016, 117, 170405,
doi:10.1103/PhysRevLett.117.170405.

  * [20] Armaitis J., Ruseckas J., Juzeliũnas G., Phys. Rev. A, 2017, 95, 033635, doi:10.1103/PhysRevA.95.033635.
  * [21] Sun X.-Q., Zhang S.-C., Wang Z., Phys. Rev. Lett., 2015, 115, 076802, doi:10.1103/PhysRevLett.115.076802.
  * [22] Castro E.V., Grushin A.G., Valenzuela B., Vozmediano M.A.H., Cortijo A., de Juan F., Phys. Rev. Lett., 2011, 107, 106402, doi:10.1103/PhysRevLett.107.106402.
  * [23] Sun K., Liu W.V., Hemmerich A., Das Sarma S., Nat. Phys., 2012, 8, 67, doi:10.1038/nphys2134.
  * [24] Chen M.N., Su W., Deng M.X., Ruan J., Luo W., Shao D.X., Sheng L., Xing D.Y., Phys. Rev. B, 2016, 94, 205429, doi:10.1103/PhysRevB.94.205429.
  * [25] Zhai X., Jin G., Phys. Rev. B, 2014, 89, 235416, doi:10.1103/PhysRevB.89.235416.
  * [26] Ezawa M., Phys. Rev. Lett., 2013, 110, 026603, doi:10.1103/PhysRevLett.110.026603.
  * [27] Gómez-León Á., Delplace P., Platero G., Phys. Rev. B, 2014, 89, 205408, doi:10.1103/PhysRevB.89.205408.
  * [28] Yan Z., Wang Z., Phys. Rev. Lett., 2016, 117, 087402, doi:10.1103/PhysRevLett.117.087402.
  * [29] Xue Y., Zhang J.Y., Zhao B., Wei X.Y., Yang Z.Q., Nanoscale, 2018, 10, 8569, doi:10.1039/C8NR00201K.
  * [30] Liang Q.-F., Yu R., Zhou J., Hu X., Phys. Rev. B, 2016, 93, 035135, doi:10.1103/PhysRevB.93.035135.
  * [31] Oka T., Aoki H., Phys. Rev. B, 2009, 79, 081406, doi:10.1103/PhysRevB.79.081406.
  * [32] Kitagawa T., Oka T., Brataas A., Fu L., Demler E., Phys. Rev. B, 2011, 84, 235108,
doi:10.1103/PhysRevB.84.235108.

  * [33] Lindner N.H., Refael G., Galitski V., Nat. Phys., 2011, 7, 490, doi:10.1038/nphys1926.
  * [34] Qi X.-L., Hughes T.L., Zhang S.-C., Phys. Rev. B, 2008, 78, 195424, doi:10.1103/PhysRevB.78.195424.
  * [35] Morell E.S., Foa Torres Luis E.F., 2012, 86, 125449, doi:10.1103/PhysRevB.86.125449.
  * [36] Abergela D.S.L., Chakraborty T., 2009, 95, 062107, doi:10.1063/1.3205117.
  * [37] Iorsh I.V., Dini K., Kibis O.V., Shelykh I.A., 2017, 96, 155432, doi:10.1103/PhysRevB.96.155432.
  * [38] Qu C., Zhang C., Zhang F., 2018, 5, 011005, doi:10.1088/2053-1583/aa9471.
  * [39] Onoda S., Sugimoto N., Nagaosa N., Phys. Rev. Lett., 2006, 97, 126602, doi:10.1103/PhysRevLett.97.126602.
  * [40] Kovalev A.A., Sinova J., Tserkovnyak Y., Phys. Rev. Lett., 2010, 105, 036601,
doi:10.1103/PhysRevLett.105.036601.

Ukrainian 0

Êðîñîâåð ñïí-õîëëâñüêîãî åôåêòó äî êâàíòîâîãî àíîìàëüíîãî åôåêòó Õîëëà â
PbC/MnSe ãåòåðîñòðóêòóðàõ Ë. åí

Âèùà øêîëà ôçèêè òà åëåêòðîíêè, óíâåðñèòåò Òàéøàíü, 271000 Øàíüäóí, Êèòàé