# Effects of metallic electrodes on the thermoelectric properties of zigzag
graphene nanoribbons with periodic vacancies

David M T Kuo  Department of Electrical Engineering and Department of Physics,
National Central University, Chungli, 320 Taiwan

###### Abstract

We theoretically analyze the thermoelectric properties of graphene quantum dot
arrays (GQDAs) with line- or surface-contacted metal electrodes. Such GQDAs
are realized as zigzag graphene nanoribbons (ZGNRs) with periodic vacancies.
Gaps and minibands are formed in these GQDAs, which can have metallic and
semiconducting phases. The electronic states of the first conduction (valence)
miniband with nonlinear dispersion may have long coherent lengths along the
zigzag edge direction. With line-contacted metal electrodes, the GQDAs have
the characteristics of serially coupled quantum dots (SCQDs) if the armchair
edge atoms of the ZGNRs are coupled to the electrodes. By contrast, the GQDAs
have the characteristics of parallel QDs if the zigzag edge atoms are coupled
to the electrodes. The maximum thermoelectric power factors of SCQDs with
line-contacted electrodes of Cu, Au, Pt, Pd, or Ti at room temperature were
similar or greater than $0.186nW/K$; their figures of merit were greater than
three. GQDAs with line-contacted metal electrodes have much better
thermoelectric performance than surface contacted metal electrodes.

##  I Introduction

Charge transport through individual quantum dots with discrete levels exhibits
an interesting phenomenon such as the Kondo effect and Coulomb
blockade.[HaugH, 1] Quantum dots (QDs) have also been suggested as
implementations of low-power devices. These low-power devices include single-
electron transistors[GuoLJ, 2,Postma, 3], single-photon sources[Michler,
4-ChangWH, 6], single-photon detectors[Gustavsson, 7] and single-electron heat
engines[Josefsson, 8]. Some applications of QD devices require both high
efficiency and considerable output power. These applications require QD solids
that can retain the size-tunable properties of QDs and possess the band
transport characteristic of bulk semiconductors [Kagan, 9]. Although
substantial effort has been devoted to producing such QD solids, experimental
studies of the thermoelectric properties of such one-dimensional QD arrays are
lacking [Harman, 10,Talgorn, 11]. Dot-size fluctuation is a challenge for
thermoelectric devices made from 1D silicon QD arrays (QDAs) [Lawrie, 12].

Recently, numerous studies have focused on two-dimensional (2D) materials such
as graphenes [Novoselovks, 13], and 2D transition metal dichalcogenides (TMDs)
and oxides [Geim, 14–Desai, 16]. Quasi-1D graphene nanoribbons (GNRs) can now
be fabricated with atomic precision by using the bottom-up approach [Cai, 17];
this approach can yield complex graphene-based nanostructures on-demand for
quantum device applications [LiuJ, 18–Almeida, 28]. The superlattices of
zigzag GNRs (ZGNRs) behave as graphene QDAs (GQDAs) [Topsakal, 29]. The phonon
thermal conductance of GQDAs can be dramatically reduced by at least one order
of magnitude compared with that of ZGNRs [Cuniberti, 30]. Therefore, GQDAs are
expected to have promising applications in nanoscale energy harvesting.

Nevertheless, previous studies have demonstrated that the contact types of
grapehene and GNRs significantly influence electron transport between the
electrodes [Darancet, 31–Matsuda, 33]. To reduce the number of electron
backward scattering components, GNR devices can be produced using normal metal
electrodes instead of graphene electrodes to improve the transmission
coefficient [GLiang, 32]. In addition, the contact resistance of a graphene
with line-contacted metal electrodes can be much smaller than that of a
graphene with surface-contacted metal electrodes [Matsuda, 33]. Other 2D
materials also have this feature [Shen, 34, ChenRS, 35]. These results all
indicate that the interface properties arising from the presence of a contact
junction is a key factor affecting the implementation of graphene-based or
TMD-based electronic and thermoelectric devices.

Few studies have examined how the performance of GQDA thermoelectric devices
is affected by the material properties and geometries of various contacting
metals. In this study, we theoretically investigate the ballistic transport
and thermoelectric properties of GQDAs coupled to various metal electrodes, as
shown in Fig. 1; GQDAs are realized by ZGNRs with periodic vacancies. We have
focused on the channel length between thermal contact shorter than the
electron mean free path ($\lambda_{e}$) but longer than the phonon mean free
path ($\lambda_{ph}$). We calculated the electrical conductance ($G_{e}$),
Seebeck coefficient ($S$) and electron thermal conductance ($\kappa_{e}$) for
different line-contacted metal electrodes and surface-contacted metal
electrodes using the tight-binding model and Green’s function. We find that
the maximum thermoelectric power factor values of GQDAs with Cu, Au, Pt, Pd or
Ti line-contacted electrodes can reach up to $79\%$ of the power factor of the
1D theoretical limit at room temperature.

Figure 1: Diagram of ZGNR without and with vacancies for $N_{a}=8$ and
$N_{z}=35$. (a) Armchair edge atoms of ZGNRs coupled to metallic electrodes.
$\Gamma_{L}$ ($\Gamma_{R}$) denotes the tunneling rate of the electrons
between the left (right) electrode and the leftmost (rightmost) carbon atoms
of the armchair edges. $T$ denotes the equilibrium temperature of the left
(right) electrode. (b, c) Graphene quantum dot arrays (GQDAs) are realized by
periodically removing one carbon and two carbons from the interior sites of
ZGNRs. The distance between the nearest vacancies is $D=4a$, where the lattice
constant of graphene is $a=2.46\AA$. (d) The zigzag edge atoms of GQDAs are
coupled to the metal electrodes. (e) GQDAs with surface-contacted metal
electrodes.

##  II Calculation method

A good approximation of the electronic states of a GQDA realized by ZGNRs with
periodic vacancies is a tight-binding model with one $p_{z}$ orbital per
atomic site [Topsakal, 29, Wakabayashi2, 36]. The Hamiltonian of the GQDA,
written, $H_{GQDA}$, is

$\displaystyle H_{GQDA}$ $\displaystyle=$
$\displaystyle\sum_{\ell,j}E_{\ell,j}d^{\dagger}_{\ell,j}d_{\ell,j}$
$\displaystyle-$
$\displaystyle\sum_{\ell,j}\sum_{\ell^{\prime},j^{\prime}}t_{(\ell,j),(\ell^{\prime},j^{\prime})}d^{\dagger}_{\ell,j}d_{\ell^{\prime},j^{\prime}}+h.c,$

where $E_{\ell,j}$ is the on-site energy for the $p_{z}$ orbital in the
${\ell}$th row and $j$th column. $d^{\dagger}_{\ell,j}(d_{\ell,j})$ creates
(destroys) one electron at the atom site labeled by ($\ell$,$j$) where $\ell$
and $j$ are the row and column indices, respectively, as depicted in Fig. 1.
$t_{(\ell,j),(\ell^{\prime},j^{\prime})}$ describes the electron hopping
energy from site ($\ell$,$j$) to site ($\ell^{\prime}$,$j^{\prime}$). The
tight-binding parameters used for GQDAs are $E_{\ell,j}=0$ for the on-site
energy and $t_{(\ell,j),(\ell^{\prime},j^{\prime})}=t_{pp\pi}=2.7$ eV for the
nearest-neighbor hopping strength. Here, the electron–electron Coulomb
interactions are neglected. Their effect is discussed in a subsequent section.

Thermoelectric coefficients, namely the electrical conductance ($G_{e}$),
Seebeck coefficient ($S$) and electron thermal conductance ($\kappa_{e}$) were
calculated using $G_{e}=e^{2}{\cal L}_{0}$, $S=-{\cal L}_{1}/(eT{\cal L}_{0})$
and $\kappa_{e}=\frac{1}{T}({\cal L}_{2}-{\cal L}^{2}_{1}/{\cal L}_{0})$ with
${\cal L}_{n}$ ($n=0,1,2$) defined as

${\cal L}_{n}=\frac{2}{h}\int d\varepsilon~{}{\cal
T}_{LR}(\varepsilon)(\varepsilon-\mu)^{n}\frac{\partial
f(\varepsilon)}{\partial\mu}.$ (2)

Here, $f(\varepsilon)=1/(exp^{(\varepsilon-\mu)/k_{B}T}+1)$ is the Fermi
distribution function of electrodes at equilibrium temperature $T$ and
chemical potential $\mu$. ${\cal T}_{LR}(\varepsilon)$ denotes the
transmission coefficient of a GQDA connected to electrodes, which can be
calculated using the formula ${\cal
T}_{LR}(\varepsilon)=4Tr[\Gamma_{L}(\varepsilon)G^{r}(\varepsilon)\Gamma_{R}(\varepsilon)G^{a}(\varepsilon)]$
[KuoDM, 37–KuoD, 39], where ($\Gamma_{L}(\varepsilon)$ and
$\Gamma_{R}(\varepsilon)$) denote the tunneling rate at the left and right
leads, respectively, and ${G}^{r}(\varepsilon)$ and ${G}^{a}(\varepsilon)$ are
the retarded and advanced Green’s functions of the GQDA, respectively. The
tunneling rates are described by the imaginary part of the self energy
resulting from the coupling between the left (right) electrode with its
adjacent GQDA atoms. Based on the tight-binding orbitals,
$\Gamma_{\alpha}(\varepsilon)$ and Green’s functions are matrices. For
simplicity, $\Gamma_{\alpha}$ for interface carbon atoms has diagonal entries
given by the same constant $\Gamma_{t}$. The magnitudes of $\Gamma_{t}$ depend
on the line- or surface-contacted metal electrodes [Matsuda, 33]. The
thermoelectric efficiency of graphene quantum dot arrays (GQDAs) is determined
by the dimensionless figure of merit
$ZT=S^{2}G_{e}T/(\kappa_{e}+\kappa_{ph})$, where $\kappa_{ph}$ denotes the
phonon thermal conductance of GQDAs. First-principle calculations of
$\kappa_{ph}$ are beyond the scope of this article, so for simplicity, we
calculate $\kappa_{ph}$ using an empirical method described in ref. [KuoD,
39]. The phonon mean free path $\lambda_{ph}$ has been measured to be reduced
from $300-600$ nm in a single-layer graphene to $10$ nm in graphene
nanostructures (see [XuY, 40] and references therein). In this study, we chose
zigzag graphene nanoribbons (ZGNRs) with $N_{z}=127$ ($L_{z}=15.5$ nm) because
the channel length between thermal contacts satisfies the condition of being
shorter than $\lambda_{e}$ but longer than $\lambda_{ph}$.

##  III Results and discussion

###  III.1 Electronic energy levels

Graphene with gapless characteristics has limited applications in devices
[Nakada, 41, Wakabayashi, 42]. However, graphene with periodic vacancies has
minibands, which enable its applications in thermoelectric devices at low
temperatures [Gunst, 43,ChangPH, 44]. To enable the graphene-based devices
operating at room temperature, some studies have considered ZGNRs with
periodic vacancies [ZhangYT, 45]. However, few studies have investigated how
line- and surface-contacted metal electrodes affect the thermoelectric
properties of ZGNRs with vacancies. The electron band structures of infinite
ZGNRs without and with periodic vacancies are presented in Fig. 2. The band
structures of infinite ZGNRs with various widths ($N_{a}$) have been
extensively studied [Nakada, 41]. The localized zero-energy flat-band modes of
the first subband are from $k=\frac{2\pi}{3a}$ to $\frac{\pi}{a}$ as
$N_{a}\rightarrow\infty$. In Fig. 2(a), the states of the first subband are
the localized states for $k$ within $k_{c}=0.738\frac{\pi}{a}\leq
k\leq\frac{\pi}{a}$, where $k_{c}$ is determined by
$k_{c}/2=\textrm{arccos}(N_{a}/(2N_{a}+4))$ and $N_{a}=8$ [Wakabayashi2, 36].
The edge state with the shortest decay length along the armchair direction
occurs at $k=\frac{\pi}{a}$. When $k$ deviates from $\frac{\pi}{a}$, the
zigzag edges states form bonding and antibonding states, and the zero-energy
modes disappear [Fig. 2(a)] . When $|E|\geq 2.16$ eV, we observe the second
conduction and valence subbands. Because these are bulk states, this study
focuses on the edge states of the first subband with $|E|\leq 2.16$ eV. For
ZGNRs with periodic vacancies [Fig. 2(b) and Fig. 2(c)], the gapped energy
regions are open and minibands are formed. ZGNR with periodic vacancies can be
in the semiconducting phase as illustrated by the sizable central gap (0.959
eV) in Fig. 2(c). To reflect this phase difference, they are denoted metal
GQDAs (M-GQDAs) and semiconductor GQDAs (S-GQDAs). Because the finite channel
length is shorter than $\lambda_{e}$ but longer than $\lambda_{ph}$, the
length-dependent energy levels of GQDAs should be further clarified.

Figure 2: Electronic band structures of infinite ZGNRs (a) without and (b, c)
with periodic vacancies. The distance between the nearest vacancies is always
$L=4a$ in this article. The size of a super unit cell is characterized by
$N_{a}=8$ and $N_{z}=8$.

To illustrate the characteristics of the energy levels in each miniband of
Fig. 2(b) and 2(c), we present the calculated eigenvalues as functions of QD
number $N_{QD}$ (or $N_{z}=8N_{QD}-1$) in Fig. 3. The spectra reveal the
transformation from the molecule-like state to band-like state as $N_{QD}$
increases from 3 toward infinity (Fig. 2). In Fig. 3, the number of energy
levels in each miniband is given by $N_{QD}$, demonstrating that finite-size
ZGNRs with periodic vacancies exhibit the characteristics of GQDAs. A single
GQD can be characterized by the size of $N_{a}=8$ and $N_{z}=7$, which results
in four energy levels of $\varepsilon_{h,2}=-1.9024$ eV,
$\varepsilon_{h,1}=-0.47905$ eV, $\varepsilon_{e,1}=0.47905$ eV and
$\varepsilon_{e,2}=1.9024$ eV in the range of $|E|\leq 2.16$ eV. These four
energy levels do not vary with $N_{QD}$. To understand this behavior, we
calculate the eigenvalues of ZGNRs without periodic vacancies as a function of
$N_{z}=8N_{QD}-1$ for $N_{a}=8$ and obtain $N_{z}$-independent
$\varepsilon_{e(h),1}$ and $\varepsilon_{e(h),2}$. The wave functions of
$\varepsilon_{e(h),1}$ and $\varepsilon_{e(h),2}$ of the ZGNRs are vanishingly
small at $j=8m$ (so-called nodes), where $m$ is an integer. This explains the
presence of $N_{QD}$-independent $\varepsilon_{e(h),1}$ and
$\varepsilon_{e(h),2}$ in Fig. 3. The sinusoidal wave of
$\varepsilon_{e(h),1}$ ($\varepsilon_{e(h),2}$) along the zigzag edges has a
much longer coherence length than away from these edges. Although the periodic
vacancies do not affect the energy levels of $\varepsilon_{e(h),1}$ and
$\varepsilon_{e(h),2}$, they restrict and perturb some energy levels of ZGNRs.
Therefore, they are called antidots [Gunst, 43–ZhangYT, 45].

Figure 3: Spectra of (a) M-GQDA and (b) S-GQDA as functions of $N_{QD}$. Each
GQD in the GQDA structures has size $N_{a}=8$ and $N_{z}=7$. The length of the
GQDA is $N_{z}=8N_{QD}-1$.

###  III.2  Line-contacted metal electrodes

When a graphene is coupled to metal electrodes, contact properties such as the
Schottky barrier or ohmic contact and contact geometries can significantly
affect the electron transport in graphene [Matsuda, 33, ChuT, 46–Hancheng,
57]. Although many theoretical studies have attempted to clarify this key
behavior from first principles, only qualitative results regarding
$\Gamma_{t}$ arising from the contact junction can be obtained due to
theoretical limitations [YangL, 58, LeeG, 59]. In Fig. 4, we present that the
calculated transmission coefficient of interface $T_{c}(\varepsilon)$ based on
the configuration considered in [Matsuda, 33], in which the authors introduced
“super unit cells ” with eight unit cells and investigated the armchair edge
atoms of graphene coupled to various metal electrodes. Contact resistance
$R_{c}=1/G_{e,c}$ per unit cell is determined by
$G_{e,c}=G_{0}T_{c}(\varepsilon=0)$, where
$G_{0}=2e^{2}/h=1/(12.9k\Omega)=77.5\mu S$ is the quantum conductance and
$T_{c}(\varepsilon=0)$ is the contact transmission coefficient at the Fermi
energy. In [Matsuda, 33], the authors reported that
$R_{c}=13.3,17.8,18.6,23.3$ and $31.7$ $k\Omega$ for Ti, Pd, Pt, Au and Cu,
respectively. The $\Gamma_{t}$ of each carbon atom at the interface can be
obtained by fitting these contact resistances. We obtain the $\Gamma_{t}$
values of Cu, Au, Pt, Pd, and Ti as $\Gamma_{t}=0.495$ eV, $\Gamma_{t}=0.603$
eV, $\Gamma_{t}=0.693$ eV, $\Gamma_{t}=0.72$ eV, and $\Gamma_{t}=0.9$ eV,
respectively.

Figure 4: Transmission coefficient at the metal–graphene interface per unit
cell for various $\Gamma_{t}$ values. Armchair edge atoms of graphene are
coupled to metal electrodes. A super unit cell comprises 16 carbon atoms.

####  III.2.1 Armchair edge atoms coupled to electrodes

Based on the aforementioned $\Gamma_{t}$ values, Fig. 5 presents the
calculated ${\cal T}_{LR}(\varepsilon)$ of M-GQDA with $N_{a}=8$ and
$N_{z}=127$ for various metal electrodes. The maximum ${\cal
T}_{LR}(\varepsilon)=1$ indicates that GQDAs behave as serially coupled QDs
(SCQDs). Note that two energy levels $\varepsilon_{e(h)}=\pm 0.47905eV$ and
$\varepsilon_{e(h)}=\pm 0.47179eV$ are too close to be resolved in the spectra
of ${\cal T}_{LR}(\varepsilon)$. From Eq. (2), the area below the ${\cal
T}_{LR}(\varepsilon)$ curve has a considerable effect on the thermoelectric
coefficients. As illustrated in Fig. 5, the area below the ${\cal
T}_{LR}(\varepsilon)$ curve increases with the tunneling rate. Moreover, these
areas maintain a near-right triangle shape that exhibits a steep change with
respect to $\varepsilon$ on the sides of the $\varepsilon_{e,1}$ and
$\varepsilon_{h,1}$ edges. For energy harvesting applications at room
temperature, we must design a ${\cal T}_{LR}(\varepsilon)$ spectrum with a
square form (SF) to obtain the optimal figure of merit and electrical power
output.[Whitney, 60] However, no method has been developed for realizing a
${\cal T}_{LR}(\varepsilon)$ curve area with SF in a finite channel length
($L_{z}<\lambda_{e}$) between thermal contacts [KuoD, 39].

Figure 5: Transmission coefficient ${\cal T}_{LR}(\varepsilon)$ of a finite
M-GQDA as functions of $\varepsilon$ for various line-contacted metal
electrodes at $N_{a}=8$ ( $L_{a}=0.71$nm ) and $N_{z}=127$ ($L_{z}=15.5$nm ).
Diagrams (a), (b), (c), and (d) correspond to Cu, Au,Pt,and Ti, respectively.

We then calculate the electrical conductance ($G_{e}$), Seebeck coefficient
($S$), power factor ($PF$), and figure of merit ($ZT$) of M-GQDAs at room
temperature (Fig. 6). $G_{e}$, $S$, and $PF$ are in units of $G_{0}=2e^{2}/h$,
$k_{B}/e=86.25\mu V/K$, and $2k^{2}_{B}/h=0.575pW/K^{2}$, respectively. Unlike
$G_{e}$, $S$ is less sensitive to variations in $\Gamma_{t}$ [Fig.6(b)]. The
behavior of $S$ can be approximated by
$S=-\frac{\pi^{2}k^{2}_{B}T}{3e}\frac{\partial
ln(G_{e}(\mu,T))}{\partial\mu}$. The maximum $PF$ and $ZT$ values are given by
$\mu$ in the gap region. The trend of $ZT$ with respect to various metal
electrodes is determined by the power factor ($PF=S^{2}G_{e}$) because the
thermal conductance of GQDAs is dominated by phonon thermal conductance
$\kappa_{ph}$. We calculated the dimensionless
$ZT=S^{2}G_{e}T/(\kappa_{e}+\kappa_{ph})$ with
$\kappa_{ph}=F_{s}\kappa_{ph,ZGNR}$, where
$\kappa_{ph,ZGNR}=\frac{\pi^{2}k^{2}_{B}T}{3h}$ is the phonon thermal
conductance of an ideal ZGNR [Zhengh, 61] and $F_{s}=0.1$ denotes a reduction
factor resulting from periodic vacancies in ZGNRs [Cuniberti, 30, ChangPH,
44]. The results in Fig. 6 imply that $ZT>3$ can be realized if GQDAs are
considered to be coupled to metallic electrodes. The four investigated metal
electrodes have similar maximum $PF$ and $ZT$ values; the maximum $PF$ can
reach $79\%$ of the theoretical limit for 1D systems of
$PF_{QB}=1.2659(\frac{2k^{2}_{B}}{h})$ [Whitney, 60]. ZGNRs are graphene-based
1D topological insulators (TIs), and $ZT$ of 1D TIs greater than 3 has been
theoretically demonstrated in the absence of periodic vacancies [Xu, 62]. The
$ZT$ values of 1D TIs tend to be at their maximum when $\mu$ is near the
minimum value of the second subband, as shown in Fig. 2(a) ($|E|=2.16$ eV).

Figure 6: (a) Electrical conductance $G_{e}$,(b) Seebeck coefficient $S$, (c)
power factor ($PF=S^{2}G_{e}$), and (d) figure of merit $ZT$ as functions of
$\mu$ for various metal electrodes at $k_{B}T=27~{}meV$. Figure 7:
Transmission coefficient ${\cal T}_{LR}(\varepsilon)$ of a finite S-GQDA as
functions of $\varepsilon$ for different line-contacted metal electrodes at
$N_{a}=8$ ( $L_{a}=0.71$nm ) and $N_{z}=127$ ($L_{z}=15.5$nm ). Diagrams (a),
(b), (c), and (d) correspond to Cu, Au, Pt, and Ti, respectively.

In Fig. 6(b), the gaps between the minibands substantially enhance the maximum
$S$ values, suggesting that the thermoelectric properties of S-GQDAs with a
large central band gap should be investigated. The calculated transmission
coefficients of S-GQDAs with various line-contacted metal electrodes are
presented in Fig. 7\. The miniband widths of S-GQDAs are much smaller than
those of M-GQDAs, resulting in a greatly reduced area of their ${\cal
T}_{LR}(\varepsilon)$ curves. Fig. 8 presents the calculated $G_{e}$, $S$,
power factor $PF$, and $ZT$ as functions of $\mu$ at room temperature. Due to
electron-hole symmetry, we considered the range of chemical potential $\mu\geq
0$. The $G_{e}$ spectra of the first miniband remains triangular as in Fig.
6(a), but the steep change with respect to $\mu$ is toward the central gap.
Because of large gaps, the maximum Seebeck coefficients in Fig. 8(b) are
almost twice as large as those in Fig. 6(b). The maximum power factor
$PF_{max}(\mu=0.45eV)=0.9682$ and figure of merit
$ZT_{max}(\mu=0.427eV)=3.442$ occur at the chemical potentials near
$\varepsilon_{e,1}=0.47905eV$ for Cu electrodes. The $\kappa_{ph}$ value used
in the analysis is the same as that in Fig. 6(d). A comparison of the maximum
$ZT$ in Fig. 8(d) and 6(d) reveals that it did not increase despite large
increase in the maximum value of $S$. For $\kappa_{ph}\gg\kappa_{e}$, the
increase of $G_{e}$ also has a substantial effect. If the miniband width was
small, $G_{e}$ was suppressed in the gap region. In Fig. 8(c) and 8(d),
$PF_{max,2}$ and $ZT_{max,2}$ are also smaller than $PF_{max,1}$ and
$ZT_{max,1}$, respectively. These results indicate that increasing electrical
power output requires not only large gaps but also relatively wide miniband
widths.

Figure 8: (a) Electrical conductance $G_{e}$, (b) Seebeck coefficient $S$,(c)
power factor $PF=S^{2}G_{e}$ , and (d) figure of merit ($ZT$) of S-GQDSLs as
functions of $\mu$ for various metal electrodes at $k_{B}T=27meV$.

####  III.2.2 Zigzag edge atoms coupled to electrodes

Our previous study demonstrated that the electron transport of finite ZGNRs is
significantly affected by the armchair or zigzag edge sides coupled to
electrodes [KuoDMT, 38]. According to [LeeG, 59], most edge-oxidized
modifications dramatically change the electron band structures of ZGNRs.
Therefore, we do not consider edge-oxidized GNRs in the following discussion.
The contact resistance of Cr/X/graphene was calculated in [GaoQ, 49], in which
the zigzag edge atoms of graphene were coupled to Cr by various atomic bridges
(H, F, O, or S). Due to the large binding distance, a high contact resistance
was observed in the Cr/H/graphene structure. We adopted $\Gamma_{t}=9$ meV and
$\Gamma_{t}=45$ meV for Cr/H/GQDAs and Cr/F/GQDAs, respectively, to calculate
$G_{e}$ as functions of $\mu$ at zero temperature (Fig. 9) for M-GQDAs and
S-GQDAs with $N_{a}=8$ and $N_{z}=127$. In Fig. 9, the spectra of $G_{e}$ are
larger than quantum conductance $G_{0}=\frac{2e^{2}}{h}$, demonstrating that
GQDAs behave as parallel QDs when the zigzag edge atoms of M-GQDAs and S-GQDAs
are coupled to the metal electrodes. In the parallel QD configuration, the
spectra resulting from electronic states near the zero-energy modes of M-GQDAs
and the high-energy modes of the first minibands of S-GQDAs substantially
differ from the spectra in Fig. 5 and Fig. 7.

Figure 9: Electrical conductance as functions of $\mu$ for $\Gamma_{t}=9$ meV
and $\Gamma_{t}=45$ meV for (a,b) M-GQDAs and (c,d) S-GQDAs. Their sizes are
$N_{a}=8$ ($L_{a}=0.71$ nm) and $N_{z}=127$ ($L_{z}=15.5$ nm).

Fig. 10 presents the calculated $S$ and $PF$ of M-GQDAs and S-GQDAs at room
temperature ($k_{B}T=27$ meV). The maximum $S$ is severely degraded for a
parallel QD configuration with a short channel length of $L_{a}=0.71$ nm
compared with the maximum $S$ values of Fig. 6 and Fig. 8. This is
attributable to the highly enhanced $G_{e}$ in the gap region, which is
similar to the metal-induced gap states in metal-semiconductor Schottky
junctions [Roksana, 63]. When the zigzag edge atoms are coupled to the metal
electrodes, the miniband edge states are readily extended into the gap regions
if the channel length is short. Although the maximum $PF$ in Fig. 10 is much
larger than $PF_{QB}=1.2659(\frac{2k^{2}_{B}}{h})$, their $ZT$ values are less
than one due to remarkable increases in electron thermal conductance
($\kappa_{e}$) and phonon thermal conductance ($\kappa_{ph}$).

Figure 10: (a) Seebeck coefficient and (b) power factor of M-GQDAs and (c,d)
S-GQDAs as functions of $\mu$ for two tunneling rates at $k_{B}T=27$meV. Other
physical parameters are the same as those in Fig. 9.

###  III.3 Surface-contacted metal electrodes

Few experimental studies have reported the transport and thermoelectric
properties of GNRs due to the difficulty of applying the line-contacted
technique [CustiT, 64]. However, surface-contacted metal electrodes are
straightforward to lay out in an experiment. Therefore, experimental methods
are preferable for revealing the thermoelectric properties of GNRs with
surface-contacted electrodes. Many theoretical studies have demonstrated that
the surface contact resistance is one to three orders of magnitude greater
than the line contact resistance [Matsuda, 33, ChuT, 46– Hancheng, 57],
indicating that $\Gamma_{t}$ may be extremely small. The calculated
transmission coefficient of M-GQDA as functions of $\mu$ for various
$\Gamma_{t}$ values is presented in Fig. 11; the two outer GQDs are coupled to
metal electrodes. The structure has 28 carbon atoms, including 6 carbon atoms
each located at zigzag edge sites under the surface-contacted left (right)
metal electrode. Due to the small $\Gamma_{t}$ values, the peaks of the ${\cal
T}_{LR}(\varepsilon)$ spectra are very narrow. These molecule-like spectra
indicate degraded electrical conductance at room temperature.

Figure 11: Transmission coefficient of an M-GQDA with surface-contacted metal
electrodes and $N_{a}=8$ and $N_{z}=127$ as functions of $\varepsilon$ for
various $\Gamma_{t}$ values. (a), (b), and (c) correspond to $\Gamma_{t}=4.5$,
$9$, and $13.5$ meV, respectively.

Fig. 12 presents the calculated $G_{e}$, $\kappa_{e}$, $PF$, and $ZT$ as
functions of $\mu$ at different temperatures with surface-contacted metal
electrodes. The maximum electrical conductance $G_{e}=0.25G_{0}$ at room
temperature is much smaller than that in Fig.6(a). Unlike the behavior of
$G_{e}$ with respect to temperature, electron thermal conductance $\kappa_{e}$
is enhanced as the temperature increases. Here, $\kappa_{e}$ is in units of
$\kappa_{0}=0.62nW/K$. The increase in the maximum $PF$ with decreasing
temperature is attributable to the enhancement of $S$, which is a highly
nonlinear function of temperature. In the calculation of $ZT$, $\kappa_{ph}$
is the same as that in Fig. 5. At room temperature, the maximum value of
$ZT(\mu=0.535eV)=1.9927$ is less than three. However, the maximum $PF$ and
$ZT$ values increase with decreasing temperature. If the spectra of $G_{e}$
have molecule-like characteristics, the thermoelectric performance of GQDAs
remains acceptable in the low-temperature region [Kuo1, 65].

Figure 12: (a) Electrical conductance, (b) electron thermal conductance, (c)
power factor, and (d) figure of merit as functions of $\mu$ for various
temperatures for M-GQDA with surface-contacted metal electrodes.
$\Gamma_{t}=14.5$ meV. Other physical parameters are the same as in Fig. 11.

The molecule-like or band-like spectra of GQDAs are determined by not only the
lengths of the GQDAs but also the contact geometry. In the line-contacted
metal electrodes, the spectra of finite GQDAs have band-like characteristics
due to the large $\Gamma_{t}$ values. By contrast, the spectra of finite GQDAs
with surface-contacted metal electrodes have molecule-like characteristics. We
did not investigate the effect of electron–electron Coulomb interactions in
this study. Although electron Coulomb interactions are strong in small GQDAs
($L_{a}=0.71$ nm and $L_{z}=15.5$ nm), their effects on the thermoelectric
quantities can be neglected if electron transport is dominated by the
thermionic-assisted tunneling procedure (TATP), which occurs at a distance of
$\mu$ away from the minibands. The thermionic electron population in the
minibands is very dilute if $\mu$ is located inside the gap. Hence, electron
correlation functions arising from electron Coulomb interactions are small
[Kuo1, 65–Kuo3, 67]. The calculated maximum $PF$ and $ZT$ are therefore valid
in the TATP.

##  IV Conclusion

Although ZGNRs have metallic phases, the minibands and gaps are formed in
ZGNRs with periodic vacancies, which behave as GQDAs. Such GQDAs can host
metallic and semiconductor phases. The edge states of the first conductance
and valence minibands, such as $\varepsilon_{e(h),1}$, are localized edge
states along the armchair edge direction, whereas they have very long coherent
lengths along the zigzag edge direction. The designed GQDAs have sufficiently
large gaps and miniband widths to be applied in thermoelectric devices
operating at room temperature. With line-contacted metal electrodes, the
maximum power factor (figure of merit) of the M-GQDAs at room temperature is
larger than that of the S-GQDAs due to the wider miniband widths of the
M-GQDAs. The GQDAs behave as serially coupled QDs (SCQDs) if the armchair edge
atoms are coupled to the metal electrodes. By contrast, the GQDAs behave as
parallel-coupled QDs if zigzag edge atoms are coupled to metal electrodes. For
SCQDs with line-contacted electrodes, the maximum power factor values of GQDAs
with different metal electrodes (Cu, Au, Pt, Pd, or Ti) at room temperature
can reach $79\%$ of the theoretical limit of 1D systems, and their figures of
merit are greater than three. If surface-contacted metal electrodes are used,
the spectra of electrical conductance of the finite M-GQDAs have molecule-like
characteristics because of high contact-surface resistance. The thermoelectric
power factor and figures of merit of the M-GQDAs at room temperature are
smaller than those of the line-contacted M-GQDAs; however, their ZT values can
reach three at low temperatures.

Acknowledgments

We thank Yia-Chung Chang for their valuable intellectual input.

This work was supported by the Ministry of Science and Technology (MOST),
Taiwan under Contract No. 110-2119-M-008-006-MBK.

E-mail address<EMAIL_ADDRESS>

## References

  * (1) Haug H, and Jauho A P, Quantum Kinetics in Transport and Optics of Semiconductors (Springer, Heidelberg, 1996)
  * (2) Guo L J, Leobandung E and Chou S Y, A silicon single-electron transistor memory operating at room temperature, 1997 Science 275 649
  * (3) Postma H W C, Teepen T, Yao Z, Grifoni M and Dekker C, Carbon nanotube single-electron transistors at room temperature, 2001 Science 293 76
  * (4) Michler P, Imamoglu A, Mason M D, Carson P J, Strouse G F and Buratto S K, Quantum correlation among photons from a single quantum dot at room temperature, 2000 Nature 406 968
  * (5) Santori C, Fattal D, Vuckovic J, Solomon G S and Yamamoto Y, Indistinguishable photons from a single-photon device, 2002 Nature 419 549
  * (6) Chang W H, Chen W Y, Chang H S, Hsieh T P, and Hsu T M, Efficient single-photon sources based on low-density quantum dots in photonic-crystal nanocavities, 2006 Phys. Rev. Lett. 96 117401
  * (7) Gustavsson S, Studer M, Leturcq R, Ihn T, Ensslin K, Driscoll D C and Gossard A C, Frequency-selective single-photon detection using a double quantum dot, 2007 Phys. Rev. Lett. 99 206804
  * (8) Josefsson M, Svilans A, Burke A M, Hoffmann E A, Fahlvik S, Thelander C, Leijnse M and Linke H, A quantum-dot heat engine operating close to the thermodynamic efficiency limits, 2018 Nature Nanotechnology 13 920
  * (9) Kagan C R and Murry C B, Charge transport in strongly coupled quantum dot solids, 2015 Nature Nanotechnology 10 1013
  * (10) Harman T C, Taylor P J, Walsh M P and LaForge B E, Quantum dot superlattice thermoelectric materials and devices, 2002 Science 297 2229
  * (11) Talgorn E, Gao Y, Aerts M, Kunneman L T, Schins J M, Savenije T J, Marijn A. van Huis, Herre S. J. van der Zant, Houtepen Arjan J and Siebbeles Laurens D A, Unity quantum yield of photogenerated charges and band-like transport in quantum-dot solids, 2011 Nature Nanotechnology 6733
  * (12) Lawrie W I L, Eenink H G J, Hendrickx N W, Boter J M, Petit L, Amitonov S V, Lodari M, Wuetz Paquelet B, Volk C. and Philips S G J et al, Quantum dot arrays in silicon and germanium, 2020 Appl. Phys. Lett. 116 080501
  * (13) Novoselov K S, Geim A K, Morozov S V, Jiang D, Zhang Y, Dubonos S V, Grigorieva I V and Firsov A A, Electric field effect in atomically thin carbon films, 2004 Science 306 666
  * (14) Geim A K and Grigorieva I V, Van der Waals heterostructures, 2013 Nature 499 419
  * (15) Novoselov K S, Mishchenko A, Carvalho A and Neto A H C, 2D materials and van der Waals heterostructures, 2016 Science 353 aac9439
  * (16) Desai S B, Seol G, Kang J S, Fang H, Battaglia C, Kapadia R, Ager J W, Guo J and Javey A, Strain-Induced Indirect to Direct Bandgap Transition in Multi layer WSe2, 2014 Nano Lett. 14 4592
  * (17) Cai J, Ruffieux P, Jaafar R, Bieri M, Braun T, Blankenburg S, Muoth M, Seitsonen A P, Saleh M, Feng X, Mullen K and Fasel Roman, Atomically precise bottom-up fabrication of graphene nanoribbons, 2010 Nature 466 470
  * (18) Liu J Z, Li B W, Tan Y Z, Giannakopoulos A, Sanchez-Sanchez C, Beljonne D, Ruffieux P, Fasel R, Feng X L and Mullen K, Toward Cove-Edged Low Band Gap Graphene Nanoribbons, 2015 J. Am. Chem. Soc. 137 6097
  * (19) Chen Y C, Cao T, Chen C, Pedramraz Z, Haberer D, Oteyza D G de, F. Fischer R, Louie S G and Crommie M F, Molecular bandgap engineering of bottom-up synthesized graphene nanoribbon heterojunctions, 2015 Nature Nanotechnology 10 156
  * (20) Ruffieux P, Wang S, Yang B, Sachez-Sachez C, Liu J, Dienel T, Talirz L, Shinde P, Pignedoli C A, Passerone D, Dumslaff T, Feng X, Mullen K and Fasel R, On-surface synthesis of graphene nanoribbons with zigzag edge topology, 2016 Nature 531 489
  * (21) Groning O, Wang S, Yao X, Pignedoli C A, Barin G B, Daniels C, Cupo A, Meunier V, Feng X, Narita A, Muellen K, Ruffieux P and Fasel R, Engineering of robust topological quantum phases in graphene nanoribbons, 2018 Nature 560 209
  * (22) Rizzo D J, Veber G, Cao T, Bronner C, Chen T, Zhao F, Rodriguez H, Louie S G, Crommie M F and Fischer F R, Topological band engineering of graphene nanoribbons, 2018 Nature 560 204
  * (23) Yan L H and Liljeroth P, Engineered electronic states in atomically precise artificial lattices and graphene nanoribbons, 2019 ADVANCES IN PHYSICS: X 4 1651672
  * (24) Rizzo D J, Veber G, Jiang J W, McCurdy R, Bronner T, Cao T, Chen T, Louie Steven G, Fischer F R and Crommie M F, Inducing metallicity in graphene nanoribbons via zero-mode superlattices, 2020 Science 369 1597
  * (25) Sun Q, Yan Y, Yao X L, Mullen K, Narita A, Fasel R and Ruffieux P, Evolution of the Topological Energy Band in Graphene Nanoribbons , 2021 J. Phys. Chem. Lett. 12 8679
  * (26) Rizzo D J, Jiang J W, Joshi D, Veber G, Bronner C, Durr R A, Jacobse P H, Cao T, Kalayjian A, Rodriguez H, Butler P, Chen T, Louie Steven G, Fischer F R and Crommie M F, Rationally Designed Topological Quantum Dots in Bottom-Up Graphene Nanoribbons, 2021 ACS Nano 15 20633
  * (27) Wang X, Ma J, Zheng W H, Osella S, Arisnabarreta N, Droste J, Serra J, Ivasenko O, Lucotti A, Beljonne D, Bonn M, Liu X Y, Hansen M R, Tommasini M, Feyter S De, Liu J Z, Wang H I and Feng X L, Cove-Edged Graphene Nanoribbons with Incorporation of Periodic Zigzag-Edge Segments, 2022 J. Am. Chem. Soc. 144 228
  * (28) Almeida P A and Martins G B, Thermoelectric transport properties of armchair graphene nanoribbon heterostructures, 2022 J. Phys: Condens. Matter 34 335302
  * (29) Topsakal M, Sevincli H and Ciraci S, Spin confinement in the superlattices of graphene ribbons, 2008 Appl. Phys. Lett. 92 173118
  * (30) Sevincli H and Cuniberti G, Enhanced thermoelectric figure of merit in edge-disordered zigzag graphene nanoribbons, 2010 Phys. Rev. B 81 113401
  * (31) Darancet P, Olevano V and Mayou D, Coherent Electronic Transport through Graphene Constrictions: Subwavelength Regime and Optical Analogy, 2009 Phys. Rev. Lett. 102 136803
  * (32) Liang G C, Neophytou N, Lundstrom M S and Nikonov D E, Contact effects in graphene nanoribbon transistors, 2008 Nano. Lett. 8 1819
  * (33) Matsuda Y, Deng W Q and Goddard III W A, Contact Resistance for ”End-Contacted” Metal-Graphene and Metal-Nanotube Interfaces from Quantum Mechanics, 2010 J. Phys. Chem. C. 114 17845
  * (34) Shen P C, Su C, Lin Y X, Chou A S, Cheng C C, Park J H, Chiu M H, Lu A Y, Tang H L, Tavakoli M M, Pitner G, Ji X, Cai Z Y, Mao N N, Wang J T, Tung J V, Li C, Bokor J, Zettl A, Wu C I, Palacios T, Li L J and Kong J, Ultralow contact resistance between semimetal and monolayer semiconductors, 2021 Nature 593 212
  * (35) Chen R S, Ding G L, Zhou Y and Han S T, Fermi-level depinning of 2D transition metal dichalcogenide transistors, 2021 J. Mater. Chem. C, 9 11407
  * (36) Wakabayashi K, Sasaki K, Nakanishi T and Enoki T, Electronic states of graphene nanoribbons and analytical solutions, 2010 Sci. Technol. Adv. Mater. 11 054504
  * (37) Kuo D M T, Thermoelectric and electron heat rectification properties of quantum dot superlattice nanowire arrays, 2020 AIP Advances 10 045222
  * (38) Kuo D M T, Effects of zigzag edge states on the thermoelectric properties of finite graphene nanoribbons, 2022 Jpn. J. Appl. Phys. 61 075001
  * (39) Kuo D M T and Chang Y C, Contact Effects on Thermoelectric Properties of Textured Graphene Nanoribbons, 2022 Nanomaterials 12 3357
  * (40) Xu Y, Li Z Y, and Duan W H, Thermal and thermoelectric properties of graphene, 2014 Small 10 2182
  * (41) Nakada K, Fujita M, Dresselhaus G and Dresselhaus M S, Edge state in graphene ribbons: Nanometer size effect and edge shape dependence, 1996 Phys. Rev. B 54 17954
  * (42) Wakabayashi K, Fujita M, Ajiki H and Sigrist M, Electronic and magnetic properties of nanographite ribbons, 1999 Phys. Rev. B 59 8271
  * (43) Gunst T, Markussen T, Jauho A P, and Brandbyge M, Thermoelectric properties of finite graphene antidot lattices, 2011 Phys. Rev. B 84 155449
  * (44) Chang P H, Bahramy M S, Nagaosa N and Nikolic B K, Giant Thermoelectric Effect in Graphene-Based Topological Insulators with Heavy Adatoms and Nanopores, 2014 Nano Lett. 14 3779
  * (45) Zhang Y T, Li Q M, Li Y C, Zhang Y Y and Zhai F, Band structures and transport properties of zigzag graphene nanoribbons with antidot array, 2010 J. Phys: Condens. Matter 22 315304
  * (46) Chu T and Chen Z, Understanding the Electrical Impact of Edge Contacts in Few-Layer Graphene, 2014 ACS Nano 8 3584
  * (47) Gong C, McDonnell S, Qin X Y, Azcatl A, Dong H, Chabal Y J, Cho K, and Wallace R M, Realistic Metal-Graphene contact structures, 2014 ACS Nano 8 642
  * (48) Song S M, Kim T Y, Sul O J, Shin W C and Cho B J, Improvement of graphene-metal contact resistance by introducing edge contacts at graphene under metal, 2014 Appl. Phys. Lett. 104 183506
  * (49) Gao Q and Guo J, Role of chemical termination in edge contact to graphene, 2014 APL Mater. 2 056105
  * (50) Matsuda Y, Deng W Q, and Goddard III W A, Contact resistance properties between nanotubes and various metals from quantum mechanics, 2007 J. Phys. Chem. C 111 11113
  * (51) Massen J, Ji W, and Guo H, First principles study of electronic transport through a Cu(111)/graphene junction. 2010 Appl. Phys. Lett. 97 142105
  * (52) Barraza-Lopez S, Vanevic M, Kindermann M and Chou M Y, Effects of Metallic Contacts on Electron Transport through Graphene, 2010 Phys. Rev. Lett. 104 076807
  * (53) Shen C, Liu J, Jiao N, Zhang C X, Xiao H, Wang R Z and Sun L Z, Transport properties of graphene/metal planar junction, 2014 Phys. Lett. A 378 1321
  * (54) Ran Q, Gao M, Guan X, Wang Y and Yu Z, First principles investigation on bonding formation and electronic structure of metal graphene contacts, 2009 Appl. Phys. Lett. 94 103511
  * (55) Stokbro K, Engelund M and Blom A, Atomic-scale model for the contact resistance of the nickel-graphene interface, 2012 Phys. Rev. B 85 165442
  * (56) Ma B, Cheng G, Wen Y W, Chen R, Cho K G and Shan B, Modulation of contact resistance between metal and graphene by controlling the graphene edge, contact area, and points defects: An ab initio study, 2014 J. Appl. Phys. 115 183708
  * (57) Qin H C, Lu W C and Bernholca J, Ab initio simulations of metal contacts for graphene-based devices, 2022 J. Appl. Phys. 131 214301
  * (58) Yang L, Park C H, Son Y W, Cohen Marvin L, and Louie Steven G, Quasiparticle Energies and Band Gaps in Graphene Nanoribbons, 2007 Phys. Rev. Lett. 99 186801
  * (59) Lee G and Cho K, Electronic structures of zigzag graphene nanoribbons with edge hydrogenation and oxidation, 2009 Phys. Rev. B 79 165440
  * (60) R. S. Whitney, Most Efficient Quantum Thermoelectric at Finite Power Output, 2014 Phys. Rev. Lett.112 130601
  * (61) Zheng H, Liu H J, Tan X J, Lv H Y, Pan L, Shi J and Tang X F, Enhanced thermoelectric performance of graphene nanoribbons, 2012 Appl. Phys. Lett. 100 093104
  * (62) Xu Y, Gan Z and Zhang S C, Enhanced Thermoelectric Performance and Anomalous Seebeck Effects in Topological Insulators, 2014 Phys. Rev. Lett. 112 226801
  * (63) Golizadeh-Mojarad R and Datta S, Effect of contact induced states on minimum conductivity in graphene, 2009 Phys. Rev. B 79 085410
  * (64) Cusati T, Fiori G, Gahoi A, Passi V, Lemme M C, Fortunelli A, and Iannaccone G, Electrical properties of graphene-metal contacts, 2017 Scientific Reports 7 5109
  * (65) Kuo D M T and Chang Y C, Thermoelectric and thermal rectification properties of quantum dot junctions, 2010 Phys. Rev. B 81 205321
  * (66) Kuo D M T, Shiau S Y and Chang Y C , Theory of spin blockade, charge ratchet effect, and thermoelectrical behavior in serially coupled quantum dot system, 2011 Phys. Rev. B 84 245303
  * (67) Kuo D M T, Chen C C and Chang Y C, Large enhancement in thermoelectric efficiency of quantum dot junctions due to increase of level degeneracy, 2017 Phys. Rev. B 95 075432


# A Nonlinear Formulation of Radiation Stress and Applications to Cnoidal
Shoaling

Martin O. Paulsen $\cdot$ Henrik Kalisch

###### Abstract

In this article we provide formulations of energy flux and radiation stress
consistent with the scaling regime of the Korteweg-de Vries (KdV) equation.
These quantities can be used to describe the shoaling of cnoidal waves
approaching a gently sloping beach. The transformation of these waves along
the slope can be described using the shoaling equations, a set of three
nonlinear equations in three unknowns: the wave height $H$, the set-down
$\bar{\eta}$ and the elliptic parameter $m$. We define a numerical algorithm
for the efficient solution of the shoaling equations, and we verify our
shoaling formulation by comparing with experimental data from two sets of
experiments as well as shoaling curves obtained in previous works.

g

### Keywords:

Surface waves $\cdot$ KdV equation $\cdot$ Conservation laws $\cdot$ Radiation
stress $\cdot$ Set-down

##  1 Introduction

In the present work, the development of surface waves across a gently sloping
bottom is in view. Our main goal is the prediction of the wave height and the
set-down of a periodic wave as it enters an area of shallower depth. This
problem is known as the shoaling problem and has a long history, with
contributions from a large number of authors going back to Green and
Boussinesq.

In order to reduce the problem to the essential factors, we make a number of
simplifying assumptions. First, the bottom slope is small enough to allow the
waves to continuously adjust to the changing depth without altering their
basic shape or breaking up. Second, reflections are assumed to be negligible
so that the energy flux of a wave is conserved as it transforms on the slope.
Third, the period $T$ of a wave is assumed to be constant as it progresses.
This last point is sometimes termed the conservation of waves as the number of
waves in a group will remain constant during the shoaling process.

Rather than following a wave in space and time as it propagates over the
slope, we estimate wave properties as functions of the depth using
conservation of period $T$ and energy flux, and a prescribed change in
momentum due to the effect of the radiation stress. These conservation
equations are due to the assumptions stated above, and they form the basis for
the description of the shoaling process as first envisioned by Rayleigh [34],
and subsequently used by a number of authors. The crucial factor which
guarantees that the shoaling assumptions are valid is that the relative change
in water depth over a wavelength is smaller than the wave steepness. This
point is explained in more detail in [32, 43].

In the context of linear wave theory, the process outlined above is the
classical approach to wave shoaling and can be found in textbooks on coastal
engineering such as [12, 40]. This approach depends on the assumption that the
waves can be described by a sinusoidal function of the form
$\frac{H}{2}\cos(kx-\omega t)$, where $H$ is the wave height, $k$ is the
wavenumber, $\omega$ is the circular frequency, $t$ is the time and $x$ is a
spatial coordinate along which the bottom is sloping. Recall that the
wavelength is $\lambda=\frac{2\pi}{k}$ and the wave period is
$T=\frac{2\pi}{\omega}$, and these are related by the dispersion relation
$\omega^{2}=gk\tanh{(kh)}$ at local depth $h$. In the linear shoaling theory,
the conservation of energy flux is central. The energy flux is generally
formulated as the average over one period of the energy density $E$ times the
group speed $C_{g}$. The energy density is given in terms of the fluid density
$\rho$, the gravitational acceleration $g$ and the wave height $H$ as
$E=\frac{1}{8}\rho gH^{2}$, and the group speed is given by
$C_{g}=\frac{d\omega}{dk}$.

If the bathymetry is sloping only gently upwards, then it may be assumed that
the waves adjust adiabatically to the changing conditions, and that
reflections and distortions are negligible. The wave-height transformation of
a shoaling wave is obtained by imposing conservation of the energy flux
$EC_{g}$ across the shoaling region. If a wave measurement $H_{0}$, $T_{0}$ at
some point offshore is given, then the wave height at a point closer to shore
may be computed from the conservation of energy flux given in the form

$HC_{g}=H_{0}C_{g0},$ (1)

where the group speed $C_{g}$ at local depth $h$ can be found from the
conservation of the wave period $T$. Similar considerations using the linear
definition of the radiation stress yield a formula for the set-down [28].

The authors of [44] used cnoidal functions in connection with the linear
formulation of the energy flux defined above to obtain a hybrid theory of
shoaling. Unfortunately, this approach led to a discontinuity in wave height
at the matching point between the linear shoaling equation (1) and the
nonlinear theory based on the cnoidal functions. The problem was remedied to
some degree in [45] by requiring continuity in wave height effected through
the use of conversion tables. This approach led to good agreement with the
experimental data, but was cumbersome to implement, and as already noted in
[45], the continuity of the wave height at the matching point in the shoaling
curve led to a discontinuity in energy flux.

The problem can be resolved by defining a shoaling equation based on the full
water-wave problem using the streamfunction method [12] for example. Such an
approach has been documented in [46]. While accurate, the streamfunction
approach depends on discretizing the steady Euler equations leading to large
systems of equations which need to be solved numerically. If the approach is
to be used in connection with ocean-wave statistics, the approach based on the
steady Euler equations might be too computationally demanding. Moreover, the
agreement of the combined linear-cnoidal approach of [45] already gave good
agreement with experiments while being much less expensive.

In present work we will show that the discrepancy between conservation of wave
height and conservation of energy flux can be resolved in the context of the
KdV theory by using a definition of the energy flux consistent with the KdV
equation [4, 15], and without using the fully nonlinear approach of [46].
Recall that the KdV equation is given in the form

$\eta_{t}+c_{0}\eta_{x}+\frac{3}{2}\frac{c_{0}}{h_{0}}\eta\eta_{x}+\frac{c_{0}h_{0}^{2}}{6}\eta_{xxx}=0,$
(2)

where $\eta(x,t)$ is the deflection of the free surface from rest at a point
$x$ and a time $t$, $g$ is the gravitational acceleration, $h_{0}$ is the
local water depth, and $c_{0}=\sqrt{gh_{0}}$ is the long-wave speed.

Shoaling in the context of steady solutions of the KdV equation is then
described as follows. Given a steady wavetrain of wave height $H$, wavelength
$\lambda$ and period $T$, the variation of the wave as it propagates from
depth $h^{A}$ to $h^{B}$ is obtained by imposing conservation of $T$,
conservation of the energy flux $q_{E}$, and balance of forces using the
bottom forcing by the pressure force, and the radiation stress.

In order to execute this plan, one needs to have in hand formulations for the
energy flux $q_{E}$ and the radiation stress $S_{xx}$ in the context of the
KdV equation. These expressions can be developed using ideas first put forward
in [3, 4]. Using the methods in these papers, it can be seen that we have

$q_{E}=c_{0}^{3}\Big{(}\frac{1}{h_{0}}\eta^{2}+\frac{5}{4h_{0}^{2}}+\frac{h_{0}}{2}\eta\eta_{xx}\Big{)}$

and

$\overline{S}_{xx}=\rho
g\Big{(}\overline{\eta}+\frac{3}{2}\overline{\eta^{2}}+\frac{h_{0}^{3}}{3}\overline{\eta}_{xx}\Big{)}.$

Note that the energy flux is defined in a pointwise sense while the radiation
stress is defined in an average sense. However, for the shoaling problem
$q_{E}$ will also be averaged over one wave period.

The plan of the paper is as follows. In the next section, the formulation of
momentum and energy balance in the context of the KdV equation will be
recalled, leading to the expression for flow force $q_{I}$ and energy flux
$q_{E}$. In Section 3, a formulation of the radiation stress consistent with
the KdV equation will be found. Section 4 contains the formulation of the
nonlinear shoaling equations, details on the numerical implementation and a
comparison with wave tank experiments with particular focus on the set-down.
Section 6 details the comparison with other nonlinear shoaling theories
without set-down, and in particular with the shoaling curves obtained by
Svendsen and Buhr Hansen [45]. In Section 7, we provide a comparison with the
experimental results of [45]. In the Conclusion we put our work into context
and mention possibilities for further work.

##  2 Momentum and energy balance in the KdV approximation

We start this section with a brief description of the water-wave problem of an
inviscid, incompressible and homogeneous fluid with a free surface. Due to the
assumption that the relative change in water depth over one wavelength is less
than the wave steepness, the problem can be locally formulated with a constant
undisturbed depth $h_{0}$ (see Figure 1).

Figure 1: A periodic wave of wave height $H$ and wavelength $\lambda$
propagating at the free surface of a fluid of undisturbed depth $h_{0}$. The
free surface is given by $z=\eta(x,t)$.

If the flow is assumed to be irrotational, the problem can be formulated in
terms of a velocity potential $\phi(x,z,t)$ in addition to the surface
deflection $\eta(x,t)$. The velocity field $(u(x,z,t),v(x,z,t))$ is then given
as the gradient of $\phi$, and the pressure $P$ can be found using the
Bernoulli equation. In terms of $\phi$ and $\eta$, the problem is described by
the Laplace equation

$\hskip 99.58464pt\phi_{xx}+\phi_{zz}=0\hskip 14.22636pt\text{for}\hskip
5.69046pt-h_{0}<z<\eta(x,t)$ (3)

in the fluid domain, and the boundary conditions

$\displaystyle\eta_{t}+\phi_{x}\eta_{x}-\phi_{z}$ $\displaystyle=0\hskip
14.22636pt\text{on}\hskip 5.69046ptz=\eta(x,t),$
$\displaystyle\phi_{t}+\frac{1}{2}(\phi_{x}^{2}+\phi_{z}^{2})+g\eta$
$\displaystyle=0\hskip 14.22636pt\text{on}\hskip 5.69046ptz=\eta(x,t),$
$\displaystyle\phi_{z}$ $\displaystyle=0\hskip 14.22636pt\text{on}\hskip
5.69046ptz=-h_{0}.$

It is well known that this problem is difficult to treat both mathematically
and numerically. In particular, it is not known whether solutions exist on
relevant time scales [22], and numerical discretization of the full problem on
the field scale remains completely out of reach. Thus in practical situations,
an asymptotic approximation of the Euler equations is usually required. The
traditional asymptotic theories are the linearized (Airy) theory and the
hyperbolic (Saint-Venant) shallow-water theory such as outlined in [42]. The
work of Boussinesq [8] and Korteweg and de Vries [18] led to another
asymptotic regime, the so-called Boussinesq scaling which can be used for long
waves of small to moderate amplitude. In the present work, we consider waves
in the Boussinesq regime, and since we consider waves propagating in the
shoreward direction only, we restrict attention to the Korteweg-de Vries (KdV)
equation. Indeed, the KdV equation describes unidirectional waves in the
presence of weakly nonlinear and dispersive effects in the case when these
corrections are approximately balanced. In order to make this statement
precise, it is useful to define the relative amplitude $\alpha=a/h_{0}$ and
the shallowness parameter $\beta=h_{0}^{2}/\lambda^{2}$, where $a=H/2$ denotes
a representative amplitude and $\lambda$ a representative wavelength of the
wavefield to be studied. Using an asymptotic expansion and a dimensional
argument it can be shown that the KdV equation is a good model for long waves
of small amplitude if terms of order
$\mathcal{O}(\alpha^{2},\alpha\beta,\beta^{2})$ are neglected.

In order to derive the KdV equation, one defines the non-dimensional variables
$x=\lambda\tilde{x}$, $z=h_{0}(\tilde{z}-1)$,
$t=\frac{\lambda}{c_{0}}\tilde{t}$, and
$\phi=\frac{ga\lambda}{c_{0}}\tilde{\phi}$, then assumes that the velocity
potential takes the form

$\displaystyle\tilde{\phi}=\tilde{f}-\beta\frac{\tilde{z}^{2}}{2}\tilde{f}_{\tilde{x}\tilde{x}}+\mathcal{O}(\beta^{2}),$
(4)

where $\tilde{f}$ is the velocity potential evaluated at the bed. Following
the description in [48], the free-surface boundary conditions can then be used
to obtain the KdV equation (2). In the process, it comes to light that the
horizontal velocity at the bottom for a right-moving wave is given by

$\tilde{u}=\tilde{\eta}-\frac{1}{4}\alpha\tilde{\eta}^{2}+\frac{1}{3}\beta\tilde{\eta}_{\tilde{x}\tilde{x}\tilde{x}}+\mathcal{O}(\alpha^{2}+\beta^{2}).$
(5)

Moreover, as noted for example in [4], combining $\eqref{horizontal velocity}$
and the derivative of $\eqref{non dim potential}$ with respect to $x$ one
finds that the horizontal component of the velocity field is given by

$\tilde{\phi}_{\tilde{x}}(\tilde{x},\tilde{z},\tilde{t})=\tilde{\eta}-\frac{1}{4}\alpha\tilde{\eta}^{2}+\beta\Big{(}\frac{1}{3}-\frac{\tilde{z}^{2}}{2}\Big{)}\tilde{\eta}_{\tilde{x}\tilde{x}}+\mathcal{O}(\alpha^{2},\alpha\beta,\beta^{2})$
(6)

in the KdV approximation. This relation will be needed later in the derivation
of the the momentum and energy balance, and the formulation of the radiation
stress. In addition, we need the vertical component of the velocity field and
the pressure. The vertical velocity is given by

$\displaystyle\tilde{\phi}_{\tilde{z}}(\tilde{x},\tilde{z},\tilde{t})=-\beta\tilde{z}\tilde{\eta}_{\tilde{x}}+\mathcal{O}(\alpha\beta,\beta^{2}).$

In order to express the pressure in terms of the surface excursion $\eta$,
assuming unit density for the moment, one first defines the dynamic pressure
in the usual way:

$P^{\prime}=P-P_{\text{atm}}+gz,\quad P^{\prime}=ag\tilde{P}^{\prime}.$

Then using the Bernoulli equation and following the computations outlined in
[4] leads to the non-dimensional perturbation pressure in the KdV
approximation:

$\tilde{P}^{\prime}=\tilde{\eta}-\frac{1}{2}\beta(\tilde{z}^{2}-1)\tilde{\eta}_{\tilde{x}\tilde{x}}+\mathcal{O}(\alpha\beta,\beta^{2}).$
(7)

Approximations of the momentum and energy balance laws in the KdV scaling
regime will now be derived. In the context of the full Euler equations with
surface boundary conditions, the momentum balance is expressed by

$\frac{\partial}{\partial
t}\int_{-h_{0}}^{\eta}\phi_{x}dz+\frac{\partial}{\partial
x}\int_{-h_{0}}^{\eta}\big{\\{}\phi_{x}^{2}+P\big{\\}}dz=0.$ (8)

Using the previously defined scaling, the momentum balance is

$\alpha\frac{\partial}{\partial\tilde{t}}\int_{0}^{1+\alpha\tilde{\eta}}\tilde{\phi}_{\tilde{x}}h_{0}d\tilde{z}+\frac{\partial}{\partial\tilde{x}}\int_{0}^{1+\alpha\tilde{\eta}}\big{\\{}\alpha^{2}\tilde{\phi}^{2}_{\tilde{x}}+\alpha\tilde{P}^{\prime}+(1-\tilde{z})\big{\\}}h_{0}d\tilde{z}=0.$

Substituting the two expressions $\tilde{\phi}_{x}$, $\tilde{P}^{\prime}$ and
evaluating the integral one finds

$\Big{(}\alpha\tilde{\eta}+\frac{3}{4}\alpha^{2}\tilde{\eta}^{2}+\frac{1}{6}\alpha\beta\tilde{\eta}_{\tilde{x}\tilde{x}}\Big{)}_{\tilde{t}}+\Big{(}\frac{1}{2}+\alpha\tilde{\eta}+\frac{3}{2}\alpha\tilde{\eta}^{2}+\frac{1}{3}\beta\tilde{\eta}_{\tilde{x}\tilde{x}}\Big{)}_{\tilde{x}}=\mathcal{O}(\alpha^{2},\alpha\beta,\beta^{2}).$

From this relation we identify the non-dimensional momentum density

$\tilde{I}=\alpha\tilde{\eta}+\frac{3}{4}\alpha^{2}\tilde{\eta}^{2}+\frac{1}{6}\alpha\beta\tilde{\eta}_{\tilde{x}\tilde{x}},$

and the non-dimensional momentum flux

$\tilde{q_{I}}=\frac{1}{2}+\alpha\tilde{\eta}+\frac{3}{2}\alpha\tilde{\eta}^{2}+\frac{1}{3}\beta\tilde{\eta}_{\tilde{x}\tilde{x}}.$

Returning the expression to its dimensional forms through the scaling
$I=c_{0}h_{0}\tilde{I}$ and $q_{I}=c_{0}^{2}h_{0}\tilde{q}_{I}$ yields

$I=c_{0}\Big{(}\eta+\frac{3}{4h_{0}}\eta^{2}+\frac{h_{0}^{2}}{6}\eta_{xx}\Big{)},$

and

$q_{I}=c_{0}^{2}\Big{(}\frac{h_{0}}{2}+\eta+\frac{3}{2h_{0}}\eta^{2}+\frac{h_{0}^{2}}{3}\eta_{xx}\Big{)}.$
(9)

Note that the expression for $q_{I}$ combines the momentum flux and the
pressure force. Following Benjamin and Lighthill [6], we will use the term
flow force for this quantity. Similarly, we can give the energy balance in the
full Euler equations as

$\frac{\partial}{\partial
t}\int_{-h_{0}}^{\eta}\Big{\\{}\frac{1}{2}|\nabla\phi|^{2}+g(z+h_{0})\Big{\\}}dz+\frac{\partial}{\partial
x}\int_{-h_{0}}^{\eta}\Big{\\{}\frac{1}{2}|\nabla\phi|^{2}+g(z+h_{0})+P\Big{\\}}\phi_{x}dz=0,$

and following the same procedure as above will lead to the expressions

$E=c_{0}^{2}\Big{(}\frac{1}{h_{0}}\eta^{2}+\frac{1}{4h_{0}^{2}}\eta^{3}+\frac{h_{0}}{6}\eta\eta_{xx}+\frac{h_{0}}{6}\eta_{x}^{2}\Big{)},$

and

$q_{E}=c_{0}^{3}\Big{(}\frac{1}{h_{0}}\eta^{2}+\frac{5}{4h_{0}^{2}}\eta^{3}+\frac{h_{0}}{2}\eta\eta_{xx}\Big{)},$
(10)

for the energy density and energy flux respectively.

Note that $q_{E}$ is of second order, while the energy flux in the linear
approximation is of first order. Indeed, it will be instructive to compare
(10) to the well known expression for the energy flux in the linear theory.
Denoting the amplitude of a linear wave by $a=H/2$ and the group velocity by
$C_{g}$ as before, it can shown (see for example [21]) that the energy flux
averaged over one period is

$\overline{q_{E}^{\mathrm{linear}}}=g\overline{\eta^{2}}C_{g}=\frac{1}{2}ga^{2}C_{g}.$

On the other hand, for cnoidal waves of very large wavelength and very small
amplitude (10) can be approximated by

$q_{E}\sim c_{0}^{3}\frac{1}{h_{0}}\eta^{2}=c_{0}g\eta^{2}.$

Averaging over one period and using the fact that the cnoidal wave is similar
to a sinusoidal function for very small amplitudes, we obtain

$\overline{q_{E}}\sim\frac{1}{2}ga^{2}c_{0}.$

Finally, the two expressions are seen to be approximately equal if it is
recognized that $c_{0}$ is the group velocity for long linear waves in shallow
water.

Note that in previous works on cnoidal shoaling such as [44, 45], the linear
formulation of the energy flux was used. In addition to being obviously
inconsistent, this approach also had practical drawbacks. Indeed as already
mentioned, the method put forward in [45, 44] led to a discontinuity in the
wave height at the matching point between linear and nonlinear theory. As will
be shown momentarily, keeping the nonlinear terms in the expression for
$q_{E}$ resolves this problem, and indeed is essential in order to obtain the
correct shoaling behavior for larger wave amplitudes.

In the next section, we will use the definition of the flow force $q_{I}$, and
ideas outlined above to derive an expression for the radiation stress in the
context of the KdV equation. Incorporating the radiation stress into the
shoaling equations will then enable us to identify the set-down in addition to
the wave height of a shoaling wave.

##  3 Radiation stress in the KdV approximation

In linear water-wave theory, the radiation stress is well understood as a
second-order response to a periodic wave train which is akin to the energy
flux. Wave radiation stress is an important ingredient in the study of wave-
current interactions, and also plays a major role in the understanding of wave
setup in the context of breaking waves. The reader may consult [27] for the
mathematical development of an expression for the radiation stress in the
context of linear waves. A discussion of the radiation stress based on a more
heuristic physical understanding is given by Longuet-Higgins and Stewart [28].
In the present work, we aim to expand the definition of the radiation stress
to the nonlinear case, in particular in the context of the KdV equation. Using
the momentum balance equation derived in the previous section, a nonlinear
version of the radiation stress will be formulated. With this expression in
hand, some consequences of the radiation stress on wave shoaling will be
investigated.

Analyzing the definition given in [28], one can simply think of radiation
stress as the total flow force of a progressive wave averaged over one period
minus the hydrostatic pressure force at rest. As previously discussed we
consider a wave propagating solely in the $x-$direction, and neglect all
transverse effects. Then the definition of the principal component of the
radiation stress is given by

$\overline{S}_{xx}=\overline{\int_{-h_{0}}^{\eta}(\rho
u^{2}+P)dz}+\int_{-h_{0}}^{0}\rho gzdz.$

The first term expresses the total flux of momentum across a plane integrated
from the bottom to the free surface and with unit width, and as usual, the
overbar denotes that the average over one wave period with respect to time is
taken. The second term expresses the flow force in the absence of any wave
motion.

In the context of the KdV equation, the flow force is given by equation (9)
when rescaled for a fluid with unit density. Consequently, the $x-$component
of the radiation stress is obtained in the KdV theory as

$\displaystyle\overline{S}_{xx}$
$\displaystyle=\overline{q_{I}}+\int_{-h_{0}}^{0}\rho gzdz,$
$\displaystyle=\rho
gh_{0}\Big{(}\frac{h_{0}}{2}+h_{0}\overline{\eta}+\frac{3}{2h_{0}}\overline{\eta^{2}}+\frac{h_{0}^{2}}{3}\overline{\eta}_{xx}\Big{)}-\frac{1}{2}\rho
gh_{0}^{2},$ $\displaystyle=\rho
g\Big{(}h_{0}\overline{\eta}+\frac{3}{2}\overline{\eta^{2}}+\frac{h_{0}^{3}}{3}\overline{\eta}_{xx}\Big{)}.$
(11)

This formulation represents an improvement over the formulation given in [43]
which only recorded the middle term $\rho g\frac{3}{2}\overline{\eta^{2}}$ in
the final expression for the radiation stress. A preliminary version of this
quantity was also presented in [19]. In order to apply the radiation stress to
the formulation of the shoaling problem, let us consider a wave encountering a
gently sloping beach. As indicated in Figure 2, the momentum flux is reduced
in the onshore direction due to an opposing horizontal force exerted by the
bed, and opposing the fluid pressure.

Figure 2: Schematic of the forces acting on a control volume given by a
differential interval of length $dx$ and reaching through the whole fluid
column.

When accounting for the set-down of the mean surface level, the flow force on
the offshore side of the boundary of the control volume is given in terms of
radiation stress by

$q_{I}=S_{xx}+\int_{-h}^{\overline{\eta}}\rho
g(\overline{\eta}-z)dz=S_{xx}+\frac{1}{2}\rho g(\overline{\eta}+h)^{2}.$

This relation can be found for example in [12, 28]. As already mentioned
above, the momentum flux is reduced by the reaction force exerted by the sea
bed on the fluid due to the weight of the fluid. This force is equal in
magnitude to the horizontal component of the pressure force at the bottom
$\rho g(h+\overline{\eta})dz$, where $dz$ denotes the vertical variation in
the sea bed over the horizontal distance $dx$. Thus the horizontal component
of the pressure force exerted on the bottom can be written as

$\rho g(h+\overline{\eta})dz=\rho g(h+\overline{\eta})h_{x}dx.$ (12)

As indicated in Figure 2, the difference in flow force between the shoreward
face and the offshore face of the control volume is given by

$\frac{d}{dx}\Big{[}S_{xx}+\frac{1}{2}\rho g(h+\overline{\eta})^{2}\Big{]}dx.$
(13)

Evidently, momentum balance requires (12) and (13) to be equal, and
simplifying yields the relation

$\frac{dS_{xx}}{dx}+\rho g(\overline{\eta}+h)\frac{d\overline{\eta}}{dx}=0.$
(14)

In order to develop the shoaling equation, we integrate over the control
interval $[x,x+dx]$ to get

$\int_{x}^{x+dx}\frac{dS_{xx}}{dx}+\int_{x}^{x+dx}\rho
g\overline{\eta}\frac{d\overline{\eta}}{dx}+\int_{x}^{x+dx}\rho
gh\frac{d\overline{\eta}}{dx}=0.$

The first two integrals are straightforward to compute. The third integral is
evaluated by using integration by parts and then approximated using the
trapezoidal rule. Defining the difference of any quantity $F$ as $\Delta
F=F|_{x+dx}-F|_{x}$ we get, the relation

$\Delta S_{xx}=-\frac{\rho g}{2}\Delta\overline{\eta}^{2}-\frac{\rho
g}{2}(2h+h_{x}dx)\Delta\overline{\eta}.$ (15)

Recall that we are interested in the averaged behavior of the wave height for
a wave-train approaching the beach. We therefore find it useful to define the
changes in radiation stress averaged over a period by

$\Delta\overline{S}_{xx}=-\frac{1}{T}\int_{0}^{T}\Big{\\{}\frac{\rho
g}{2}\Delta\overline{\eta}^{2}+\frac{\rho
g}{2}(h_{0}+h)\Delta\overline{\eta}\Big{\\}}dt,$ (16)

where $\overline{S}_{xx}$ is defined by (11). With this identity in hand, it
will be possible to include the development of the set-down $\overline{\eta}$
in the shoaling problem. The formulation of the complete shoaling equation
will be the subject of the next section.

##  4 The nonlinear shoaling equations

An explicit set of equations describing the shoaling of long waves on a gently
sloping beach will now be given. The idea is to use the well known cnoidal
wave solution for periodic waves in the KdV equation together with the
assumption that the bottom slope is gentle enough so that reflections can be
neglected, and the waves are able to adjust adiabatically to the new local
depth. The wave profile $\eta$ will distort, but the underlying cnoidal shape
and periodicity are preserved.

A solution of the KdV equation for periodic waves of constant form was first
discovered by Korteweg and de Vries in 1895 [18]. Assuming a traveling-wave
solution of constant shape, one may make the ansatz $\eta(x,t)=f(x-ct)$, and
reduce the KdV equation to an ordinary differential equation of the form

$-\frac{h_{0}^{2}}{3}(f^{\prime})^{2}=F(f)=\frac{1}{4}\frac{c_{0}}{h_{0}}f^{3}+\Big{(}\frac{c-c_{0}}{2}\Big{)}f^{2}+Af+B,$
(17)

with $A$ and $B$ being constants of integration. This differential equation
has a number of solutions, but for practical applications only periodic real-
valued and bounded solutions correspond to realistic wave profiles. Since
$F(f)$ is a third-order polynomial it has three roots. Periodic solutions
occur only in the case where the roots are distinct, so that they can be
labeled $f_{3}<f_{2}<f_{1}$, and this convention allows us to write

$F(f)=(f-f_{1})(f-f_{2})(f-f_{3}).$

Figure 3: Graph of the function $-F(f)$. The roots are labelled by
$f_{3}<f_{2}<f_{1}$. The red lines denote the values of $f$ such that $-F\geq
0$.

Requiring the solution to be real, it can be seen from (17) that $-F(f)$ must
be positive. Examining the phase plane plot in Figure 3, it is clear that
periodic solutions will oscillate between $f_{2}$ and $f_{1}$. Since
$f_{2}<f_{1}$ by assumption, $f_{1}$ will denote the wavecrest while $f_{2}$
will be the trough. Consequently, the wave height $H$ is given by the
difference $H=f_{1}-f_{2}$. As is well known, the periodic solutions of (17)
are given in terms of the Jacobian elliptic function $cn$. The solution of the
KdV equation is then

$\eta(x,t)=f_{2}-(f_{2}-f_{1})\text{cn}^{2}\Big{(}2K(m)\big{(}{\textstyle\frac{t}{T}}-{\textstyle\frac{x}{\lambda}}\big{)},m\Big{)},$
(18)

where $m$ is the elliptic modulus defined by
$m=\frac{f_{1}-f_{2}}{f_{1}-f_{3}}$, and $K(m)$ is the complete elliptic
integral of first kind [24]. The wavespeed $c$ and the wavelength $\lambda$
are then given by

$\hskip
28.45274ptc=c_{0}\Big{(}1+\frac{f_{1}+f_{2}+f_{3}}{2h_{0}}\Big{)},\hskip
28.45274pt\lambda=K(m)\sqrt{\frac{16h_{0}^{3}}{3(f_{1}-f_{3})}}.$ (19)

In order to use these formulas in practice, it is convenient to take the wave
height $H$, the mean surface level $\overline{\eta}$ and the elliptic
parameter $m$ as given parameters. The roots of $F$ can then by computed
explicitly and are given as follows.

$\begin{dcases}f_{3}&=\overline{\eta}-\frac{HE(m)}{mK(m)},\\\
f_{1}&=f_{3}+\frac{H}{m},\\\ f_{2}&=f_{1}-H.\end{dcases}$ (20)

The shoaling problem can be formulated as follows. Suppose we know that a
periodic wave train of wave height $H$, wavelength $\lambda$ and mean-surface
level $\overline{\eta}$ is given at a local depth $h^{A}$. Assuming that the
waves can be approximately described by the the KdV equation, a unique cnoidal
wave solution can be found, and the frequency $\nu$, the energy flux $q_{E}$
and the radiation stress $S_{xx}$ be found. The shoaling problem now consists
of finding an appropriate cnoidal solution at a smaller depth $h^{B}$. Relying
on the conservation of frequency and energy flux, and a prescribed change in
the radiation stress, the shoaling equations can be written as

$\displaystyle\nu^{B}(f_{1},f_{2},f_{3},h^{B})\ =$ $\displaystyle\ \ \nu^{A}\
=\ \frac{c^{A}}{\lambda^{A}},$ (21)
$\displaystyle\overline{q}_{E}^{B}(f_{1},f_{2},f_{3},h^{B})\ =$
$\displaystyle\ \ \overline{q}_{E}^{A}\ =\
\frac{1}{T}\int_{0}^{T}\Big{\\{}\frac{1}{h^{A}}\eta^{2}+\frac{5}{4(h^{A})^{2}}\eta^{3}+\frac{h^{A}}{2}\eta\eta_{xx}\Big{\\}}\hskip
2.84526ptdt,$ (22)
$\displaystyle\overline{S}^{B}_{xx}(f_{1},f_{2},f_{3},h^{B})\ =$
$\displaystyle\ \
\overline{S}^{A}_{xx}-\frac{1}{T}\int_{0}^{T}\Big{\\{}\frac{\rho
g}{2}\Delta\eta^{2}+\frac{\rho g}{2}(h^{B}+h^{A})\Delta\eta\Big{\\}}dt.$ (23)

The first equation is conservation of frequency $\nu$, while the second
equation expresses conservation of energy flux integrated over a period $T$
[4]. Finally, the third equation is (16), and is indeed a formulation of
conservation of momentum expressed in terms of radiation stress. In previous
work [20], a similar system was found, and then solved for the parameters
$f_{1}$, $f_{2}$ and $f_{3}$. However, the method used in [20] had several
drawbacks. First, due to a lack of a definition of radiation stress, the
equation for momentum conservation was replaced by conservation of the mean
fluid level. This approach was therefore not able to predict the wave set-
down. Moreover, solving for the parameters $f_{1}$, $f_{2}$ and $f_{3}$ was
not optimal from a numerical point of view, and the algorithm terminated long
before the highest wave was reached. In the present work, we use the explicit
form of the parameters to solve directly for $H$, $m$ and $\overline{\eta}$.
This approach is more expedient from a numerical point of view, and allows us
to go much higher up on the shoaling curve. In terms of the wave height $H$,
the conservation of frequency turns into the following cubic equation for $H$:

$\frac{-3g(\frac{3E(m)}{mK(m)}-\frac{2}{m}+1)^{2}}{64K(m)^{2}h^{4}m}H^{3}+\frac{3g(\frac{3\overline{\eta}}{2h}+1)(\frac{3E(m)}{mK(m)}-\frac{2}{m}+1)}{16K(m)^{2}h^{3}m}H^{2}+-\frac{3g(\frac{3\overline{\eta}}{2h}+1)^{2}}{16K(m)^{2}h^{2}m}H+\frac{1}{T^{2}}=0.$
(24)

Similarly we can manipulate the equation describing the changes in radiation
stress to find an expression for the set-down. Indeed, (23) reduces to the
quadratic equation

$\displaystyle-\overline{\eta}^{2}+\mathcal{A}\overline{\eta}+\mathcal{B}=0,$
(25)

with

$\displaystyle\mathcal{A}=\Big{(}3H-\frac{3h^{B}}{2}+\frac{h^{A}}{2}-3H\overline{\text{cn}^{2}}(m)-\frac{3H}{m}+\frac{3HE(m)}{mK(m)}\Big{)}$

and

$\displaystyle\mathcal{B}=\overline{S}^{A}_{xx}(m)-\frac{\overline{\eta^{2}}^{A}(m)}{2}-\frac{3}{2}H^{2}\overline{\text{cn}^{4}}(m)-\frac{(h^{B})^{3}\overline{\eta_{xx}}^{B}(m)}{3}+\frac{3}{2}h^{A}\overline{\eta}^{A}(m)-\frac{1}{2}h^{B}\overline{\eta}^{A}(m)$
$\displaystyle-\frac{3}{2}(H-\frac{H}{m}+\frac{HE(m)}{mK(m)})^{2}+3H\overline{\text{cn}^{2}}(m)(H-\frac{H}{m}+\frac{HE(m}{mK(m)}).$

Note that we had to integrate various powers and derivatives of
$\text{cn}^{2}({\xi;m})$. These formulas are based on calculations that can be
found in [24] and [1], and are given in explicit form in the appendix. Having
this representation in hand, we can in principle write (24) as
$H=F(m,\overline{\eta})$ and take (25) to be $\overline{\eta}=G(m,H)$ in
explicit terms as two coupled equations. To solve these two equations
simultaneously, we may iterate between them at current local depth $h^{B}$ as
follows; first initialize the procedure with $\overline{\eta}^{0}(m)$ and then
find $H$ by

$H(m)^{i+1}=F(m,\overline{\eta}^{i}(m)),$

solving the cubic equation. This can in turn be used to update the set-down by

$\overline{\eta}(m)^{i+1}=G(m,H^{i+1}(m)).$

Repeating this process, we continue to approximate $H$ and $\overline{\eta}$
until a stopping criterion has been reached. Using this reduction allows us to
solve the nonlinear system (21),(22),(23), by the following procedure. From
the value of $m$ at depth $h^{A}$, we increment up, at each step iterating the
above equations for $H$ and $\overline{\eta}$, and then checking whether (22)
is satisfied to a specified tolerance. If that is the case, we consider the
system solved at depth $h^{B}$, and move on to the next step with depth
$h^{C}$, thereby moving up the slope.

##  5 Implementation of the shoaling equations

Having explained the nonlinear shoaling equations, and the strategy for
finding approximate solutions, we will now implement the equations and produce
shoaling curves for various deep-water data. The method used here proceeds in
three stages, such as originally developed in [44, 45]. Suppose an incoming
wave of wavelength $\lambda_{0}$ is given at a starting depth $h_{0}$. While
this depth may not be deep water, it is assumed that $h_{0}>0.28\lambda_{0}$,
so that we label this starting value the deep-water wavelength. First, the
linear shoaling equation is used up to the point $h/\lambda_{0}=0.1$. At this
point the cnoidal theory is valid [39], and we propose a matching technique to
obtain the fundamental parameters of the nonlinear wave. Lastly, we use the
nonlinear KdV theory to follow the shoaling curve into the nonlinear region.

Stage 1. Linear theory: The essential ingredient in this step is the linear
dispersion relation

$\omega^{2}-gk\tanh{(kh)}=0,$ (26)

and the conservation of energy. Since the circular frequency $\omega$ is
conserved, the nonlinear equation $\omega^{2}-gk^{B}\tanh{(k^{B}h^{B})}=0$
needs to be solved numerically for the wave number $k^{B}$ using a Newton
solver. The wave height is then found from the conservation of energy as

$H^{B}=H^{A}\sqrt{\frac{C_{g}^{A}}{C_{g}^{B}}},$ (27)

where $C_{g}=\frac{d\omega}{dk}$ is the group velocity depending on the wave
number and depth. The derivation of equation (27) is described in detail by
for example in [12]. The set-down can be solved similarly using the procedure
outlined in [28] using the radiation stress.

Stage 2. Matching linear and nonlinear theory: Due to increasing waveheight,
the nonlinear algorithm must be utilized on the final stretch of the shoaling
curve. The transition must be made in a region where both theories are valid,
and some general conditions on the transition depth are given in [45]. In
order to initialize the nonlinear solver at the transition depth, three
parameters need to be specified. We choose to match the wave height $H$ and
the set-down $\overline{\eta}$ with the linear theory. In addition one of the
parameters $\lambda$, $c$, $\nu$, $\overline{q}_{E}$ can be be matched. To do
this we choose one of these parameters and then keeping $H$ and
$\overline{\eta}$ constant, we find an elliptic modulus $m\in(0,1)$ which
leads to equality in the chosen parameter. In principle we can only match one
parameter, and it is not clear which one will be the most convenient. To
investigate this issue, we define the root problems

$\lambda^{lin}-\lambda^{nonlin}(m)=0,\hskip 28.45274pt\hskip
14.22636ptc^{lin}-c^{nonlin}(m)=0,$ $\nu^{lin}-\nu^{nonlin}(m)=0,\hskip
39.83368pt\overline{q}_{E}^{lin}-\overline{q}_{E}^{nonlin}(m)=0.$

Here the linear parameter is a fixed constant while the nonlinear quantity is
a function of $m$ given that $H$ and $\overline{\eta}$ are given. Plotting the
root problems we observe that the wavelength is most sensitive to varying $m$,
and therefore is the best candidate for a matching parameter. Figure 4 shows a
particular example.

Figure 4: Root problems defined by the parameters
$\lambda,c,\nu,\overline{q}_{E}$ as functions of $m$.

Stage 3. Cnoidal shoaling: The final step in the shoaling procedure is solving
for wave height and set-down using the scheme defined in Section 6. First,
define $H$ and $\overline{\eta}$ as functions of $m$ as given by formula (24)
and (25). Then use a nonlinear solver to find $m$ from equation (22). With $m$
in hand, one can determine the wave height $H(m)$ at the current local depth.
Repeating this procedure will then determine changes of wave height and set-
down at consecutive points, allowing us to move up the slope.

Figure 5: Profile of the mean water level $\overline{\eta}$ and the wave
height $H$ as functions of the depth, compared to data points taken from [9].
In this case, we have wave period $T=1.14$s incident wave height
$H_{0}=6.45$cm and breaking wave height $H_{b}=8.55$cm. In the experiment, the
beach slope is $\mathrm{tan}(\beta)=0.082$. The lower panel depicts the same
scenario in a different aspect ratio with a zoom-in around the still water
line (S.W.L.).

Having defined the shoaling equations and the numerical approach, we may plot
the development of the wave height and set-down for a practical shoaling
problem. A comparison of our numerical results with the wave tank data
obtained by Bowen et al. [9] and the classical linear theory is shown in
Figure 5. In this figure, the red circles represent the wave height of the
shoaling wave according to the experimental data of [9], and the red asterics
represent the experimentally determined set-down and set-up, i.e. the
deviation of the mean depth from the still water line (S.W.L.). The dashed
blue curves show the prediction of the linear theory as already presented in
[9]. Finally, the solid curves represent the shoaling model put forward in the
present article. The waves under consideration had an initial wave length of
$\lambda_{0}=202$cm and an initial wave height $H_{0}=6.45$cm before they
reached the toe of the slope. The beach had a slope of $1:12$.

The linear shoaling model agrees well with the experimental data up to
moderate wave heights. On the other hand, the benefit of the nonlinear
formulation can be clearly seen. Indeed, the nonlinear theory clearly yields a
better fit of both wave height and set-down close to the breaking point of the
wave. One interesting detail in the lower panel of Figure 5 is that the mean
depth appears to rise in the nonlinear theory while it continues on downwards
in the linear theory. This nonlinear mid-elevation rise was acknowledged in
[25] has been found previously with higher-order methods. For example, a
combination of third-order hyperbolic waves near the shore and Stokes waves
further offshore was used in [17], and Cokelet’s extension of Stokes’
approximation of periodic waves was used in [41]. In both works, a similar
rise of the nonlinear set-down was observed. In the present work, the computed
set-down matches the experimental data even well beyond the breaking point,
but we have to acknowledge that the KdV model ceases to be valid for breaking
waves.

##  6 Shoaling without set-down

It was argued in [43] that the wave set-down is negligible in many practical
cases. With this assumption, the shoaling problem can be simplified since only
two equations need to be solved. This approach has been taken in many works on
wave shoaling.

The particular aim of this section is to improve the method put forward in
[44]. As already mentioned, the use of a first-order approximation of the
energy flux [44] gave rise to a discontinuity in wave height at the matching
point between the linear and their nonlinear theory due to the inconsistency
of the approximation of the energy flux. A partial fix was presented in [45]
where a continuous shoaling curve is obtained by use of conversion tables but
at the cost of not conserving the energy flux. More recently, the shoaling
theory was extended in [20] using the nonlinear definition of the energy flux
developed in [3, 4]. However, due to the specific way the problem was
formulated, it was not possible to reach all the way up the shoaling curve due
to numerical instabilities.

In the present section, it will be shown that with the definition of the
shoaling equation from Section 4, and the implementation explained in Section
5, we can find a continuous shoaling curve that reaches all the way up to the
highest wave and agrees well with the more accurate methods of [36] using
Cokelet’s theory and that of [46] using the fully nonlinear steady Euler
equations.

The algorithm is the same as in Section 5, except that in equation (25) we set
$\overline{\eta}=0$ at each step. We then use the same three stages as
explained in Section 5 to determine the development of the wave height in the
shoaling region. First formulate $H=H(m)$ according to formula (24) given in
Section 6. Then use the assumption of zero set-down to find the roots given by
(20). Having the roots as functions of $m$ we may use energy conservation to
define a nonlinear equation as done in equation (22). Solving for $m$ we are
free to determine the wave height at a specified depth $h$.

Figure 6: Shoaling curves based on present theory in blue compared to the
theory presented by Svendsen and Brink-Kjær [44] in red. Deep water values
$H_{0}/\lambda_{0}=\\{0.001,0.002,0.004,0.006,0.01\\}$.

We note that the present implementation is able to determine the wave height
further into the shoaling region as compared to the curves presented by [20].
This is due to the simplicity of the implementation, solving two nonlinear
equations rather than three. We also see in Figure 6 that our theory is in
fairly good agreement with the shoaling curves presented in [44] which are in
reasonable agreement with higher-order theories. Moreover, it is clearly
visible in Figure 6 that our curves do not feature a discontinuity in wave
height.

Finally, comparisons of shoaling curves computed using the method at hand with
data provided by wave tank experiments are presented. The experimental data
are taken from [45] where considerable pains were taken to stay within the
correct scaling regime. The waves were generated with a piston-type wavemaker,
then propagated over a flat bottom with still water depth of $36$ cm before
shoaling on a $1:35$ beach.

Figure 7: Comparisons of shoaling curves from six experiments with a range of
initial wave steepness.

We initialize the code with wave height and wavelength at depth $h=36$cm at
the toe of the slope. The code then computes wave height $H$ with the linear
model up to $h/\lambda_{0}>0.1$ at which time the code switches to the cnoidal
theory.

We observe in Figure 7 that there is very good agreement between the numerical
model (shown in blue) and the experimental data (indicated by red crosses).
The only exception is the experiment with deep water steepness
$H_{0}/\lambda_{0}=0.064$, shown in the uppermost panel in Figure 7. This is a
rather steep wave and it can be observed that even the linear theory fails to
get good agreement. This error is then propagated and compounded by the
nonlinear theory. On the other hand, as long as the deep-water steepness is
moderate enough, the linear theory does an adequate job, and we observe fairly
good agreement in the plots between experimental data and numerical
simulations for both linear and nonlinear regimes. One should note that the
higher-order theory by Cokelet was used numerically for the same data set in
[36], but no better agreement was found than with the present theory. We have
also done some preliminary computations with a full Boussinesq model [35], and
similar agreement with the wavetank data was found after some experimentation
with grid size and friction coefficients.

##  7 Conclusion

In this paper we have shown how to derive expressions for the energy flux and
radiation stress of a wavetrain in the context of the well known KdV equation.
The derivations are based on a general approach developed in [3, 4], and yield
approximations of the same order as the KdV equation itself. Compared to
previous definitions of such quantities as for example in [43, 44], these
expressions feature additional terms which are required by the theoretical
underpinnings of the method used here.

The energy flux and radiation stress are used to develop shoaling equations
for long waves of small to moderate amplitude, and since such waves can be
described by cnoidal functions, the shoaling equations can be posed as a
$3\times 3$ system of equations for the three parameters of the cnoidal
functions. The shoaling equations are then formulated in a numerically
convenient way, and it is shown that they can be solved to yield shoaling
curves up to the breaking point. In particular, our formulation resolves a
problem encountered in earlier work [20] where the curves terminated
prematurely due to numerical instability. The curves are compared with
experimental data and higher-order theories, and are found to be quite
accurate. In addition, since the shoaling equations are defined using
approximations of energy flux and radiation stress which are consistent with
the approximations made in the KdV equation itself, the shoaling curves are
continuous, and no ad-hoc fixes such as advocated for in [45] are necessary to
get a continuous development of the wave height,

Finally, since the radiation stress is included in the formulation of the
shoaling equations, it is possible to compute the set-down as part of the
shoaling problem. Comparisons with experimental data from [9] show that the
set-down can be computed accurately. The nonlinear set-down defined here may
also potentially be useful in the study of waves interacting with currents
[43].

To put our study into context, note that due to the ready availability of
computational resources, shoaling is nowadays usually computed by utilizing
nearshore wave models which are used to forecast wave conditions in the
coastal zone. These models are generally known as Boussinesq models, and
probably the first such model was developed by Peregrine [33]. These models
have become fairly sophisticated in recent years, being able to treat even
fairly short waves due to extending the dispersion properties using ideas of
Nwogu [31] and Witting [49]. Some of the models currently in use are described
in [10, 30, 35]. In particular, as shown for example in [29], shoaling can be
computed accurately with such models. Following ideas of [38], some models
have been extended to be able to treat larger-amplitude waves, and it is now
possible to simulate waves with fully nonlinear and highly dispersive models
[23, 47] or with the full Euler equations [13]. Nevertheless, for practical
purposes, these models need to be initialized with data which are generally
given by stochastic wave models, and in general the transition between
stochastic and deterministic models is still poorly understood. The method of
computing wave height and set-down put forward here may have some potential if
coupled with stochastic input data (see [14]) since the computational
complexity of our method is by far smaller than for any phase-resolving
nearshore model and no tuning is necessary with the current model.

## Acknowledgments

The authors thank Volker Roeber for helpful discussions. This work has
received funding from the European Union’s Horizon 2020 research and
innovation programme under grant agreement No. 763959.

##  Appendix A Integrals of cnoidal functions

In order to define the two shoaling models we need to handle integrals of
various powers of derivatives of $\eta$. In particular we need to determine
$\overline{\eta}$, $\overline{\eta^{2}}$, $\overline{\eta^{3}}$,
$\overline{\eta_{xx}}$ and $\overline{\eta\eta_{xx}}$ in order to evaluate
(25) in terms of $m$ and the current depth $h$. First note that we can write
the time-averaged $\eta$ given by (18) in the more convenient form

$\overline{\eta}=\frac{1}{T}\int_{0}^{T}\eta\Big{(}2K(m)(\frac{t}{T}-\frac{x}{\lambda})\Big{)}dt=f_{2}+H\int_{0}^{1}\text{cn}^{2}(2K(m)\xi;m)d\xi,$

(for more details see [43]). Similarly considerations apply to the powers of
$\eta$ which involve integrals of the form

$\displaystyle\int_{0}^{1}\text{cn}^{2}(2K\xi;m)d\xi=\frac{1}{4mK}\Big{(}E-(1-m)K\Big{)},$
$\displaystyle\int_{0}^{1}\text{cn}^{4}(2K\xi;m)d\xi=\frac{1}{3m^{2}}\Big{(}3m^{2}-5m+2(4m-2)\frac{E}{K}\Big{)},$
$\displaystyle\int_{0}^{1}\text{cn}^{6}(2K\xi;m)d\xi=\frac{1}{5m^{2}}\Big{(}4(2m^{2}-1)\int_{0}^{1}\text{cn}^{4}(2K\xi;m)d\xi+3(1-m^{2})\int_{0}^{1}\text{cn}^{2}(2K\xi;m)d\xi\Big{)}.$

These expressions can be found in [24]. The terms $\overline{\eta_{xx}}$ and
$\overline{\eta\eta_{xx}}$ also include such terms due to the identity

$\eta_{xx}(2K\xi;m)=4K^{2}H(2-2m+(8m-4)\text{cn}^{2}(2K\xi;m)-6m\text{cn}^{4}(2K\xi;m)),$

and similar relations which can be found in [1]. Combining the results above
we deduce that

$\displaystyle\overline{\eta}\hskip
14.22636pt=f_{2}+H\int_{0}^{1}\text{cn}^{2}(2K\xi;m)d\xi,$
$\displaystyle\overline{\eta^{2}}\hskip
10.243pt=f_{2}^{2}+2Hf_{2}\int_{0}^{1}\text{cn}^{2}(2K\xi;m)d\xi+H^{2}\int_{0}^{1}\text{cn}^{4}(2K\xi;m)d\xi,$
$\displaystyle\overline{\eta^{3}}\hskip
10.243pt=H^{3}\int_{0}^{1}\text{cn}^{6}(2K\xi;m)d\xi+3H^{2}f_{2}\int_{0}^{1}\text{cn}^{4}(2K\xi;m)d\xi+3Hf_{2}^{2}\int_{0}^{1}\text{cn}^{2}(2K\xi;m)d\xi+f_{2}^{3},$
$\displaystyle\overline{\eta_{xx}}\hskip
5.69054pt=\frac{3H^{2}}{2m}(-3m^{2}\int_{0}^{1}\text{cn}^{4}(2K\xi;m)d\xi+4m^{2}\int_{0}^{1}\text{cn}^{2}(2K\xi;m)d\xi-2\int_{0}^{1}\text{cn}^{2}(2K\xi;m)d\xi-m^{2}+1),$
$\displaystyle\overline{\eta\eta_{xx}}=3H^{2}f_{3}\int_{0}^{1}\text{cn}^{4}(2K\xi;m)d\xi-3H^{2}f_{1}\int_{0}^{1}\text{cn}^{4}(2K\xi;m)d\xi+\frac{3}{2}Hf_{1}f_{2}-\frac{3}{2}Hf_{2}f_{3}$
$\displaystyle\hskip
28.45274pt+\frac{3}{2}H^{2}f_{1}\int_{0}^{1}\text{cn}^{2}(2K\xi;m)d\xi-\frac{3}{2}H^{2}f_{3}\int_{0}^{1}\text{cn}^{2}(2K\xi;m)d\xi+6H^{2}f_{1}m^{2}\int_{0}^{1}\text{cn}^{4}(2K\xi;m)d\xi$
$\displaystyle\hskip
28.45274pt-6H^{2}f_{3}m^{2}\int_{0}^{1}\text{cn}^{4}(2K\xi;m)d\xi-\frac{3}{2}Hf_{1}f_{2}m^{2}+\frac{3}{2}3Hf_{2}f_{3}m^{2}-\frac{3}{2}H^{2}f_{1}m^{2}\int_{0}^{1}\text{cn}^{2}(2K\xi;m)d\xi$
$\displaystyle\hskip
28.45274pt+\frac{3}{2}H^{2}f_{3}m^{2}\int_{0}^{1}\text{cn}^{2}(2K\xi;m)d\xi-3Hf_{1}f_{2}\int_{0}^{1}\text{cn}^{2}(2K\xi;m)d\xi+3Hf_{2}f_{3}\int_{0}^{1}\text{cn}^{2}(2K\xi;m)d\xi$
$\displaystyle\hskip
28.45274pt-\frac{9}{2}H^{2}f_{1}m^{2}\int_{0}^{1}\text{cn}^{6}(2K\xi;m)d\xi+\frac{9}{2}H^{2}f_{3}m^{2}\int_{0}^{1}\text{cn}^{6}(2K\xi;m)d\xi$
$\displaystyle\hskip
28.45274pt+6Hf_{1}f_{2}m^{2}\int_{0}^{1}\text{cn}^{2}(2K\xi;m)d\xi-6Hf_{2}f_{3}m^{2}\int_{0}^{1}\text{cn}^{2}(2K\xi;m)d\xi$
$\displaystyle\hskip
28.45274pt-\frac{9}{2}Hf_{1}f_{2}m^{2}\int_{0}^{1}\text{cn}^{4}(2K\xi;m)d\xi+\frac{9}{2}Hf_{2}f_{3}m^{2}\int_{0}^{1}\text{cn}^{4}(2K\xi;m)d\xi.$

## References

  * [1] M. Abramowitz and I.A. Stegun, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, National Bureau of Standards Applied Mathematics series 55, (1965).
  * [2] A. Ali and H. Kalisch, Energy balance for undular bores, Compt. Rend. Mecanique 338 (2010), 67–70.
  * [3] A. Ali and H. Kalisch, Mechanical balance laws for Boussinesq models of surface water waves, J. Nonlinear Sci. 22 (2012), 371–398.
  * [4] A. Ali and H. Kalisch, On the formulation of mass, momentum and energy conservation in the KdV equation, Acta Applicandae Mathematicae 133 (2014), 113–131.
  * [5] J.A. Battjes and H.W. Groenendijk, Wave height distribution on shallow foreshores, Coastal Engineering 40 (2000), 161–182. no. 3
  * [6] B.T. Benjamin and J. Lighthill, On cnoidal waves and bores, Proc. R. Soc. A 224 (1954), 448–460.
  * [7] M. Bjørkavåg and H. Kalisch, Wave breaking in Boussinesq models for undular bores, Phys. Lett. A 375 (2011), 1570–1578.
  * [8] J. Boussinesq, Théorie des ondes et des remous qui se propagent le long d’un canal rectangulaire horizontal, en communiquant au liquide contenu dans ce canal des vitesses sensiblement pareilles de la surface au fond, J. Math. Pures Appl. 17 (1872), 55–108.
  * [9] A.J. Bowen, D.L. Inman and V.P. Simmons, Wave ‘set-down’ and set-up, J. Geophysical Research 73 (1968), 2569–2577.
  * [10] M. Brocchini, A reasoned overview on Boussinesq-type models: the interplay between physics, mathematics and numerics, Proc. R. Soc. A 469 (2013), 1–27.
  * [11] J. Buhr Hansen and I.A. Svendsen, Laboratory generation of waves of constant form, Proc. Coast. Eng. Conf., 14th, Cope, 1974.
  * [12] R.G. Dean and R.A. Dalrymple, Water Wave Mechanics for Scientists and Engineers, (World Scientific, Singapore, 1984).
  * [13] S. Grilli, R. Subramanya, I.A. Svendsen and J. Veeramony, Shoaling of solitary waves on plane beaches, J. Wtrwy., Port, Coast., and Oc. Engrg. 120 (1994), 609–628.
  * [14] L.H. Holthuijsen, Waves in Oceanic and Coastal Waters. (Cambridge University Press, Cambridge, 2010).
  * [15] S. Israwi and H. Kalisch, Approximate conservation laws in the KdV equation, Phys. Lett. A 383 (2019), 854–858.
  * [16] S. Israwi, H. Kalisch, A mathematical justification of the momentum density function associated to the KdV equation, in press.
  * [17] I.D. James, 1974. Non-linear waves in the nearshore region: shoaling and set-up, Estuarine Coastal Marine Science 2 (1974), 207–234.
  * [18] D.J. Korteweg and G. de Vries, On the change of form of long waves advancing in a rectangular channel and on a new type of long stationary wave, Philos. Mag. 39 (1895), 422–443.
  * [19] Z. Khorsand, Flow Properties of Fully Nonlinear Model Equations for Surface Waves, The University of Bergen, (2017).
  * [20] Z. Khorsand and H. Kalisch, On the shoaling of solitary waves in the KdV equation, Coastal Engineering Proceedings, 34(2014).
  * [21] P. Kundu and I. Cohen, Fluid Mechanics (Academic Press, San Diego, 2004).
  * [22] D. Lannes, The Water Waves Problem. Mathematical Surveys and Monographs, vol. 188 (Amer. Math. Soc., Providence, 2013).
  * [23] D. Lannes and P. Bonneton, Derivation of asymptotic two-dimensional time-dependent equations for surface water wave propagation, Phys. Fluids 21 (2009), 016601.
  * [24] D.F. Lawden, Elliptic Functions and Applications (Springer, New York, 2013).
  * [25] B. Le Méhauté, An Introduction to Hydrodynamics and Water Waves (Springer, New York, 1976).
  * [26] J. Lighthill, Waves in Fluids (Cambridge University Press, Cambridge, 1978).
  * [27] M.S. Longuet-Higgins and R.W. Stewart, Radiation stress and mass transport in gravity waves, with application to ‘surf beats’, J. Fluid Mech. 13 (1962), 481–504.
  * [28] M.S. Longuet-Higgins and R.W. Stewart, Radiation stresses in water waves; a physical discussion, with applications, Deep Sea Research 11 (1964), 529–562.
  * [29] P.A. Madsen D.R. Fuhrman and B. Wang, A Boussinesq-type method for fully nonlinear waves interacting with rapidly varying bathymetry, Coast. Eng. 53 (2006), 487–504.
  * [30] P.A. Madsen and H.A. Schäffer, Higher-order Boussinesq-type equations for surface gravity waves: derivation and analysis, Phil. Trans. Roy. Soc. Lond. Ser. A 356 (1998), 3123–3184.
  * [31] O. Nwogu, Alternative form of Boussinesq equations for nearshore wave propagation, J. Waterway, Port, Coastal and Ocean Engineering 119 (1993), 618–638.
  * [32] L.A. Ostrovsky and E.N. Pelinovsky, Wave transformation on the surface of a fluid of variable depth, Atmospheric and Oceanic Physics 6 (1970), 552–555.
  * [33] D.H. Peregrine, Long waves on a beach, J. Fluid Mech. 27 (1967), 815–827.
  * [34] Lord Rayleigh, Hydrodynamical notes, Phil. Mag. 21 (1911), 177–195.
  * [35] V. Roeber, K.F. Cheung and M.H. Kobayashi, Shock-capturing Boussinesq-type model for nearshore wave processes, Coastal Engineering 57 (2010), 407–433.
  * [36] T. Sakai and J.A. Battjes, Wave shoaling calculated from Cokelet’s theory, Coastal Engineering, 4 (1980), 65–84.
  * [37] T. Saville, Experimental determination of wave set-up, Proc. 2nd Tech. Conf. on Hurricanes, Miami Beach, FL, (1961).
  * [38] F. Serre, Contribution á l’ étude des écoulements permanents et variables dans les canaux, La Houille Blanche, 8 (1953), 374–388.
  * [39] O. Skovgaard, Sinusoidal and cnoidal gravity waves formulae and tables, Institute of Hydrodynamics and Hydraulic Engineering, Technical University of Denmark, (1974).
  * [40] R.M. Sorensen, Basic Wave Mechanics for Coastal and Ocean Engineers (John Wiley & Sons, New York, 1993).
  * [41] M. Stiassnie and D.H. Peregrine, Shoaling of finite-amplitude surface waves on water of slowly-varying depth, J. Fluid Mech. 97 (1980), 783–805.
  * [42] J.J. Stoker, Water Waves: The Mathematical Theory with Applications, (Interscience Publishers, New York, 1957).
  * [43] I.A. Svendsen, Introduction to Nearshore Hydrodynamic, (World Scientific, 2006).
  * [44] I.A. Svendsen and O. Brink-Kjær, Shoaling of cnoidal waves, Proc. 13th Intern. Conf. Coastal Engrg., Vancouver (1972), 365–383.
  * [45] I.A. Svendsen and J. Buhr Hansen, The wave height variation for regular waves in shoaling water, Coastal Engineering, 1 (1977), 261–284.
  * [46] R.H. Swift and J.C. Dixon, Transformation of regular waves, Proceedings of the Institution of Civil Engineers 83 (1987), 359–380.
  * [47] G. Wei, J.T. Kirby, S.T. Grilli and R. Subramanya, A fully nonlinear Boussinesq model for surface waves. Part 1. Highly nonlinear unsteady waves, J. Fluid Mech. 294 (1995), 71–92.
  * [48] G. B. Whitham, Linear and Nonlinear Waves (Wiley, New York, 1974).
  * [49] J.M. Witting, A unified model for the evolution of nonlinear water waves, J. Comp. Phys. 56 (1984), 203–236.