\section{Introduction}

It has now become our basic
knowledge that the extratropical
atmosphere is driven strongly by the horizontal
potential temperature gradient that arises from the
differential solar heating (e.g., Lorenz 1955). The horizontal
gradient in the potential temperature, often referred
to as baroclinicity, is a measure of upper-level wind
steering via thermal wind and a measure of baroclinic
instability in the atmosphere. Baroclinicity in the lower
atmosphere in classic theories of atmospheric stability is
measured by a combination of the static stability and
horizontal temperature gradient, the latter of which is
equivalent to vertical shear in the horizontal wind through
the thermal wind balance (Charney 1947; Eady 1949). In
its original form, the Eady's maximum growth rate for
baroclinic instability $B_{GRMax}$ is defined by $B_{GRMax} =
0.31(|f|/N)(\partial U/\partial_z)$ in a zonally homogeneous steady
mean state, where $U$ is the mean zonal flow, $f$ is the
Coriolis parameter, and $N$ is the Brunt--V\"ais\"al\"a frequency.
Charney's formula is slightly different from the Eady's,
but still incorporates the same effects.

Lindzen and Farrell (1980) first applied the Eady's parameter
to atmospheric data to successfully estimate the
maximum growth rate of baroclinic instability in the troposphere.
Hoskins and Valdes (1990) used its localized
version (i.e., $U$, $N$, and $f$ are all local Eulerian mean values)
as the central parameter in their study of the Northern
Hemispheric storm tracks. This local version, or its simplified
version, has been used successfully as an indicator
of baroclinic wave generation in diagnostic studies of
stormtracks in recent years as well (Nakamura and Sampe
2002; Nakamura and Shimpo 2004; Nakamura et al. 2004).
In our study, the North Atlantic part of which was reported
in Nakamura and Yamane (2009, hereafter Part I), we
define the near-surface baroclinic vector, $\mathbf{B}=B^x\mathbf{i} +
B^y\mathbf{j}$,
where $B^x= -(g/\theta N)(\partial\theta/\partial y)$
and $B^y = (g/\theta N)(\partial\theta/\partial x)$  with
$\theta$ being the monthly-mean potential temperature at 2 m
above the surface, and use it as the central quantity of the
diagnoses. Unless stated otherwise, ``anomalies'' refer to
deviations from the climatology hereafter. Though its meridional
component does not appear in any classic theory of
baroclinic instability, a theory that does incorporate the
effect of $B^y$ shows its important role in enhancing baroclinic
wave generation locally to the east of the mean trough
(Niehaus 1980). In the North Atlantic storm track region,
we indeed found that the substantial zonal gradient in the
surface temperature in and around the Labrador Sea plays
a major role in the large-scale atmospheric state.

\section{The sea surface temperature}
The sea surface temperature
 (SST) is an important factor
in determining $\mathbf{B}$ in the storm-track regions (e.g.,
Hoskins and Valdes 1990; Nakamura et al. 2004; Part I).
SST anomalies (SSTAs) around an oceanic front along the
Gulf Stream (GS), Kuroshio Extension (KE), or Oyashio
Extension (OE) can have a profound impact on $\mathbf{B}$ along
the storm tracks.

\subsection{A subtle point}
A subtle but important point that has to
be considered carefully in this regard is the spatial scale
and the location of SSTAs with respect to the climatology,
since it is the anomalous surface temperature
gradient whose structure has a spatial scale of the atmospheric
Rossby deformation radius that can exert significant
influence on the large-scale atmospheric flow.

\subsubsection{Changes in the temperature}
The
high sensitivity of $\mathbf{B}$ to changes in the temperature contrast
across the front and changes in the width of the front, and
the uncertainty in the impact of SSTAs of small spatial
scales on $\mathbf{B}$ make it difficult to assess the effective B
anomalies that are attributable to the SSTAs from the
available data. Moreover, it is uncertain exactly how the
SSTAs in the presence or absence of the land surface
temperature anomalies may or may not produce $\mathbf{B}$ anomalies
that are significant to the atmosphere.

\paragraph{Complicating factor}The complicating
factor introduced by the land surface must be taken
into account when studying potential roles of extratropical
SSTAs in the extratropical atmospheric anomalies.

Lau (1988) investigated patterns of anomalous storm
track activity and associated low-frequency flow anomalies
by computing empirical orthogonal functions (EOFs)
for high-frequency 500-hPa geopotential height for the
Northern Hemisphere winters. He found that both North
Atlantic and North Pacific storm tracks have a pattern of
meridional shift and a pattern of increased or decreased
eddy activity in the first two EOFs.He also found that these
changes in the storm tracks have symbiotic relationships
with the background flows and have substantial impacts on
the mean flow. Part I approached the issue of the stormtrack
and low-frequency flow variability in connection with
SSTAs in the extratropics, focusing on $\mathbf{B}$ as the key parameter
of diagnoses, and found similar patterns of variability
in the eddy activity and low-frequency flow in the
North Atlantic basin. Much of this variability was connected
to SSTAs in the vicinity of the Gulf Stream in cold
months (Part I). Since the winter North Pacific basin has a
storm track andmean flow that appear to be related to the
oceanic fronts, Kuroshio--Oyashio Extensions (KOE) in
this case, in a manner essentially the same as those in the
North Atlantic basin related to the GS,we have attempted
to find similar results for the North Pacific basin. In this
regard, we have chosen not to project our results onto the
major mode of variability in the extratropical North Pacific
basin, the North Pacific decadal variability (PDV), so
that our presentation and discussion are mostly confined
to the wave--mean flow dynamics of monthly time scale or
shorter.

Our approach to the search for a link between anomalies
in the KOE and the overlying atmosphere is as follows:
(i) identify dominant patterns in anomalous $\mathbf{B}$ in the
storm track for each calendarmonth and identify years in
which the anomaly fits the pattern well, (ii) composite
anomalies in the monthly-mean circulation and highfrequency
transients in the atmosphere to obtain a typical
atmospheric state that accompanies the patterns of
anomalous $\mathbf B$, (iii) composite SSTA to obtain a typical
oceanic state that accompanies and precedes the patterns
of anomalous $\mathbf B$, and (iv) composite anomalous net surface
heat flux that accompanies and precedes the pattern of
anomalous $\mathbf B$. With this approach, we obtain typical pictures
of anomalous states in the atmosphere and oceans
with anomalous $\mathbf{B}$ as their connecting interface.
Section 2 describes the data and procedure to compute
$\mathbf B$. Section 3 describes the climatology and variance of $\mathbf B$.
Dominant patterns of $\mathbf B$ are shown in section 4, followed by
composited anomalies in various atmospheric fields and
SST in section 5. Finally, we present our discussion on the
results, examining a potential cause--effect relationship
between anomalies in the SST and atmosphere in section 6.

Figure \ref{fig1} shows the climatology of $Bx$, $U^{200}$,
$\overline{V'\theta'}^{850}$, $\overline{V'\theta'}^{200}$,
$MR^{z850}$, $U^{1000}$, SST, and $F_h$ for February and
August as examples of the reference state in the winter
and summer. The numeric superscript indicates the pressure
level in hPa. The seasonal mean is visibly more
diffused in structure, particularly for Bx, than those
shown in Fig.~\ref{fig1}. There are no surprises in the overall
picture of Bx. The regions of large land--sea temperature
contrast and the oceanic fronts show very large Bx in cold
months. The position of Bx maximum in the storm track
region is found along theKOEregion throughout the year.
The zonally elongated band of large Bx in the storm track
is generally wider in cold months than in warm months.
In fact, large Bx values that presumably accompany the
Kuroshio and its extension in the cold months vanish in
the summer (Figs. 1a,b). The seasonal variation in the Bx
values in the storm track basically follows that of the
north--south differential heating--largest in the winter and
smallest in the summer. As reported by Nakamura et al.
(2004), the $U^{200}$ maximum is displaced southward from the
Bx maximum visibly in winter months (Fig.~1a), although
U200 is generally large over the area of large Bx in the core
of the storm track. The southward displacement of the
U200 maximum from the Bx maximum in the storm track
is less pronounced or even reversed in warmer months
(Fig.~1b). We note that the structure of the climatological
Bx generally reflects those of ($\partial{\theta^{2{\textrm m}}}/\partial{y}$) with some exceptions where
the surface slope contributes significantly to Bx
in isolated areas over the land, most notably around the
Himalayas.

The position and structure of the storm track as indicated
by $\overline{V'\theta'}^{850}$
and $\overline{V'V'}^{200}$
are, at least for the winter,
essentially the same as those reported in earlier studies on
the storm tracks (e.g., Chang et al. 2002). The maxima in
$\overline{V'\theta'}^{850}$
and $\overline{V'V'}^{200}$ are located in the band of large Bx.

\section{Data and calculation procedures}
The data used to calculate $\mathbf{B}$ are the monthly-mean
temperature at 2 m above the surface ($T^{2m}$) and temperature
at pressure levels available from the 40-yr European
Centre for Medium-Range Weather Forecasts (ECMWF)
Re-Analysis (ERA-40; Uppala et al.~2005). We chose
the ERA-40 $T^{2m}$ data rather than the National Centers
for Environmental Prediction--National Center for Atmospheric
Research (NCEP--NCAR) reanalysis products
for its explicit inclusion of the observed near-surface temperature
in producing the $T^{2m}$ data. The monthly-mean
surface pressure data from the NCEP--NCAR reanalyses
(Kalnay et al. 1996) were used to determine the pressure
levels to be used for $\mathbf B$ calculation, and to calculate $\theta$
at 2 m
above the surface from $T^{2m}$. We used the NCEP--NCAR
surface pressure data for convenience, since we had already
compiled the dataset for calculating transient eddy
fluxes to be mentioned later and the ERA-40 surface
pressure data are not readily available. We later compared
the NCEP--NCAR monthly-mean sea level pressure with
that of ERA-40, and found the difference between the two
products to be immaterial for the purpose of the current
study. We also used ERA-40 monthly-mean horizontal
wind and geopotential height at pressure levels, net surface
heat flux, Fh (the sum of latent heat flux, sensible heat flux,
solar radiation, and the thermal radiation), and Hadley
Centre sea surface temperature data (Rayner et al. 2003)
to compile anomaly composites accompanying anomalous
patterns in $\mathbf B$. In addition, we used 6-hourly temperature
and wind data from the NCEP--NCAR reanalyses to
compute various eddy fields. The accuracy of the Fh data
used here is, as true for other reanalyses surface heat flux
products, may not be so high to produce reliable anomaly
composites.

We computed $\mathbf{B}$ near the surface by calculating the
horizontal gradient in $\theta^{2{\textrm m}}$, using the centered finite differencing,
and calculating $N$ from the lowest three vertical
pressure levels that are location dependent because of
topography. Both $\nabla\theta^{2{\textrm m}}$ and $N$ were calculated locally as in
Hoskins and Valdes (1990) and Nakamura and Shimpo
(2004). The entire 45 yr from September 1957 to August
2002 were used for the Northern Hemisphere. To resolve
the dominant modes in $\mathbf{B}$ arising from the land--sea temperature
contrast, one may need much higher horizontal
resolution in the data. The relatively coarse horizontal
resolution of the data may artificially suppress the significance
of the variability associated with the land--sea temperature
contrast.One should keep this limitation in mind.
The 6-hourly bandpassed (period of 2--7 days) eddy
fields and ultra-low-frequency (period of 30 days and longer)
background fields were computed from the NCEP--
NCAR reanalyses, using simple time filters (Lau and Lau
1984) first. The filtered time series were then visually examined
against the raw time series and, then, used to calculate
the slowly evolving bandpassed meridional velocity
variance ($\overline{V'V'}$), meridional temperature flux ($\overline{V'\theta'}$), and
the three-dimensional transient wave activity flux defined
on a zonally varying basic state by Plumb (1986). The wave
activity flux consists of the zonal and meridional advective
fluxes (MU and MV), the zonal and meridional
radiative fluxes (MRx and MRy), and the radiative vertical
flux (MRz). The flux is essentially the Eliassen--Palm
flux (Eliassen and Palm 1961) in a zonally inhomogeneous
mean flow (Plumb 1986). The wave activity flux was calculated
from February 1948 to November 2004 only for
the extratropics poleward of 208 latitude. Also, it was
calculated only from 850 to 30 hPa because of the double
differentiation with respect to pressure required for the
calculation. The flux of particular interest in this study is
the vertical component. Here MR$^z$ is defined by
\[
MR^z=
\frac{
pf \cos\phi}
{p^0|\nabla h\bar q|(d\theta_0/dz)}
\left(\frac{\partial\bar q}{\partial x}
\overline{U'\theta'} +
\frac{\partial\bar q}{\partial y}
\overline{V'\theta'}\right),
\]
where $q$ is quasigeostrophic potential vorticity, $p$ is the
pressure, $p^0$ is the reference surface pressure set to
1000 hPa here, $\theta_0$ is the area-weighted ultra-low-frequency
hemispheric mean potential temperature at each height,
u9 is the bandpassed potential temperature, and z is the
geopotential height. An overbar denotes an ultra-low frequency
component and a prime denotes bandpassed
component. The 6-hourly time series of wave fluxes was
computed by using the time series of ultra-low-frequency
fields as the basic-state and high-frequency fields as eddies.
In short, the time series was calculated by changing the
meaning of an overbar from the time mean state to an
ultra-low-frequency state, and changing the meaning of a
prime from a departure from the mean to a high-frequency
state. The 6-hourly eddy time series was averaged over each
month to produce monthly-mean time series. This dataset
allows us to examine anomalous eddy fields accompanying
anomalous $\mathbf B$ in specific months. The climatology for
the eddy fields was computed from 46 yr, January 1958 to
December 2003. The calculation of the wave activity and
its flux is described in detail by Nakamura et al. (2010).

\section{Climatology and variance}
The climatology and variance of $Bx$ and $By$ were computed
for each calendar month and examined closely for
their spatial and temporal structures. The monthly climatology,
rather than the seasonal climatology, is used as the
reference in our study, to avoid contamination of the diagnostic
results arising from differences in the climatology
between two successive months. Unlike in the North Atlantic
basin reported in Part I, we found the impact of $By$
variations on the large-scale atmospheric state in the
North Pacific basin much weaker than that of Bx variations,
presumably because of themore zonal orientation of
the KOE in comparison to the Gulf Stream and North
Atlantic Current. In the following, thus, we focus our
presentation on the climatology and variations of $Bx$ and
their impact on the large-scale atmospheric state.

Sample citations: \citet{Becker+Schmitz2003}, \citet{Knutti2008},
and \citep{MeixnerEA2002,Kuji_Nakajima2002,EmeryEA1986}.

\acknowledgments
We thank two anonymous reviewers for their comments, which helped
to improve the manuscript.

\datastatement
Start data availability statement here..

\section{Introduction}
Oceanic internal tides (ITs) are waves generated in a stratified ocean (internal waves) through the interaction of the astronomically induced barotropic tidal flow with the seafloor. They are ubiquitous oscillatory baroclinic perturbations of the fluid velocity and density at the tidal frequency, which propagate away from bottom irregularities such as seamounts, continental shelves, or ridges \citep{Garrett_Kunze_2007}. ITs are considered to be one of the main sinks of energy for the barotropic tide, and an important contributor to ocean mixing at a global scale \citep{Garrett2003, Wunsch_Ferrari, 2016_JPO_deLavergneMLNN}. An accurate description of the IT flow and a good prediction of barotropic-to-baroclinic energy conversion are necessary for the reliable ocean and climate modelling. Climate models often implement internal-wave driven mixing through a parametric formula involving an idealized estimation of the local energy conversion rate at generation sites. Moreover, the dissipation of ITs is also parametrized by using a modal description of the flow. We refer to \citet{CMP2017} for a detailed overview of these considerations.
Various analytical or semi-analytical approaches have been developed in the past decades to model ITs, mostly in two-dimensions. As discussed in \citet{Garrett_Kunze_2007}, these approaches may be classified in two categories in terms of the type of the barotropic forcing: $(i)$ the boundary-forcing approach and $(ii)$ the body-forcing approach. Topographic features are commonly classified into two regimes, subcritical or supercritical, depending on whether the maximum slope value is smaller or larger than the characteristic slope of internal waves at the tidal frequency. Both of the aforementioned approaches can be used to study ITs generated in any of the two topographical regimes.

In the boundary-forcing approach,  a homogeneous internal wave equation is solved with a non-homogeneous bottom boundary condition, which implies an oscillating barotropic velocity of constant amplitude in the vertical direction. It was first developed in conjunction with either of two strong and opposing assumptions. The first one is the ``weak topography approximation" (WTA) which allows the linearization of the bottom boundary condition around a flat bottom and therefore  approximates the flow over a topography of small height but arbitrary horizontal shape \citep{Bell75, Llewellyn2002, Khatiwala2003, Grisouard2012}.
The second one is the  ``knife-edge" approximation \citep{Robinson1969, llewellyn_smith_young_2003}, where the topography is assumed of zero elevation everywhere except at a single point, where it is arbitrarily large. These idealizations make the problem analytically tractable, i.e., analytical methods such as Fourier/eigenfunction expansions and Green's function methods can lead to useful expressions for the far-field energy flux or conversion rate \citep{Llewellyn2002, llewellyn_smith_young_2003, StLaurent2003, Nycander2005, Falahat2014}.  Semi-analytical methods that consider the full bottom boundary condition have also been developed based on the Fourier representation of a periodic topography in an infinitely deep ocean \citep{Balmforth_Ierley_Young_2002} and on the Green's function method \citep{Petrelis_et_al_2006, Balmforth_Neil_Peacock_2009, echeverri_peacock_2010, iTides, Mathur2016}.

In the body-forcing formulation, the governing wave equation has a homogeneous bottom boundary condition but is forced by a source term that represents a non-uniform or topographically adjusted barotropic flow. This formulation was introduced by \cite{Baines1973}, who also proposed a technique based on ray theory and an integral-equation reformulation for the calculation of the far-field energy flux \citep[see also][]{BAINES1982}. \cite{Garett_Gerkemma_2007} noted that the body forcing term in  \cite{Baines1973} is inconsistent with non-hydrostatic conditions and derived a consistent formulation.  \cite{maas_2011} used appropriate transformations of the underlying equations, establishing the existence of a class of topographies that do not lead to far-field internal-tide generation (non-radiating topographies). To the best of our knowledge, general solution procedures that take into account the non-uniformity of the barotropic flow over varying topography are less developed in comparison with the uniform barotropic forcing case, a gap we propose to address here.

In this work, we adopt the body-forcing formulation and develop a semi-analytical model that describes the internal-tide generation process. The analytical step is the exact reformulation of the problem as a Coupled-Mode System (CMS) of equations that can be solved numerically in an efficient way. This is accomplished by means of a modal or eigenfunction expansion of the unknown stream function, in terms of appropriate basis functions and unknown modal amplitudes. This approach, also called coupled-mode theory, has been applied to the solution of various non-uniform waveguide problems in acoustics \citep{Brekhovskikh, Desaubies, Maurel}, geophysics \citep{MAUPIN1988, 8805310}, water waves \citep{porter_staziker_1995, athanassoulis_belibassakis_1999, PAPOUTSELLIS2018199} among other disciplines. In the context of IT-generation within the body-forcing formulation, \cite{griffiths2007} derived a CMS from the velocity formulation of the problem, with the hydrostatic approximation, and calculated 2D ITs over continental shelves. The CMS we develop here differs from this work in two ways. First we remove the hydrostatic approximation, and, second, our modal decomposition is applied to the stream function formulation with a different set of vertical basis functions allowing for faster convergence of the numerical solution.

To our knowledge, although several approaches to study IT generation are available, some questions remain that we seek to address in this paper. First, how does the nature of the barotropic forcing influence the nature of the  generated ITs? Indeed, the boundary-forcing approach is associated with a barotropic flow that is uniform with depth, whereas the body forcing approach allows for a topographically adjusted barotropic flow which satisfies the correct bottom boundary-condition. The comparison of the two approaches has not been clearly addressed in the literature.
Second, how valid are the analytical predictions of the conversion rate based on the WTA or the knife-edge approximation? This type of predictions are often used in global ocean models \citep{CLEMENT,lahaye2020b}, or for comparisons with estimates made at large supercritical topographic features when investigating specific locations \citep{Gerkema2004, Qian2010, Rainville}. To our knowledge, the limit of validity of these predictions in terms of topographic features (height, slope) has not been sufficiently quantified \citep{Zhang2017}.
A third question relates to the recent results of \cite{maas_2011} concerning specially constructed non-radiating topographies. They also require more investigations in order to assess if such a peculiar phenomenon is more generally observed for any type of topography.
Here we discuss all these aspects, based on a description of the barotropic forcing that consistently takes into account large bottom variations and a modal description of the baroclinic response.

The paper is organized as follows. In section 2, we begin with a review of the basic governing equations, in which we incorporate an arbitrary volume forcing, and we present an asymptotic solution for the barotropic flow. We then formulate the internal tide generation problem in section 3 and derive the associated energy equation and conversion rates. The coupled-mode reformulation of this problem is presented in section 4. In section 5, we present numerical convergence results of our method and give some examples of calculated flow fields. In section 6, we present our results on the conversion rates for various topographies in comparison with existing approaches. Finally, in section 7, we present our conclusions.

\section{Governing Equations\label{sec:governing}}

\subsection{Posing the problem\label{subsec:posing}}
We consider a two-dimensional Cartesian coordinate system $Oxz$, with the vertical $z$ axis pointing upwards,  describing a horizontally infinite layer of a density stratified fluid bounded from above by a ``rigid lid'' $\{z=0\}$, and from below by the bottom surface $\{z=-h_0+h(x)\}$ with $0<h<h_0$. Specifically, we consider bottom topographies $h(x)$ of characteristic height $\Lambda$ and characteristic horizontal scaling length $L$. We do not limit ourselves to small heights, i.e., $\Lambda$ does not have to be much smaller than $h_0$ or any other characteristic vertical length scale. However, we require the topography to be isolated or asymptotically flat, that is, its slope must vanish at infinity, $\lim_{x\to\pm\infty}[\partial_x h] =0$, and we define  $\lim_{x\to\pm\infty}h = h_\pm$.
The latter requirement allows us to take into account fluid domains with different depths at infinity, e.g., continental slopes and shelves.

In static equilibrium the fluid velocity is zero and the background density stratification is $\rho_0(z)$, weakly departing from a constant reference density $\overline\rho_0$ such that the Boussinesq approximation applies. Furthermore, the static density profile decreases linearly with height, such that the Brunt-V\"ais\"al\"a frequency $N = \sqrt{(-g/\overline\rho_0)\text d\rho_0/\text dz}$ is constant. The hydrostatic pressure $p_0(z)$ is then defined by its vertical gradient $p_{0,z} \equiv \text d p_0/\text dz= -\rho_0 g$. The hydrodynamic problem
is posed on the $f$-plane, with $f$ the Coriolis parameter. Our aim is to find perturbations of this state driven by the interaction of an ambient tidal flow and the bottom topography.
This barotropic tidal forcing is characterized by an angular frequency $\omega \in (f, N)$ for $N>f$ and a constant volume flux $Q$.

\begin{figure}[h]
\captionsetup{width=.9\linewidth}
    \centering
    \includegraphics[width=5.5in]{Figures/sketch.eps}
    \caption{Sketch of the set-up and summary of several parameters, used in this article. We use a supercritical topography ($\varepsilon>1$) for purposes of illustration, although our model can also be applied to subcritical ones.}
    \label{fig:sketch}
\end{figure}

Under the Boussinesq approximation, our setup is characterized by seven dimensional parameters ($h_0$, $L$, $\Lambda$, $f$, $N$, $\omega$ and $Q$, and we defer discussion of $h_\pm$ to later), which we summarize in figure \ref{fig:sketch}, measured in combinations of meters and seconds.
A complete dynamical description of our system therefore requires five non-dimensional numbers.
\cite{Garrett_Kunze_2007} highlight a procedure to choose them, from which we deviate somewhat.
First, we introduce a ``funnelling'' ratio,
\begin{equation}
    \label{eq:delta}
    \delta = \frac{\max\{h\}}{h_0} \sim \frac\Lambda{h_0},
\end{equation}
that measures the reduction in cross-section of the flow.
The second and third parameters are the non-dimensional frequency $\omega/f$ and the internal wave characteristic slope,
\begin{equation}
    \mu^{-1} = \sqrt{\frac{\omega^2 - f^2}{N^2 - \omega^2}} = \tan\theta,
    \label{eq:mu}
\end{equation}
$\theta$ being the angle of the free internal wave group velocity with respect to the horizontal plane,
which derives from the hyperbolic operator that describes the linear internal-wave dynamics (see section \ref{sec:governing}\ref{subsec:StreamFctFormulation}). Our fourth non-dimensional parameter is the relative steepness,
\begin{equation}
    \label{eq:epsilon}
    \varepsilon = \mu\max\{|\partial_x h|\} \sim \frac{\mu\Lambda}L,
\end{equation}
and we shall use the first, more accurate version of its definition above (see Section \ref{sec:Validation}).
This parameter's first purpose is to measure the criticality of the topography, with $\varepsilon < 1$ ($>1$) corresponding to the subcritical (supercritical) regime. Fifth and finally, we have the tidal excursion,
\begin{equation}
    \tau = \frac{Q}{(h_0 - \Lambda)\omega L},
\end{equation}
which compares the typical displacement amplitude of a water parcel above the topography, $Q/[(h_0-\Lambda)\omega]$, with the horizontal scale $L$.
Should this parameter be finite, the curvature of the particle trajectories at the bottom would generate internal waves of frequencies, other that $\omega$ \citep{Bell75}.
To ensure monochromatic disturbances, we therefore assume $\tau \ll 1$.
In the ocean, for the lunar semi-diurnal tide M$_2$, the tidal excursion over flat bottom $Q/(\omega h_0)$ is $O(100\text{ m})$ \citep{Bell75}, and therefore, even a moderate topographic width $L=O(10\text{ km})$ would satisfy this so-called ``acoustic limit''.

Of the five parameters defined above, $\omega/f= O(1)$ corresponds to a tidal component at mid-latitudes, and we will hold $\mu^{-1}$ fixed to a small value (for purposes of illustration, $f=10^{-4}$ s$^{-1}$ is the value around latitude 45$^\circ$N, $\omega/f = 1.4$ for the M$_2$ component and $N = 1.5\times 10^{-3}$ s$^{-1}$, meaning that $\mu\approx 15.2$). Starting in Section \ref{sec:Validation}, the parameters $\varepsilon$ and $\delta$ will be the parameters we vary primarily.
Finally, for the seamount (respectively shelf) cases,
we adopt a standard fixed value for $Q = 120$ m$^2$\,s$^{-1}$ (resp. $Q = 134$ m$^2$\,s$^{-1}$), corresponding e.g.\ to a barotropic velocity amplitude at $x\to-\infty$ of $U_0 =Q/h_0 = 4$ cm\,s$^{-1}$ (resp. $U_0 = 4.46$ cm\,s$^{-1}$) and a depth $h_0 = 3$ km.
With the variations of $\varepsilon$ and $\delta$, the value of $\tau$ will therefore vary across calculations while remaining small.
From these five parameters, we can derive other non-dimensional parameters that will prove useful in characterizing the barotropic response.
These are
\begin{equation}
    \mu_0^{-1} = \sqrt{1 - \frac{f^2}{\omega^2}}\qquad \text{and}\qquad \sigma = \frac{h_0^2}{L^2} \sim \left(\frac{\varepsilon}{\mu\delta}\right)^2,
    \label{eq:mu0}
\end{equation}
the former being a frequency parameter appearing in the barotropic elliptical operator, and the latter representing the aspect ratio to be used for the perturbative computation of the barotropic response in subsection \ref{sec:governing}\ref{subsec:barotropic} and Appendix A.

With the parameters defined above, we can discuss the linearity of our equations.
A first way to define it is to estimate the susceptibility of the radiated internal waves to undergo instabilities, which, in the $\tau\ll 1$, $\varepsilon \ll 1$ regime, is small when  $\varepsilon\tau\ll1$ \citep[e.g.,][]{Balmforth_Ierley_Young_2002, Garrett_Kunze_2007}.
In other words, the $\tau\ll 1$, $\varepsilon \ll 1$ regime is linear by construction.
However, no matter the value of $\tau$, linearity breaks down as  $\varepsilon$ increases, implying that this parameter is a better measure of linearity, or absence thereof \citep{Garrett_Kunze_2007, buhler_muller_2007, Grisouard2012}.
In this article, we adopt the common approach of letting $\varepsilon$ be significantly supercritical in some cases, while always solving a linear set of equations \citep{Petrelis_et_al_2006,griffiths2007,Balmforth_Neil_Peacock_2009,echeverri_peacock_2010,Mathur2016}. Indeed, we are primarily interested in predicting conversion from a given large-scale barotropic forcing to a topography-scale response.
The subsequent non-linear evolution of this response, which could take the form of different instabilities \citep[see][for a review]{2018_ARFM_DauxoisJOV}, could be addressed by separate parameterization procedures such as that of \citet{2009_JPO_MullerBuehler}, but would represent a layer on top of the radiation scheme we present in this article.

\subsection{Equations of motion and stream function formulation}\label{subsec:StreamFctFormulation}

Based on the discussion above, we neglect nonlinear effects, as well as diffusion of momentum and buoyancy. The flow may be approximated by the linearized, inviscid Boussinesq equations. Introducing the buoyancy $b= -g(\rho/\overline{\rho}_0 - 1)$, where $\rho \equiv \rho(x,z,t)$ is the \emph{total} density field, the governing equations are
\begin{subequations}\label{eq:Euler_eqs}
\begin{align}
    u_t - f v			& =  -p_x,\quad v_t + f u			 =  0, \label{eq:Eulx} \\
	w_t				    & =  -p_z + b, \label{eq:Eulz} \\
    b_t + N^2 w		    & =  0, \label{eq:buoy}\\
    u_x  + w_z	        & =  0, \label{eq:incompr}
	\end{align}
\end{subequations}
where $(u,v,w)\equiv(u(x,z,t),v(x,z,t),w(x,z,t))$ are the velocity components and $p\equiv p(x,z,t)$ the pressure divided by $\overline{\rho}_0$, and subscripts denote partial derivatives. Additionally, we invoke the non-penetration boundary condition on $z=-h_0+h$ and the rigid-lid approximation on $z=0$ by requiring that
\refstepcounter{equation}
$$
-h_x u(x,-h_0+h) + w(x,-h_0+h)  = 0,\quad w(x,0)  = 0.
\eqno{(\theequation{\text{a},\text{b}})}\label{eq:BCs}
$$
Eq.\ \eqref{eq:incompr} allows to introduce a stream function $\psi(x,z,t)$ such that
\begin{align}
u = -\psi_z,\quad w = \psi_x,
\end{align}
and elimination of $v$, $p$ and $b$ from \eqref{eq:Eulx}--\eqref{eq:buoy} results in a single partial differential equation on $\psi$,
\begin{align}\label{eq:psiPDE}
\psi_{xxtt} + \psi_{zztt} + f^2\psi_{z z} + N^2\psi_{x x}=0.
\end{align}
The rigid-lid condition (\ref{eq:BCs}b) becomes $\psi_x(x,0,t)=(\psi(x,0,t))_x=0$, which implies that $\psi(x,0,t)$ depends solely on time, i.e., $\psi(x,0,t)=c_1(t)$.  Similarly, from (\ref{eq:BCs}a) we have
\begin{align*}
h_x\psi_z(x,-h_0+h,t)+\psi_x(x,-h_0+h,t)=(\psi(x,-h_0+h,t))_x=0
\end{align*}
therefore $\psi(x,-h,t) = c_2(t)$. In order to prescribe $c_1$ and $c_2$, we first note that we may take $c_1(t)=0$ since we can always redefine $\psi$ as $\psi+c_1$. Next, we observe that the volume flux at a vertical section is
\begin{align*}
\int_{-h_0+h}^0u dz = -\int_{-h_0+h}^0\psi_z dz = c_2(t)-c_1(t) = c_2(t).
\end{align*}
In our case, $\psi$ should account for an oscillating flow, hence we choose $c_2(t) = Q \exp(-i\omega t)$ where $Q$ is a constant magnitude and $\omega$ the frequency of the oscillation. Eq.\ \eqref{eq:psiPDE} is thus supplemented by the boundary conditions $\psi(x,0,t) =0$ and $\psi(x,-h_0+h,t) = Q\exp(-i\omega t)$. Looking for time harmonic solutions in the form
\begin{align}
\psi(x,z,t) = \Re\{\phi(x,z)\exp(-i\omega t)\},
\end{align}
where $\phi$ is a complex time independent amplitude and $\Re$ stands for the real part, we obtain the following boundary value problem (BVP)
\begin{align}\label{eq:BVP1_freq}
 \Lin_{\mu} \phi & = 0,\quad
\phi(x,0) =0,\quad
\phi(x,-h_0+h(x)) = Q,
\end{align}
where
\begin{align}
\Lin_{\mu} := \partial_x^2 - \mu^{-2}\partial_z^2,
\end{align}
with $\mu$ defined in \eqref{eq:mu}. The BVP  \eqref{eq:BVP1_freq} can be equivalently written as a BVP with a homogeneous boundary condition on $z=-h_0+h$ by introducing another unknown function $\phi^{\#}(x,z)$ and a known function $\chi(x,z)$, yet to be determined, satisfying the conditions
\begin{align}\label{eq:def_diese}
\begin{aligned}\phi^{\#}= \phi-\chi(x,z)\phi(x,-h_0+h)=\phi-\chi(x,z) Q, \\\quad \chi(x,0)=0,\quad\text{and}\quad \chi(x,-h_0+h) =1.\end{aligned}
\end{align}
The advantage of solving for $\phi^\#$ rather than for $\phi$ is that we are now solving for a volume-forced problem with homogeneous Dirichlet boundary conditions, i.e., $\phi^{\#}$ satisfies
\begin{align}\label{eq:BVP_wforcing}
 \Lin_{\mu} \phi^{\#} &= - Q\Lin_\mu\chi,\quad
\phi^{\#}(x,0) =0,\quad
\phi^{\#}(x,-h_0+h) = 0,
\end{align}
and a straightforward representation of $\phi^\#$ in terms of basis functions (eigenfunctions) that vanish on $z=0, -h_0+h$ can be invoked (Section \ref{sec:CMS}). In order to eventually find the baroclinic response governed by the BVP \eqref{eq:BVP_wforcing}, the function $\chi$ must be determined so that it accounts for an appropriate barotropic flow. Such a flow is considered in the following subsection.

\subsection{Barotropic flow}\label{subsec:barotropic}

We now compute the background barotropic flow that is established in the fluid domain as a response to the external tidal forcing. Later, in section \ref{sec:ITgeneration}\ref{subsec:BCsolution}, we will use this barotropic flow as a volume forcing for the baroclinic component.
The barotropic flow is solution to the system
\begin{subequations}\label{eq:Euler_Barotrop}
	\begin{align}
	U_t - f V			& =  -P_x,\quad V_t + f U=  0, \label{eq:Eulx_Barotrop} \\
	W_t				    & =  -P_z, \label{eq:Eulz_Barotrop}\\
	B_t + N^2 W		& =  0, \label{eq:Buoy}\\
	U_x  + W_z	& =  0,\label{eq:incompr_Barotrop}
    \end{align}
\end{subequations}
where $(U,V,W)$, $P$ and $B$ denote the velocity components, the scaled pressure, and the buoyancy perturbations induced by the oscillating barotropic flow, respectively (see e.g.\ \citep[Sec.\ 8]{Baines1973} or \citep{Garett_Gerkemma_2007}), together with the boundary conditions
\refstepcounter{equation}
$$
-h_x U(x,-h_0+h)+W(x,-h_0+h)  = 0,\quad W(x,0)  = 0.
\eqno{(\theequation{\text{a},\text{b}})}\label{eq:BCs_Barotrop}
$$
Introducing a barotropic stream function $\Psi$ such that $U = -\Psi_z$ and $ W = \Psi_x$, we obtain from the system\ \eqref{eq:Euler_Barotrop} a single equation on $\Psi$,
\begin{align*}
    \Psi_{xxtt}+\Psi_{zztt}+f^2\Psi_{zz} &= 0.
\end{align*}
In terms of $\Psi$, the boundary conditions \eqref{eq:BCs_Barotrop} may be written as $\Psi(x,0,t)=0$ and $\Psi(x,-h_0+h,t)=Q\exp(-i\omega t)$. Assuming that $\Psi = \Re\{\Phi e^{-i \omega t}\}$ for some time independent function $\Phi$, we obtain the following BVP for the barotropic flow
\begin{align}\label{eq:BarotropicBVP}
\Phi_{xx}+ \mu_0^{-2}\Phi_{zz}=0, \quad \Phi(x,0)=0,\quad \Phi(x,-h_0+h)=Q,
\end{align}
with $\mu_0$ defined in \eqref{eq:mu0}. Note that the partial differential equation in \eqref{eq:BarotropicBVP} is elliptic (cf. Eq.\ \eqref{eq:BVP1_freq}). The above problem is solved in a perturbative way in terms of the ratio of the water depth and the horizontal scale of the topography, $\sigma=h_0^2/L^2$. This is done in detail in Appendix A. The net result is that the barotropic flow quantities can be approximated by an expression of the form
\begin{align}
    \Xi = \Xi^{(0)} + \sum\limits_{i=1}^{K}\Xi^{(i)},
\end{align}
for some order $K\geq 1$, where $\Xi$ is a placeholder for any of  $\Phi$, $\Psi$, $U$, $V$, $W$, $B$ and $P$. In particular, $\{\Phi^{(i)}\}_{i=0}^{K}$ are determined recursively (see Eq.\ \eqref{eq:Phi_orders_nondim}); the first two orders are
\begin{align}\label{eq:Phi_orders_dim}
\Phi^{(0)}=Q\frac{z}{-h_0+h},\quad
\Phi^{(1)}=-Q\mu_0^2\left(\frac{1}{-h_0+h}\right)_{xx}\frac{1}{6}\left(z^3-(-h_0+h)^2z\right)
\end{align}
and for $i\geq 0$ we have
\refstepcounter{equation}
$$
U^{(i)} = -\Psi^{(i)}_z = - \Phi^{(i)}_z\cos\omega t,\quad W^{(i)} = \Psi^{(i)}_x= \Phi^{(i)}_x\cos\omega t
\eqno{(\theequation{\text{a},\text{b}})}\label{eq:UW}
$$
\refstepcounter{equation}
$$
V^{(i)} = \frac{f}{\omega}\Phi^{(i)}_z\sin\omega t,\quad B^{(i)} = -\frac{N^2}{\omega} \Phi^{(i)}_x\sin\omega t.
\eqno{(\theequation{\text{a},\text{b}})}\label{eq:VBUO}
$$
As shown in Appendix A, the above solution provides a decomposition of the barotropic flow into a hydrostatic part $\Xi^{(0)}$ and a residual non-hydrostatic part, denoted by $\Xi^{\text r} \equiv \sum\limits_{i=1}^{K}\Xi^{(i)}$. Note that the hydrostatic part depends on the topography $h(x)$ and that the non-hydrostatic part vanishes when $h_x=0$.

\section{The internal tide generation problem\label{sec:ITgeneration}}

\subsection{Formulation\label{subsec:BCsolution}}

Owing to the solution of the barotropic flow developed in the previous section, we proceed by making the internal-tide generation problem \eqref{eq:BVP_wforcing} explicit. Choosing
 $\chi(x,z) Q = \Phi^{(0)} + \Phi^{\text{r}}$, the BVP \eqref{eq:BVP_wforcing} takes the form

\refstepcounter{equation}
$$\label{eq:BVP_freq0NH}
\Lin_{\mu} \phi^{\#}  = -\Lin_{\mu}\left(\Phi^{(0)}+\Phi^{\text{r}}\right),\quad
\phi^{\#}(x,0)=0\quad
\phi^{\#}(x,-h_0+h)= 0,\eqno{(\theequation{\text{a},\text{b},\text{c}})}
$$
This shows that the barotropic flow, $\Phi^{(0)}+\Phi^{\text{r}}$, forces the baroclinic response $\phi^{\#}$. Alternatively, exploiting the linearity of $\Lin_\mu$, the BVP \eqref{eq:BVP_freq0NH} can be equivalently stated in the form
\refstepcounter{equation}
$$\label{eq:BVP_freq0H}
\Lin_{\mu} \phi^{\dagger}  = -\Lin_{\mu}\Phi^{(0)},\quad
\phi^{\dagger}(x,0)=\phi^{\#}(x,0)=0,\quad
\phi^{\dagger}(x,-h_0+h)=\phi^{\#}(x,-h_0+h)= 0,\eqno{(\theequation{\text{a},\text{b},\text{c}})}
$$
with $\phi^{\dagger}=\phi^{\#} +\Phi^{\text{r}}$, implying that the hydrostatic part of the barotropic flow, $\Phi^{(0)}$, forces a baroclinic response as well as a non-hydrostatic barotropic one \citep{Garett_Gerkemma_2007}. Here, we shall proceed with the latter formulation for two reasons. First, the non-hydrostatic barotropic flow does not propagate and therefore does not influence the conversion calculations, making $\phi^\#$ and $\phi^\dagger$ equally useful for this purpose. Second, it is easier to implement numerically because the right-hand side of Eq. (\ref{eq:BVP_freq0H}a), $\Lin_{\mu}\Phi^{(0)}$, contains second order derivatives of the topography $h$, whereas  the right-hand side of Eq. (\ref{eq:BVP_freq0NH}a), $\Lin_{\mu}\Phi^{\text{r}}$, contains derivatives of $h$ of order $2(K+1)$ with $K\geq 1$.

The above set of equations must be supplemented with  radiation conditions for $|x|\rightarrow\infty$ ensuring that the waves generated in the interior of the domain propagate outward at infinity as plane waves. Plane internal waves can be expressed through separation of variables of $\Lin_{\mu} \phi^{\dagger}=0$ on a domain of constant depth $h_{\pm}$ as
\begin{align}\label{eq:radiation_conditions}
\phi^{\dagger}(x,z) &= \sum\limits_{n=1}^{\infty}c^\pm_n\exp(\pm i k^{\pm}_n x)\sin(\kappa_n^{\pm} z),\quad\text{as}\quad x\rightarrow\pm\infty,
\end{align}
where $c^{\pm}_n\in\mathbb{C}$ are constants,
\begin{align}\label{eq:wavenumbers_farfield}
\kappa^{\pm}_n= \frac{n\pi}{-h_0+h_{\pm}},\quad k^{\pm}_n = \frac{\kappa^{\pm}_n}{\mu},\quad n=1,2,\dots,
\end{align}
and $h\rightarrow h_{\pm}$ as $x\rightarrow\pm\infty$.

Assuming that the solution $\phi^{\dagger}$ has been found, the purely baroclinic response is $\phi^{\#} = \phi^{\dagger}-\Phi^{\text r}$.
In Section \ref{sec:Validation}, we calculate the flow fields by using $\Phi^{\text{r}}=\Phi^{(1)}$, where $\Phi^{(1)}$ is given in \eqref{eq:Phi_orders_dim}. This choice is justified since $\sigma<0.01$ in the presented solutions and keeping a higher order in $\Phi^{\text{r}}$ does not result in a significant correction of the baroclinic response.

It will later prove useful to write down equations for the rest of the response fields. Let us first note that, in view of the first equation in \eqref{eq:def_diese} and the analysis in subsection \ref{sec:governing}\ref{subsec:barotropic}, the total field $\phi$ is written as $\phi = \Phi^{(0)}+\Phi^{\text{r}}+\phi^{\#} = \Phi^{(0)}+\phi^{\dagger}$. Analogously, we introduce the decomposition
\begin{align}\label{eq:fdecomposition}
\begin{aligned}
    \xi &= \Xi^{(0)} + \Xi^{\text r} + \xi^\# = \Xi^{(0)} + \xi^\dagger,
\end{aligned}
\end{align}
where, like for the barotropic response, $\xi$ is a placeholder for any of $\phi$, $u$, $v$, $w$, $b$ or $p$.
Plugging the second equality of \eqref{eq:fdecomposition} into Eqs.\ \eqref{eq:Euler_eqs}--\eqref{eq:BCs} and taking into account Eqs.\ \eqref{eq:Euler_Barotrop_asympt}--\eqref{eq:BCs_Barotrop_asympt} with $i=0$, we obtain
\begin{subequations}\label{eq:Eul_Baroclin}
\begin{align}
    u^{\dagger}_t - f v^{\dagger}			 & =  -p^{\dagger}_x, \quad
	v^{\dagger}_t + f u^{\dagger}			 =  0, \label{eq:Eulxy_Baroclin} \\
	w^{\dagger}_t				    & =  -p^{\dagger}_z + b^{\dagger} + B^{(0)}  - W_t^{(0)}, \label{eq:Eulz_Baroclin}
	\end{align}
\end{subequations}
\begingroup\abovedisplayskip=0pt
\begin{align}	b^{\dagger}_t + N^2 w^{\dagger}		& = 0, \label{eq:buoy_Baroclin}
\end{align}
\endgroup
\begingroup\abovedisplayskip=-5pt
\begin{align}
 u^{\dagger}_x  + w^{\dagger}_z	& =  0,\label{eq:incompr_Baroclin}
\end{align}
\begin{align}
-h_x u^{\dagger}(x,-h_0+h)+w^{\dagger}(x,-h_0+h)  = 0,\quad w^{\dagger}(x,0)  = 0.
\label{eq:BCs_Baroclin}
\end{align}
\endgroup
Linearity of the problem allows us to assume that all daggered variables in Eqs.~\eqref{eq:Eul_Baroclin}--\eqref{eq:BCs_Baroclin} have the same time-periodicity as $\phi^\dagger$, e.g., $u^\dagger = \Re\left\{\phi^\dagger_z e^{-i\omega t}\right\}$.
We can therefore deduce from the previous statement and  Eqs.~\eqref{eq:Eul_Baroclin}--\eqref{eq:incompr_Baroclin} the standard polarization for internal waves, i.e.,
\begin{equation}\label{eq:responnse}
    \left(u^\star,\ v^\star,\ w^\star,\ b^\star\right) = \Re\left\{\left(-\phi^\star_z,\ \frac{if}\omega\phi_z^\star,\ \phi_x^\star,\ -\frac{iN^2}{\omega}\phi^\star_x \right)e^{-i\omega t}\right\},
\end{equation}
where the superscript $\star$ stands for either $\dagger$ or $\#$, i.e., the above polarization relations apply to both the combined baroclinic and non-hydrostatic response $(\dagger)$, and to the baroclinic response alone $(\#)$. For easy reference, we summarize the notations for the flow fields that have been defined in table~\ref{tab:fields_notations}.
\begin{table}
    \centering
    \begin{tabular}{c|c|c}
        Notation & Comments & Related definition(s) \\\hline
        $\xi$ & Solution to the full problem \eqref{eq:Euler_eqs}--\eqref{eq:BCs} & $\xi = \Xi + \xi^\# = \Xi^{(0)} + \xi^\dagger$ \\
        $\Xi$ & Barotropic component, solution to \eqref{eq:Euler_Barotrop}--\eqref{eq:BCs_Barotrop} & \\
        $\xi^{\#}$ & Baroclinic component & See eqns.~\eqref{eq:def_diese} and \eqref{eq:BVP_freq0NH} \\
        $\Xi^{(n)}$, $n\in \mathbb N$ & Asymptotic components of $\Xi$, ordered by $\sigma$ & $\Xi = \sum_{n=0}^\infty\sigma^n\Xi^{(n)}$ \\
        $\Xi^{(0)}$ & Leading-order barotropic response, hydrostatic & \\
        $\Xi^\text{r}$ & Higher-order barotropic response, non-hydrostatic & $\Xi^\text{r} = \Xi - \Xi^{(0)} = \sum_{n=1}^\infty\sigma^n\Xi^{(n)}$ \\
        $\xi^\dagger$ & Flow response to $\Xi^{(0)}$ as a forcing & $\xi^\dagger = \xi^\# - \Xi^\text{r}$; see eqn.~\eqref{eq:BVP_freq0H} \\
    \end{tabular}
    \caption{Various notations for the fields, we solve for.
    The symbol $\xi$ stands for $u$, $v$, $w$, $b$, $p$, $\psi$ or $\phi$, and $\Xi$ for their capitalized versions.}
    \label{tab:fields_notations}
\end{table}

\subsection{Energy equation and conversion rate}
Finally, we derive the energy equation associated with the internal-tide generation problem formulated in the previous subsection, to  be used later for validation purposes (Section \,\ref{sec:Validation}\ref{sec_convergence}) and to define conversion rates and energy fluxes. Taking the dot product of Eqs.\ \eqref{eq:Eul_Baroclin}--\eqref{eq:buoy_Baroclin} with the state vector $\left(\mathbf{u}^{\dagger},b^{\dagger}/N^2\right) \equiv (u^{\dagger}, v^{\dagger}, w^{\dagger},b^{\dagger}/N^2)$ and using Eq.\ \eqref{eq:incompr_Baroclin}, we obtain
\begin{align}
    \mathcal{E}^\dagger_t + \nabla\cdot\left(p^\dagger\mathbf{u}^\dagger\right) = S,\label{eq:en_density}
\end{align}
where $\nabla = \left(\partial_x, \partial_z \right)$, $\mathcal{E}^\dagger(x,z,t)$ is the energy density associated with the baroclinic and non-hydrostatic barotropic flows,
and $S(x,z,t)$ is a source term due to the conversion from the hydrostatic-barotropic flow, i.e.,
\begin{align*}
    \mathcal{E}^\dagger = \frac{1}{2}(\mathbf{u}^{\dagger})^2 +\frac{1}{2}\frac{(b^{\dagger})^2}{N^2
    \quad\text{and}\quad S=  \left(B^{(0)}-W^{(0)}_t\right)w^\dagger.
\end{align*}
Integrating Eq.\ \eqref{eq:en_density} over the domain $\Omega = [x_L,x_R]\times [-h_0+h(x),0]$, we introduce the instantaneous energy
\begin{align*}
E^\dagger(t) = \int_\Omega \mathcal{E}^\dagger d\Omega = \int_{x_L}^{x_R}\int_{-h_0+h}^0 \mathcal{E}^\dagger dx dz,
\end{align*}
and invoking the divergence theorem for the flux term, and using the boundary conditions in \eqref{eq:BCs_Baroclin}, we obtain
\begin{align}
    E^\dagger_t - \int_{-h_0+h(x_L)}^0 p^{\dagger} u^{\dagger} dz + \int_{-h_0+h(x_R)}^0 p^{\dagger} u^{\dagger} dz
    =  \int_{\Omega}\left(B^{(0)}  - W_t^{(0)}\right)w^{\dagger} d\Omega,
\end{align}
which expresses the energy balance over $\Omega$. Applying the phase average operator $\left\langle \,\cdot\,\right\rangle  = 1/T \int_0^T \cdot\, dt$, with $T = 2\pi/\omega$, on the above equation and taking into account the periodicity of $E^{\dagger}(t)$, yields
\begin{align}\label{eq:energy_eq_avg}
- \int_{-h_0+h(x_L)}^0 \left\langle p^{\dagger} u^{\dagger}\right\rangle dz  + \int_{-h_0+h(x_R)}^0 \left\langle p^{\dagger} u^{\dagger}\right\rangle dz
    =  \int_{\Omega} \left\langle \left(B^{(0)}  - W_t^{(0)}\right)w^{\dagger}\right\rangle  d\Omega.
\end{align}
The two terms in the left-hand side of the above equation represent the mean energy flux crossing the vertical boundaries of $\Omega$, at $x_L$ and $x_R$. At the far-field ($x_L\rightarrow\infty$ and $x_R\rightarrow\infty$), these terms represent the rate at which energy is radiated away from the topography and are called conversion rates. They are defined here by
\begin{align}\label{eq:conv_rate}
    C_{\pm} = \lim_{x\rightarrow \pm\infty}\left[\pm \int^0_{-h_0+h(x)} \left\langle p^{\dagger} u^{\dagger}\right\rangle dz\right].
\end{align}
In order to calculate $ C_{\pm}$, we first note that the integral in \eqref{eq:conv_rate} can be expressed in terms of $p^{\dagger}_z$ by using integration by parts and the fact that $\psi^{\dagger}$ vanishes on $z=0,\,\,-h_0+h(x)$:
\begin{align*}
\int_{-h_0+h(x)}^0 \left\langle  p^{\dagger} u^{\dagger}\right\rangle  dz =   -\int_{-h_0+h(x)}^0 \left\langle p^{\dagger} \psi_z^{\dagger}\right\rangle  dz =  \int_{-h_0+h(x)}^0 \left\langle p_z^{\dagger} \psi^{\dagger}\right\rangle  dz.
\end{align*}
Using the expression for $p_z^{\dagger}$ obtained from eqs.\ \eqref{eq:responnse} and
\eqref{eq:Eulz_Baroclin}, we obtain
\begin{align}
\left\langle p_z^{\dagger}\psi^{\dagger}\right\rangle  = \frac{\omega^2-N^2}{\omega} \left<\Re\{i\phi^{\dagger}_x e^{-i\omega t}\}\Re\{\phi^{\dagger} e^{-i\omega t}\} + \Phi_x^{(0)}\sin\omega t\, \Re\{\phi^{\dagger}e^{-i\omega t}\}\right>.
\end{align}
Writing $e^{-i\omega t} = \cos\omega t - i \sin\omega t $ and noting that $\left\langle\cos\omega t\sin\omega t\right\rangle = 0$, $\left\langle\cos^2\omega t\right\rangle = \left\langle\sin^2\omega t\right\rangle = 1/2$ and $\lim_{|x|\to\infty}[\Phi_x^{(0)}]= 0$ yields
\begin{align}\label{eq:Conv}
C_{\pm}  = \frac{N^2-\omega^2}{2\omega}\int_{-h_0+h_{\pm}}^0 \Im\left(\phi^{\dagger}\overline{\phi_x^{\dagger}}\right)  dz,
\end{align}
where the overline denotes the complex conjugate, $\Im$ denotes the imaginary part, and where we recall that $\lim_{|x|\to\infty}[h]= h_{\pm}$.

\section{Modal decomposition}\label{sec:CMS}
\subsection{The Coupled-Mode System}\label{subsec:CMS}
In order to solve the problem \eqref{eq:BVP_freq0H}--\eqref{eq:radiation_conditions},  we follow a modal decomposition approach. The unknown function $\phi^{\dagger}(x,z)$ is represented as a series of the form
\begin{align}
\phi^{\dagger}(x,z) = \sum\limits_{n=1}^{\infty}\phi_n(x)Z_n(z;x),\label{eq:RL_representation}
\end{align}
where $\{\phi_n(x)\}_{n=0}^{\infty}$ are the unknown complex modal amplitudes to be determined and $\{Z_n(z;x)\}_{n=0}^{\infty}$ are prescribed vertical basis functions with a parametric dependence on $x$. In order for the expansion in eq.\ \eqref{eq:RL_representation} to be exact, the set $\{Z_n(z;x)\}_{n=0}^{\infty}$ must be complete. Such a set is obtained as the set of local eigenfunctions associated with a  Sturm-Liouville problem parametrized by $x$, also called ``reference waveguide" \citep{Brekhovskikh}. In our case, the reference waveguide reads \begin{subequations} \label{eq:SLproblem}
    \begin{align}
        Z_{n,zz} + \tilde{\kappa}_n^2Z_n &= 0,\quad- h_0+h(x)\leq z\leq 0,\\
        Z_n &= 0,\quad z=0,\\
        Z_n &= 0,\quad z= - h_0+h(x).
    \end{align}
\end{subequations}
The corresponding local eigenfunctions and eigenvalues are given by
\refstepcounter{equation}
$$
Z_n(z;x) =  \sin\left(\tilde{\kappa}_n z\right), \quad \tilde{\kappa}_n= \tilde{\kappa}_n(x) = \frac{n\pi}{-h_0+h(x)}, \eqno{(\theequation{\text{a},\text{b}})}\label{eq:RL_vertical_function}
$$
and satisfy the orthogonality relation $\int_{-h_0+h}^{0}Z_n Z_m dz = (h_0-h)\delta_{nm}/2$, where $\delta_{nm}$ is the Kronecker delta. It follows that each $\phi_n$ is defined
by
\begin{align}\label{eq:def_mod}
    \phi_n(x) = \frac{2}{h_0-h(x)}\int_{-h_0+h(x)}^0\phi^{\dagger}(x,z) Z_n(z;x) dz.
\end{align}
Note that this eigenfunction expansion is consistent with the boundary conditions in \eqref{eq:BVP_freq0H}. Furthermore, provided that $\phi^{\dagger}$ is sufficiently smooth, a repeated integration by parts in \eqref{eq:def_mod} yields
 \begin{align}\label{eq:estimate}
\begin{aligned}
    \phi_n(x)=\frac{2}{(h_0-h)\tilde{\kappa}_n^3}\left( \left[\phi^{\dagger}_{zz}Y_n\right]_{-h_0+h}^0 - \int_{-h_0+h}^0 \phi^{\dagger}_{z z z}Y_n dz\right) ,
    \end{aligned}
 \end{align}
 with $Y_n=\cos(\tilde{\kappa}_n z)$.
This shows that $\|\phi_n\|_{\infty}:=\max |\phi_n|=O(n^{-3})$ and that the eigenfunction expansion
converges uniformly.
Similar estimates are obtained for the decay of $\|\phi_{n,x}\|_{\infty}$ and $\|\phi_{n,xx}\|_{\infty}$ by adapting the procedure developed in \citet[Section 4]{Ath_Pap_2017}. These estimates suffice to establish the term-wise differentiability of \eqref{eq:RL_representation} and justify its use in conjunction with Eq.\ \eqref{eq:BVP_freq0H}.

We proceed by projecting Eqs.\ \eqref{eq:BVP_freq0H}--\eqref{eq:radiation_conditions} onto \eqref{eq:RL_vertical_function}. Substituting \eqref{eq:RL_representation} into the differential equation in Eq.\ \eqref{eq:BVP_freq0H}, multiplying with $Z_n$ and integrating over the interval $[-h_0+h(x),0]$, we find that the $\phi_n(x)$ are solutions to the Coupled-Mode System (CMS)
\begin{align}\label{eq:CMS}
\sum\limits_{n=1}^{\infty}A_{mn}(h)\phi_{n,xx}+B_{mn}(h)\phi_{n,x}+C_{mn}(h)\phi_n=F_m(h),\quad m\geq 1,
\end{align}
where the matrix coefficients $A_{mn}=A_{mn}(h)$, $B_{mn}=B_{mn}(h)$, $C_{mn}=C_{mn}(h)$ and the sequence $F_m=F_m(h)$ are provided in Appendix B.  The CMS \eqref{eq:CMS} must be supplemented with appropriate radiation conditions on the modal amplitudes $\phi_n$. In view of Eq.\ \eqref{eq:RL_representation}, the radiation conditions \eqref{eq:radiation_conditions} are written
\begin{align}\label{eq:BVP2_RCleft_m}
\sum\limits_{n=1}^{\infty}\phi_n(x)Z_n(x,z) = \sum\limits_{n=1}^{\infty}c_n\exp(\pm i k^{\pm}_n x)\sin(\kappa^{\pm}_n z),\quad\text{as}\quad x\rightarrow \pm\infty.
\end{align}
Multiplying both sides by $Z_q$, $q\in \{1, 2, \dots\}$, integrating over $[-h_0+h,0]$ and taking into account that  $Z_n\rightarrow \sin(\kappa^{\pm}_n z)$ as $ x\rightarrow \pm \infty$, we obtain
\begin{align*}
\forall n\in \{1, 2, \dots\},\quad\phi_n = c_n\exp(\pm i k^{\pm}_n x),\quad\text{as}\quad x\rightarrow\pm\infty.
\end{align*}
and therefore
\begin{align}\label{eq:rad_cond_left}
\partial_x\phi_n \pm i k^{\pm}_n\phi_n= 0,\quad\text{as}\quad x\rightarrow {\pm}\infty.
\end{align}
Thus, the IT generation problem has been exactly reformulated as the CMS  \eqref{eq:CMS} together with the radiation conditions \eqref{eq:rad_cond_left}. A numerical method for its solution is presented in the following subsection.

\subsection{Numerical solution}\label{subsec:numerical_scheme}
In order to solve numerically the infinite CMS Eqs.\ \eqref{eq:CMS}-\eqref{eq:rad_cond_left},
a finite number $M$ is chosen as the truncation order and the infinite domain is replaced by a finite interval $[x_L,x_R]$. In practice, $x_L$, $x_R$ are chosen so that $h_{x}(x_L)$, $h_{x}(x_R)$ are negligible. The truncated CMS reads
\begin{subequations}\label{eq:CMS_truncated}
\begin{align}
\sum\limits_{n=1}^{M}A_{mn}\phi_{n,xx}+B_{mn}\phi_{n,x}+C_{mn}\phi_n&=F_m,\quad x\in[x_L,x_R],\\
\phi_{m,x} +ik_m\phi_m&= 0,\quad\text{at}\quad x=x_\text{L},\\
\phi_{m,x} -ik_m\phi_m&= 0,\quad\text{at}\quad x=x_\text{R}.
\end{align}
\end{subequations}
with $m= 1,\dots,M$. For the numerical solution of the above problem, we discretize the interval $[x_L,x_R]$ using a uniform grid of spacing $\delta x$, $\{x_i$, $i=1,N_X\}$, and compute the grid values $\phi^i_n$ that approximate the solution $\phi_n(x_i)$.  The derivatives in \eqref{eq:CMS_truncated} are approximated  with $4^{th}$-order central finite differences except for the boundary points and their neighboring points where one-sided and asymmetric finite-differences are used in order to preserve the $4^{th}$-order accuracy. The corresponding formulae can be found in \citep[Appendix C]{PAPOUTSELLIS2019103579}. Using these expressions, we obtain a sparse square linear system of dimension $(N_X M)^2$ for the determination of $\left\{\phi^i_n, n\in[1,M], i\in[1,N_X]\right\}$. The system's matrix is comprised by $M^2$ blocks of dimension $N_X^2$ and every block has an almost five-diagonal structure. Its solution is then obtained by means of a $LU$ decomposition.

Once the solution is found, the baroclinic fields  $u^{\#}$, $v^{\#}$, $w^{\#}$ and $b^{\#}$ are constructed using a truncated version of \eqref{eq:responnse}. Similarly, from Eq.\ \eqref{eq:Conv}  the conversion rate is computed by
\begin{align}\label{eq:Conv_modal}
C_{\pm} = \frac{N^2-\omega^2}{2\omega} \frac{h_0-h_{\pm}}{2}\sum\limits_{n=1}^{M} \Im\left[\phi_n\overline{\phi}_{n,x}\right]_{x=x_R,x_L}.
\end{align}
Finally, we also provide here an estimation of the free-surface elevation due to the internal-tide motion given by $\eta = [p^{\#}]_{z=0}/ g$, where $[p^{\#}]_{z=0}$ is the pressure induced on $z=0$ by the baroclinic flow. Details are given in Appendix C.

\section{Cases studied and numerical convergence \label{sec:Validation}}

\subsection{Topographic profiles and numerical values}

In the following section, we calculate ITs for three bathymetries, namely the ``witch of Agnesi", the ``Gaussian" and what we will call the ``bump", which are given by
\begin{align}\label{eq:topos}
    h_{\text{W}}=\frac{\Lambda}{1+\frac{x^2}{L^2}},\quad
    h_{\text{G}}=\Lambda\exp\left(-\frac{x^2}{2L^2}\right),\quad
    h_{\text{B}}=\Lambda\exp\left(1-\frac{1}{1-\frac{x^2}{L^2}}\right)\mathds{1}_{(-L,L)},
\end{align}
with $\mathds{1}_{(-L,L)}=1\,\,\text{in}\,\,(-L,L)$ and $\mathds{1}_{(-L,L)}=0$ otherwise. Note that the ``bump" profile has a compact support in contrast with the other two profiles that tend to zero at infinity; see figure \ref{Fig:seamounts}. We shall also consider the case of a shelf connecting two different constant depths. The shelf profile is the same as in \citet[Section 5]{griffiths2007}, hereafter referred to as GG07:
\begin{align*}\label{eq:shelf_topo}
    h = \left\{\begin{aligned}
        & h_-,\quad x\leq 0\\
        & h_- + (h_+-h_-) \sin^2(\pi x/2L), \quad 0\leq x\leq L\\
        & h_+,\quad x\geq L
    \end{aligned}\right. .
\end{align*}
We recall from subsection \ref{sec:governing}\ref{subsec:posing} that two important parameters will be considered, namely, the criticality of the bottom slope and the relative topography height
\begin{align}
\varepsilon = \frac{\max\{|h_x|\}}{\tan\theta}\quad\text{and}\quad \delta = \frac{\max\{h\}}{h_0}.
\end{align}
We also recall that, for illustration purposes, $(\omega,U_0)$ correspond to a typical $M_2$ tide, i.e., $(\omega,U_0) = (1.4\times10^{-4}\ \text{s}^{-1}, 0.04\,  \text{m\,s}^{-1})$ and that the Brunt-V\"ais\"{a}l\"{a} and Coriolis frequencies are  $ N = 0.0015\ \text{s}^{-1}$ and $f = 10^{-4} \ \text{s}^{-1}$ respectively. For all seamounts, the depth far from the topography is $h_0 = 3000$~m.
\begin{figure}
\captionsetup{width=.9\linewidth}
\centering
  \includegraphics[scale = 0.7]{Figures/Fig_seamounts.eps}
\caption{The three considered seamount profiles.}
\label{Fig:seamounts}
\end{figure}

\subsection{Truncation and numerical convergence}
\label{sec_convergence}
We first examine the rate of decay of the modal amplitudes with respect to mode number obtained by solving the proposed CMS \eqref{eq:CMS_truncated}. For this calculation, the ``bump" profile is considered with $\Lambda=1500$ m ($\delta=0.5$) and  $\varepsilon = 0.9$, $1$ and $1.5$. Results on the uniform norm $\|\cdot\|_{\infty}$ of $\phi_n$, $\phi_{n,x}$ and $\phi_{n,xx}$ are shown in the first row of figure \ref{fig:Error_decay}. For $\varepsilon=0.9$, $\|\phi_n\|_{\infty}$  decays like  $n^{-s}$ with $s\approx 3$ for $n>10$ while $s\approx 3.8$ and $4.0$ for $\phi_{n,x}$ and $\phi_{n,xx}$ in the same range. Similar decay rates are observed for $\varepsilon=1$, though in this case they take effect after $n>30$. For $\varepsilon=1.5$, $\|\phi_n\|_{\infty}$ decays more slowly, reaching a minimum at $n_{\min}=66$. In this case, a linear fit gives $\|\phi_n\|_{\infty}\sim n^{-s}$ with $s\approx 1.5$ for $n\in[30,n_{\min}]$.  The first derivative decays with a smaller rate, $s\approx 0.5$  for $n\in[30,n_{\min}]$, while the second derivative clearly diverges. This change in the decay rate is related with the formation of singularities for $\varepsilon>1$ and is also reported in the modal solution of GG07. In fact, as shown in Section \ref{sec:CMS}\ref{subsec:CMS}, $\|\phi_n\|=O(n^{-3})$ if $\phi^{\dagger}$ is smooth, in accordance with our numerical calculations when $\varepsilon<1$. On the other hand, a decay rate $n^{-s}$, with $s=3/2$ suggests a formation of a singularity in $z$ takes place \citep{Salem1939, Raisbeck1955}.
For all the considered topographies in this paper, the rapid decay rate $n^{-s}$, with $s\approx 3$, of $\|\phi_n\|_\infty$ holds for $\varepsilon\leq 1$ and $\delta>0.1$. We observe that for smaller $\delta$ and sufficiently large $\varepsilon$ ($\varepsilon\gtrapprox 0.2$), the decay rate approaches the one obtained with $\varepsilon>1$, namely, $s\approx 1.5$, showing that the solution becomes singular also in this case. For $\varepsilon>1$, the slow convergence $s\approx 1.5$ holds for all $\delta$.

Next, in order to assess the accuracy of the CMS approach and the numerical scheme, we compare the energy radiated away from the topography with the energy conversion in the entire domain. To do so, we solve Eqs. \eqref{eq:CMS_truncated} for a wide range of the parameters $(M,\delta x)$ and we monitor the decay of the corresponding absolute error of the averaged energy equation Eq.\ \eqref{eq:energy_eq_avg},
\begin{align}\label{eq:energy_eq_avg_error}
\text{Err.} =  \Bigl|- C_-  + C_+ -  \int_{\Omega} \left\langle B^{(0)}w^{\dagger}  - W_t^{(0)}w^{\dagger}\right\rangle  d\Omega\Bigr|,
\end{align}
where $C_{\pm}$ are given by Eq.\ \eqref{eq:Conv}. We use $M = 10$ to $120$ modes and choose to resolve each case with $\delta x = L_M/r$ where $L_M$ is the horizontal wavelength of the $M^{\text{th}}$ free internal wave mode over the summit of the seamount, i.e.\ $L_M = 2\mu(h_0-\Lambda)/M$, and $r= 4,6,8,10$ and $12$. Results for  $h=h_{\text{B}}(x)$ with $h_\text{r}=1500$ m ($\delta=0.5$), and $\varepsilon = 0.9$, $1$ and $1.5$ are shown in figure \ref{fig:Error_decay}. The expected $r^{-4}$ rate is obtained for all $\varepsilon$ verifying the fourth-order accuracy of the present spatial discretisation scheme. For $\varepsilon=0.9$ and $\varepsilon=1$, the error decays rapidly as $M^{-4}$ and $M^{-3.5}$ respectively. For $\varepsilon = 1.5$, the error decays as $M^{-1}$ showcasing the slow convergence of the modal solution when the underlying field becomes singular.

At this point, some general remarks are in order.
For the subsequent analysis, the calculations are performed with $M=50$ and $M=100$ when $s\approx3$ and $s\approx1.5$ respectively, while $\delta x$ is at least six times smaller than $L_M$ in order for the last calculated mode to be sufficiently resolved. After extensive numerical investigation, we have observed that increasing the resolution $(M,r)$ further than the aforementioned values does not significantly change the energy flux nor the first thirty modes, thus these values may be considered a good compromise between efficiency and accuracy.
It must be stressed that our rates of decay are consistent with the $n^{-2}$ and $n^{-1/2}$ rates reported by GG07 for the modal amplitudes of the horizontal velocity; indeed, with our notation, $u_n\sim \phi_n \tilde{\kappa}_n \cos(\tilde{\kappa}_n z)\sim \phi_n n \cos(\tilde{\kappa}_n z)$. Nevertheless, if 2D solutions are sought, it is more efficient to work with the stream function formulation.

\subsection{Examples of flow fields}

In the remainder of this section, we consider converged solutions and the corresponding magnitude of the baroclinic velocity field $|\mathbf{u}^{\#}| = \sqrt{(u^{\#})^2+(v^{\#})^2+(w^{\#})^2}$, calculated by Eq. \eqref{eq:responnse} and normalized by $U_0$, isodensity curves and free-surface perturbation calculated by Eq.\ \eqref{eq:fs_pert}. In figure \ref{Fig:flowmaps_bump}, we plot the results for two characteristic instants at $t=0$ and $t=T/4$, for the bump profile introduced above. We describe below the main features of the obtained solutions. The full evolution can be found in the supplemental material.

For $\varepsilon=0.9$,  the isodensity curves have a smooth wavelike form away from the seamount and follow the two low-amplitude beams which are visible in the velocity field. Over the seamount, the isopycnals are more perturbed and the beams emanating from the sides are clearly visible. The free-surface perturbation is a smooth, low-amplitude wave except over the seamount where it is steeper.

For $\varepsilon=1$, a single beam is clearly formed over the seamount and the isodensity curves assume a cusp-like structure. This beam connects the two points of the topography for which the slope becomes equal to $1/\mu$. The free-surface is steeper and has a larger wave height in comparison with the previous case.

For $\varepsilon=1.5$, two fully developed beams are visible throughout the entire domain. They emanate from the sides of the seamount, and they intersect over its summit. The isopycnals have a cusp-like form for one beam and a step-like form for the other, and the free-surface follows a similar trend. Their vanishing widths and diverging amplitudes are physical manifestations of the singularity discussed in the previous subsection.

 Results concerning the shelf profile with $h_-/h+=0.5$ and $\varepsilon=2$ are shown in figure \ref{Fig:flowmaps_shelf}. Fully developed beams are also visible. The shoreward beams appear more intense in comparison with the oceanward ones. This calculation is qualitatively similar with the one found in figure 6 of GG07. The free-surface perturbation alternates between singular cusp-like forms and smooth steep forms, assuming larger amplitudes towards the deeper direction.

It must be noted that for $\varepsilon>1$, small-amplitude high-frequency oscillations are visible. These are attributed to the Gibbs-like phenomenon occurring due to the formation of singularities in the solution and are also present in other supercritical solutions \citep[GG07;][]{echeverri_peacock_2010}. Following GG07, we attenuate this phenomenon when plotting the fields by adjusting the Fourier-like summation \eqref{eq:RL_representation}, after the solution of the CMS, with the multiplication of a filter function of the form $\sin(n\pi/M)M/(n\pi)$ when $\varepsilon>1$. Finally, we also note that the cusp-like features of the isopycnal contours as $\varepsilon$ increases physically correspond to statically unstable density profiles, which would be regularized in the ocean via gravitational or shear instabilities. This is however outside the scope of this article.
\begin{figure}
\captionsetup{width=.9\linewidth}
\centerline{\includegraphics[width = 16cm]{Figures/Convergence/Fig_convergence.eps}  }
\caption{$(a)$-$(c)$: Decay of the uniform norms of the modal amplitudes, $\|\phi_n\|_{\infty}$,  $\|\phi_{n,x}\|_{\infty}$, $\|\phi_{n,xx}\|_{\infty}$, for the case of a bump topography with $\delta = 0.5$ and $\varepsilon = 0.9$, $1$, $1.5$. $(d)$-$(f)$: Absolute error \eqref{eq:energy_eq_avg_error} as a function of $(M,r)$.}
\label{fig:Error_decay}
\end{figure}

\begin{figure}
\captionsetup{width=.9\linewidth}
\centerline{
\includegraphics[width = 39pc]{Figures/Fields/Fig_flowmaps_bump2.eps}  }
\caption{Baroclinic velocity amplitude $|\mathbf{u}^{\#}|/U_0$, isodensity curves (full black lines) and free surface perturbation $2.5\times 10^3\times\eta$ (magenta lines)  at $t=0, T/4$ for the bump profile with $\delta = 0.5$  and $\varepsilon = 0.9$, $1$, $1.5$ ($L \approx 55$, $49$, $32$ km). The dashed black lines correspond to the rigid lid. The vertical magenta lines corresponds to a 2 cm amplitude.}
\label{Fig:flowmaps_bump}
\end{figure}

\begin{figure}
\captionsetup{width=.9\linewidth}
\centerline{
\includegraphics[width = 39pc]{Figures/Fields/Fig_flowmaps_shelf2.eps}  }
\caption{Baroclinic velocity amplitude $|\mathbf{u}^{\#}|/U_0$, isodensity curves (full black lines) and free surface perturbation $2.5\times 10^3\times\eta$ (magenta lines)  at $t=0, T/4$ for the shelf profile with $h_-/h_+=0.5$  and $\varepsilon = 2$. The dashed black lines correspond to the rigid lid. The vertical magenta line corresponds to a 2 cm amplitude.}
\label{Fig:flowmaps_shelf}
\end{figure}

\section{Conversion rates and comparison with previous approaches}

\subsection{Energy conversion rate}
We first examine the far-field energy flux or energy conversion rate for the topographies presented in section \ref{sec:Validation}. Our primary objective is to compare our calculations with existing analytical predictions, or past numerical computations.
We also provide quantitative estimates of the limits of applicability for analytical results, as well as a parametric formula for estimating energy flux radiated by supercritical topographies.

\subsubsection{Subcritical seamounts}
We examine the subcritical ($\varepsilon< 1$) and critical ($\varepsilon = 1$) cases first. In the regime $\varepsilon\ll 1$ and $\delta\ll 1$ (vanishing topographic slope and height), the WTA may be combined with a uniform barotropic tidal flow along the fluid domain. Using these assumptions together with the hydrostatic approximation, \cite{Llewellyn2002} derived a simple expression for the far-field energy flux in terms of the topography and the physical parameters of the problem. In the two-dimensional non-hydrostatic case, it is written as \citep{StLaurent2003}
\begin{align}\label{eq:ConversionStL}
\mathcal{C}\left[h\right] = F_0\sum\limits_{n=1}^{\infty}\kappa_n \hat{h}(k_n)\overline{\hat{h}}(k_n) \delta k,\quad\text{with}\quad F_0 = \frac{1}{2\pi}\overline{\rho}_0\frac{\left[(N^2-\omega^2)(\omega^2-f^2)\right]^{1/2}}{\omega}U^2h_0^2,
\end{align}
where $\hat{h}(\xi) = \int_{-\infty}^{+\infty}\exp(-ix\xi)h(x) dx$ is the Fourier transform of $h$, $k_n=n\pi/(\mu h_0)$ and $\delta k = k_n/n = \pi/(\mu h_0)$. From Eq.~\eqref{eq:ConversionStL} with $h=h_{\text{W}}$ and $h=h_{\text{G}}$, we obtain
\begin{align}
    \mathcal{C}\left[h_{\text{W}}\right]=F_0^2\frac{\pi^2}{4}\delta^2\frac{c^2 e^{-c}}{(1-e^{-c})^2},\quad
     \mathcal{C}\left[h_{\text{G}}\right]= 2\pi F_0\delta^2\left(\frac{L\pi}{h_0\mu}\right)^2\sum\limits_{n=1}^{\infty}n \exp{\left(-\frac{n\pi}{\mu h_0}\right)^2}
\end{align}
with $c=(3^{3/2}\pi/4)(h_0\delta/2L)$. For $h=h_{\text{B}}$, we calculated $\mathcal{C}\left[h_{\text{B}}\right]$ numerically using the Gauss-Kronrod quadrature. The total conversion rate in Eq.\,\eqref{eq:Conv_modal} calculated by means of the present CMS for $\varepsilon = 0.1$, $0.5$ and $1$, and $\delta\in[0.05,9]$, is compared with the  predictions of Eq.\ \eqref{eq:ConversionStL} in figure \ref{Fig:CG_comparison}.

Considering first the ``witch of Agnesi" and ``Gaussian" profiles in figure \ref{Fig:CG_comparison}$(a)$ and $(b)$, we observe that for every $\varepsilon$, both the calculated and analytical conversion rates increase proportionally to $\delta^2$ up to a global maximum value attained at a value $\delta_c({\varepsilon})$. The values of the maximum conversion rate and of $\delta_c$ both increase with $\varepsilon$. The WTA predictions based on Eq.~\eqref{eq:ConversionStL} are in adequate agreement with our calculations only for small values of $\delta$.
Arguably the most striking difference between our semi-analytical calculation and the analytical prediction is the accurate prediction, in our case, of ``non-radiating topographies'' \citep{Petrelis_et_al_2006,maas_2011}. These topographies manifest themselves as abrupt drops in the conversions rate for discrete values of $\delta$, up to $\delta \approx 1$, which indicate a sharp transition from weakly radiating to non-radiating.
We must note however that while the phenomenon happens for all topographies we tested, the conversion rate never strictly drops to zero (due to numerical uncertainties).
In the case of the bump profile (figure \ref{Fig:CG_comparison}$(c)$), oscillating patterns are more regular and a similar mismatch with analytical predictions occurs as $\delta$ increases, but some local minima cannot be considered as zeros and a small amount of energy escapes to infinity.

\begin{figure}
\captionsetup{width=.9\linewidth}
\centerline{
\includegraphics[width=16cm]{Fig_Csub_all.eps}}
\caption{Energy conversion rate as a function of $\delta$ for the $(a)$ witch of Agnesi, $(b)$ Gaussian and $(c)$ bump profiles with $\varepsilon=0.1$, $0.5$, $1$. Full lines correspond to the calculated conversion rates by the present method and dashed lines to the predictions by Eq. \eqref{eq:ConversionStL}. Their corresponding maxima are shown with full circles and crosses.}
\label{Fig:CG_comparison}
\end{figure}

A more thorough comparison of our calculations with WTA predictions as a function of $\delta$ and $\varepsilon$ is shown in the left panels of figure\,\ref{Fig:limits_model} by mapping the quantity $E_{\gamma} = |(C-\mathcal{C}[h_{\gamma}])/C|$, $\gamma = W$, $G$, $B$, for the three different topographies given in Eq.\ \eqref{eq:topos}.
A visual indication (black dashed lines) indicates our suggestion for the limit of applicability of analytical predictions when the differences between our computations and predictions are greater than $100\%$. Eq.~\eqref{eq:ConversionStL} should be used only for values of $\delta$ smaller than $\alpha\varepsilon$ with $\alpha$ being $0.8$ for the ``witch", $0.7$ for the ``Gaussian", and $0.4$ for the ``bump" profiles.

\subsubsection{Supercritical seamounts}
We now turn our attention to supercritical seamounts ($\varepsilon>1$). In this case, an analytical formula for seamounts of general smooth shape is lacking. Instead, we consider here the energy conversion rate due to knife-edge ridges, i.e., finite-height, zero-width vertical barriers. Indeed, using a Green's function approach, \cite{llewellyn_smith_young_2003} derived an expression for the energy flux, namely, \citep[see also][]{Garrett_Kunze_2007}
\begin{align}\label{eq:knife1}
\mathcal{C}^{\text{LSY}}_{\text{knife}} = 2\pi F_0\int_{0}^{\delta}Z\left(\frac{1-\cos\pi Z}{\cos\pi Z-\cos\pi\delta}\right)^{1/2} dZ.
\end{align}
where the integral in Eq.\ \eqref{eq:knife1} is computed here using the Gauss-Kronrod quadrature. At the same time, \cite{StLaurent2003} using a matching of modal solutions at the knife edge, obtained
\begin{align}\label{eq:knife2}
\mathcal{C}^{\text{StL}}_{\text{knife}} = F_0\sum_{n=1}^{\infty}n^{-1}a_n^2,
\end{align}
where $a_n$ are constant modal amplitudes obtained by the solution of a linear algebraic system. In order to compare our results with the above predictions, we calculate the energy flux for the three considered seamounts with $\varepsilon=2$, $5$, $10$ and $\delta\in[0.1,0.9]$.

Figure \ref{Fig:CE_super} shows that the conversion rate for supercritical topographies has a simpler dependence on $\delta$, in comparison with the subcritical case. It is in general increasing with $\delta$ with no points of low conversion rates.
The two knife-edge predictions, Eqs.\ \eqref{eq:knife1} and \eqref{eq:knife2}, differ slightly and give a more or less accurate prediction for $\varepsilon=2$. For larger values of $\varepsilon$ our calculations differ significantly from the knife-edge predictions for all seamounts; for example, $C$ becomes up to four times as large as $\mathcal{C}^{\text{LSY}}_{\text{knife}}$ for $\varepsilon=10$ and $\delta>0.5$.
For the ``bump" profile, a small hump in the conversion rate is observed for $\varepsilon=2$ and $0.4<\delta<0.7$ and only weak changes in the slope of $\mathcal{C}(\delta)$ are visible for $\varepsilon=5$ and $10$.
It must be noted that similar observations for the supercritical radiated energy flux have been observed with other approaches for generation \citep{Nycander2006, Balmforth_Neil_Peacock_2009, Zhao2015} and scattering problems \citep{Mathur2014} .

\begin{figure}
\captionsetup{width=.9\linewidth}
    \centering
    \includegraphics[width=16cm]{Fig_Csuper_all.eps}
    \caption{Energy conversion rates versus $\delta$ for the supercritical seamounts ($\varepsilon=2$, $5$, $10$), for $(a)$ witch of Agnesi, $(b)$ Gaussian and $(c)$ ``bump". The full and pointed black lines correspond to the knife-edge predictions $\mathcal{C}^{\text{LSY}}_{\text{knife}}$ and $\mathcal{C}^{\text{StL}}_{\text{knife}}$, respectively.}
\label{Fig:CE_super}
\end{figure}
We empirically define a parametric radiated energy flux $\mathcal{C}_{param}$, prescribed as
\begin{align}
\mathcal{C}_{param}/F_0 = p_0\varepsilon^{p_1}\delta^{p_2}\,,\label{eq_Cparam}
\end{align}
for values of $\delta$ between $0$ and $1$ and $\varepsilon$ between $1.1$ and $10$. The values for the parameters $p_0$, $p_1$, $p_2$ are summarized in Table \ref{tab_Cparam}. The influence of the topographic height (via $\delta$) increases with a power law slightly larger than the WTA model ($p_2>2$), and the radiated energy flux is approximately a linear function of $\varepsilon$.
\begin{table}[]
    \centering
    \begin{tabular}{c|c|c|c}
         seamount & $p_0$ & $p_1$& $p_2$\\ \hline
         witch & $8.19\times 10^3$ & 1.02 & 2.29 \\
         Gaussian & $8.60\times 10^3$ & 0.97 & 2.16 \\
         bump & $7.16\times 10^3$ & 0.93 & 2.03 \\
         \hline
         {\it mean} & $7.8\times 10^3$ & 1.0 & 2.1 \\
    \end{tabular}
    \caption{Values for the parameters in eq.~(\ref{eq_Cparam})} for different topographies.
    \label{tab_Cparam}
\end{table}
Like before, we compare our calculations with this parametric prediction as a function of $\delta$ and $\varepsilon$ by mapping the quantity $E_{\gamma} = |(C-\mathcal{C}_{param})/C|$ in the right panels of figure\,\ref{Fig:limits_model}.
On average, conversion rates for topographies with $\delta>0.3$ are well prescribed by this parametric approach. Trends in the conversion rates are not adequately reproduced for $\varepsilon\gg 1$ and $\delta<0.3$.
Similar results are obtained with the parametric model using the values of the parameters for each topography or the mean values in Table \ref{tab_Cparam}, which are thus a convenient set of values for any type of seamount.

\begin{figure}[htb]
\captionsetup{width=.9\linewidth}
\centering
\includegraphics[width=\linewidth]{Figures/map_error_alls_v2.eps}
\caption{Relative error maps between our computations for the subcritical cases (left panels) and supercritical cases (right panels) and the analytical or parametric models for $(a)$ the 'witch', $(b)$ the Gaussian and $(c)$ the bump topographies ($E_W$, $E_G$ and $E_B$ respectively), as a function of $\delta$ and $\varepsilon$.}
\label{Fig:limits_model}
\end{figure}

\subsubsection{Shelf profile}
In order to compare with the calculations of GG07, we take $h_-/h_+=0.5$ and calculate the shoreward and oceanward non-dimensional energy fluxes given by $C^{\prime}_- = C_-/Q^2 $ and $C^{\prime}_+ = C_+/Q^2 $ for $\varepsilon \in [0.1,4]$. Our results are plotted in figure \ref{Fig:CGrif}, together with results digitized from Fig.~1 of GG07. For $\varepsilon\lessapprox 0.4$, our calculations are in very good agreement capturing the oscillating pattern of $C^{\prime}_{\pm}$ and the appearance of local extrema. In the region  $0.4<\varepsilon<0.8$, we obtain the same increasing pattern, with our calculations giving slightly smaller values. Differences are more significant for $\varepsilon>0.8$. In our calculations, the shoreward and oceanward energy fluxes start to significantly diverge from each other at about $\varepsilon = 1.2$ while in GG07's results this happens earlier, at about $\varepsilon = 0.8$. Right after these points, the shoreward energy flux for both methods attains a local maximum. In our calculation, the shoreward energy flux starts increasing after a while, whereas GG07's calculation decreases. The oceanward energy flux increase rapidly in a small region ($\varepsilon\in[0.8,1]$ for GG07 and $\varepsilon\in[1,1.2]$ for the present CMS) and then follows an increasing pattern of smaller rate. This difference is possibly due to the hydrostatic approximation invoked in GG07 for the baroclinic equations, which is absent here.
\begin{figure}
\captionsetup{width=.9\linewidth}
\centerline{\includegraphics[scale = 1]{Figures/Fig_Grif_comparison.eps}}
\caption{Non-dimensional energy fluxes $C^{\prime}_- = C_-/Q^2$ and $C^{\prime}_+ = C_-/Q^2$  for the shelf profile with $h_-/h_+ = 0.5$. Circles and crosses correspond to the calculations of \cite{griffiths2007}.}
\label{Fig:CGrif}
\end{figure}

\subsection{Comparisons with the Green's function method and discussion on non-radiating topographies}

As stated before, semi-analytical models based on the Green's function method consider a uniform barotropic flow \citep{Petrelis_et_al_2006, echeverri_peacock_2010}, whereas in the present body-forcing approach we consider a topographically adjusted barotropic flow. Therefore, it is natural to assess the effect of the body-forcing approach by comparing our calculations with the ones based on the Green's function method, implemented here by means of the tool iTides \citep{iTides}.

In figure \ref{Fig:forcing_compare_max}$(a)$, we plot the normalized baroclinic velocity amplitude ($|\mathbf{u}^{\#}|/U_0$ in our case) for a subcritical Gaussian topography with $\varepsilon = 0.5$ and $\delta=0.3350$, for which a maximum energy flux is obtained. Differences are observed in the region over the summit of the seamount, with solution from iTides resulting in larger velocity amplitudes. In both calculations, we see a single wave beam on either side of the seamount that is of the same far-field velocity amplitude and parallel to characteristic lines. The difference between the calculated conversion rates is less than 2\%.

The supercritical calculations ($\varepsilon=1.5$, $\delta=0.5$) shown in figure \ref{Fig:forcing_compare_max}(b) feature stronger differences in the velocity fields. Our CMS calculation produces two beams, respectively initially upward and downward at the critical slope point, of clearly different far-field velocity amplitudes. On the other hand, iTides yields two beams of comparable far-field velocity amplitudes. In the CMS calculation, the beams emanating over the summit do not follow a straight path towards the sea surface and seabed. They slightly deviate from the characteristics and align with them only when $h_x$ becomes negligible.
This peculiar observation is due to the spatial dependency of the eigenvectors $\tilde{\kappa}_n(x)$ in the modal decomposition of the CMS problem (see Eq.\ \eqref{eq:RL_vertical_function}).
Indeed, the deviation from straight characteristics is more pronounced when $h_x$ is large.
Consequently, beams in the two calculations follow a visibly different path away from the topography. Nevertheless, the difference between the calculated far-field energy fluxes does not exceed 10\%.

For a non-radiating case ($\delta=0.7577$, $\varepsilon=0.5$) shown in figure \ref{Fig:forcing_compare_min}, the CMS calculation yields a trapped baroclinic response,  corresponding to a three-cell structure of the stream function in accordance with the results of \citet[figure 3b]{maas_2011}. On the other hand, the response calculated by the Green's function method has a simpler shape, the corresponding velocity being almost constant along the $z$-direction. Surprisingly, the solution from iTides is non-radiating but cannot be seen as solely the baroclinic response;  the wavefield shown resembles the adjusted hydrostatic barotropic field $\sqrt{U^{(0)}+V^{(0)}+W^{(0)}}$ of section \ref{sec:governing}\ref{subsec:BCsolution} (not shown here).
In both calculations, the normalized far-field energy flux is extremely small, of the order of $10^{-4}$, and the solutions are trapped around the topography.

\begin{figure}
\captionsetup{width=.9\linewidth}
\centerline{\includegraphics[width = 39pc]{Figures/Fields/Fig_submax_super2.eps}}
  \caption{Gaussian profile with $(a)$ $\varepsilon = 0.5$, $\delta = 0.3350$ and $(b)$ $\varepsilon = 1.5$, $\delta = 0.5$. Left panel: Baroclinic velocity amplitude $|\mathbf{u}^{\#}|/U_0$ and isodensity curves calculated by the present CMS. Right panel: baroclinic velocity amplitude calculated by the iTides-Green's function method. Characteristic lines are shown in orange.}
\label{Fig:forcing_compare_max}
\end{figure}
\begin{figure}
\captionsetup{width=.9\linewidth}
\centerline{\includegraphics[width = 39pc]{Figures/Fields/Fig_submin.eps}}
  \caption{Gaussian profile with $\delta=0.7577$ and $\varepsilon = 0.5$. $(a)$ Baroclinic velocity amplitude $|\mathbf{u}^{\#}|/U_0$ calculated by the present CMS. $(b)$ Baroclinic velocity amplitude calculated by iTides (Green's function method).}.
\label{Fig:forcing_compare_min}
\end{figure}
More generally, our computations on the energy conversion rate have shown that all considered topographies exhibit non- or weakly-radiating behavior for discrete values of $\delta$, at a fixed value of $\varepsilon$.
Such a behavior was noticed by \cite{Petrelis_et_al_2006} for the cases of a triangular and polynomial ridges. For these profiles, these so-called null points of the parameter space are already present within the WTA and their calculations show that they also exist for the complete problem. However, in the case of the asymptotically flat Gaussian and witch of Agnesi profiles, the WTA formula \eqref{eq:ConversionStL} does not predict null-points in contrast with our calculations.
For the bump profile, the WTA formula yields points of low energy conversion rate (no null points) but, as we have shown, these are different from the solutions of our CMS.
\cite{maas_2011} constructed non-radiating topographies, arguing that the absence of tidal conversion may be common. He also predicted that non-radiation is only possible for subcritical topographies with sufficiently large values of $\delta$. Our calculations confirm these facts.
It must be noticed that we did not find a valid approach to predict the values of the parameters $\varepsilon$ and $\delta$ for non-radiating cases, without solving the IT-generation problem. GG07 have discussed the weakly-radiating cases for the continental shelf in terms of an approximate phase variation of the first internal wave mode along the topography, but we tested their predictions with no success.
The only robust result we observed is that a non-radiating topography is also a non-scattering one. Indeed, we can solve the problem of mode-1 scattering by a Gaussian non-radiating topography ($\delta=0.7577, \varepsilon=0.5$), and our computations have shown that the transmitted wave field is only a mode-1 of same amplitude, while the reflected field is zero.

Finally, it is worth noting that our CMS computations can treat domains of different depths at infinity (shelves) for which a treatment with Green's function is currently lacking. The numerical efficiency has not been compared in detail with the tool iTides, but our implementation provides solutions faster.

\section{Conclusions}
A new semi-analytical model is developed that describes the generation of linear internal tides due to the interaction of the barotropic tide with arbitrary topography in two dimensions. The main novelty of the present approach is that it takes into account a barotropic flow that consistently varies over the bottom topography. The hydrodynamical problem is reformulated as a Coupled-Mode System (CMS) based on a local eigenfunction expansion of the unknown baroclinic response in terms of prescribed local vertical eigenfunctions and unknown horizontal modal amplitudes. Numerically solving the CMS yields the modal amplitudes, which are then used to determine the response fields and the energy conversion rate. A method to calculate the pressure-induced free-surface elevation is also proposed. The streamfunction formulation of our new CMS shows good convergence properties, being more efficient than approaches based on the velocity formulation, such as that of \citet{griffiths2007}.

We have compared our calculations with existing simplified predictions (WTA, knife-edge) of the radiated energy flux for three seamount profiles: two commonly used asymptotically flat topographies, the Gaussian and the witch of Agnesi, and a compactly supported one, the bump, studied here for the first time. The agreement between our calculations and analytical predictions is adequate provided $\delta$ is small enough for the WTA, and $\varepsilon$ close to 1 for the knife-edge approach. For $\varepsilon\leq 1$, our calculations show that there exist specific values of $\delta$ for which tidal conversion is either almost lacking for the Gaussian and witch of Agnesi profiles, or it is very weak for the bump profile. This effect is not captured by the WTA predictions. In the supercritical regime, for fixed $\delta$ the energy flux is increasing  with $\varepsilon$, in contrast with the knife-edge predictions. The poor agreement between analytical predictions and the present calculations suggests that the use of the former in global ocean models for the parameterization of energy conversion should be reviewed. For practical purposes, we propose a simplified dependence of the energy conversion rate in the supercritical regime, approximately proportional to $\varepsilon\delta^2$. We have also calculated the energy flux radiated due to a shelf profile. In this case, we have compared our calculations with the coupled-mode approach of \citet{griffiths2007}. We have observed strong differences between the two models, especially for $\varepsilon>1$, which are likely due to the hydrostatic approximation invoked in the latter work. Another set of comparisons was performed with the Green's function method, in which a uniform barotropic flow is considered. The CMS approach provides a more realistic flow above the topography, but yields more or less the same energy conversion rate as the Green's function method.
Nevertheless, our approach is more versatile in the sense that it can treat shelf-like domains for which a Green's function tool such as iTides is currently lacking. The treatment of trenches is also possible.

The present CMS can be extended in several interesting directions. The assumption of constant stratification can be lifted, at the cost of using a set of local basis functions which are calculated numerically, given a smooth stratification $N(z)$. Similarly, the set of local basis functions can be modified in order to take into account the linearized free-surface condition, instead of using the rigid-lid approximation. In this case, the added complexity comes from the fact that the local wavenumbers and eigenfunctions are determined in terms of a local transcendental equation (dispersion relation) \citep[see e.g.][]{Kelly2016}, which can nevertheless be solved efficiently \citep{PPA18}.  Finally, we note that following \citep[Chapter 2.1]{Papoutsellis}, it can be shown that the rapid convergence of the eigenfunction expansion remains after these extensions, thus we anticipate that realistic topographies and stratification will be efficiently taken into account. For realistic oceanic configurations, comparing satellite observations with calculated free-surface elevations will provide further validations of the method.

\section*{Acknowledgments}
The authors acknowledge the support of the IRT/STAE Foundation (funding for the post-doctorate contract of Ch.P), and also L.\ R.\ M.\ Maas for fruitful discussions. M.M. thanks Lucile Planes and Martin Benebig for testing computations with iTides. N.G.\ acknowledges the support of the Natural Sciences and Engineering Research Council of Canada (NSERC) [RGPIN-2015-03684].
\section*{Data statement}
Data available upon request.