[figure]style=plain,subcapbesideposition=top

# Experimental investigations of linear and nonlinear periodic traveling waves
in a viscous fluid conduit

Yifeng Mao<EMAIL_ADDRESS>Department of Applied Mathematics,
University of Colorado, Boulder, CO 80309, USA  Mark A. Hoefer  Department of
Applied Mathematics, University of Colorado, Boulder, CO 80309, USA

###### Abstract

Conduits generated by the buoyant dynamics between two miscible, Stokes fluids
with high viscosity contrast exhibit rich nonlinear wave dynamics. However,
little is known about the fundamental wave dispersion properties of the
medium. In the present work, a pump is used to inject a time-periodic flow
that results in the excitation of propagating small and large amplitude
periodic traveling waves along the conduit interface. This wavemaker problem
is used as a means to measure the linear and nonlinear dispersion relations
and corresponding periodic traveling wave profiles. Measurements are favorably
compared with predictions from a fully nonlinear, long-wave model (the conduit
equation) and the analytically computed linear dispersion relation for two-
Stokes flow. A critical frequency is observed, marking the threshold between
propagating and non-propagating (spatially decaying) waves. Measurements of
wave profiles and the wavenumber-frequency dispersion relation quantitatively
agree with wave solutions of the conduit equation. An upshift from the conduit
equation’s predicted critical frequency is observed and is hypothesized to be
explained by incorporating a weak recirculating flow into the full two-Stokes
flow model which we observe to be operable in the experiment. When the
boundary condition corresponds to the temporal profile of a nonlinear periodic
traveling wave solution of the conduit equation, weakly nonlinear and strongly
nonlinear, cnoidal-type waves are observed that quantitatively agree with the
conduit nonlinear dispersion relation and wave profiles. This wavemaker
problem is an important precursor to the experimental investigation of more
general boundary value problems in viscous fluid conduit nonlinear wave
dynamics.

_K_ eywords Stokesian dynamics, interfacial flows

##  1 Introduction

The wavemaker problem consists of waves launched from a boundary into a medium
at rest. It is a fundamental problem for inviscid flows with physical
applications such as surface waves in reservoirs or lakes produced by the
movement of a solid body (Chwang (1983b)). The Sommerfeld radiation condition
(Sommerfeld (1912); Sommerfeld (1949)) can be used to select a particular
solution that requires only outgoing waves at infinity. This boundary
condition has been applied to many boundary value problems in fluid dynamics
and elsewhere (see, e.g. Bühler (2014)). The classical theory considers
surface water waves forced by a wavemaker, represented by a vertical
impermeable plate which oscillates horizontally with a small displacement
(Havelock (1929); Kennard (1949); Chwang (1983a)). The wavemaker problem has
also been studied experimentally, primarily in water wave tanks exhibiting
various wave structures. For example, Goring & Raichlen (1980) successfully
produced clean cnoidal waves in shallow water with a piston wavemaker. Other
examples include the features and stability of two-dimensional periodic
surface wave patterns propagating on deep water with nearly persistent form
(Kimmoun et al. (1999); Hammack et al. (2005); Henderson et al. (2006)). An
isolated rogue wave was successfully created in a water tank by Chabchoub et
al. (2011).

This paper presents an alternative fluid medium in which to investigate the
wavemaker problem. The buoyancy driven, cylindrical free interface between two
miscible, Stokes fluids with high viscosity contrast results in fluid flow
through a pipe-like conduit. This scenario models viscously deformable media
in a variety of physical environments including magma in the earth’s mantle
(Scott & Stevenson (1984)) and channelized water flow beneath glaciers
(Stubblefield et al. (2020)). Beyond applications, this medium has been shown
to exhibit a rich variety of nonlinear wave dynamics. Experiments in
established conduits have been extensively studied in consideration of
solitons (Olson & Christensen (1986); Scott et al. (1986); Helfrich &
Whitehead (1990); Lowman et al. (2014)), periodic waves (Kumagai & Kurita
(1998); Olson & Christensen (1986)), solitary wave fission (Anderson et al.
(2019); Maiden et al. (2020)), dispersive shock waves (DSWs) (Maiden et al.
(2016)), and the interactions between solitons and DSWs (Maiden et al.
(2018)). Quantitative observational agreement with the conduit equation has
been achieved in certain regimes. But such agreement has largely been confined
to the long-wave, albeit strongly nonlinear, regime. Herein, we experimentally
investigate long and short periodic traveling waves in a viscous fluid
conduit.

Periodic wavetrains have already been observed and studied in viscous fluid
conduits. Some of the earliest apparently involved periodic waves by either
quickly increasing the injection rate from one constant value to another
constant value (Scott et al. (1986); Kumagai & Kurita (1998)) or, at the onset
of injection, the generation of periodic waves trailing a rising diapir on
plume of injected fluid (Olson & Christensen (1986)). However, since the
injection rate was held constant after the initial change, both scenarios
involve different dynamics than what we consider here. In the first case, we
now understand that a step up in injection rate from a nonzero previous rate
results in a dispersive shock wave consisting of a slowly varying periodic
traveling wave that ultimately returns to a constant injection rate, i.e. it
is transient (Maiden et al. (2016)). In the latter case, it is likely that an
oscillatory (absolute) instability developed due to sufficiently large
injection rate rather than the convective instability regime in which we
operate (d’OLCE et al. (2009)). In this work, we controllably create high
quality linear and nonlinear periodic traveling waves by use of a time-
periodic injection rate. The time-series are specifically designed by
evaluating exact periodic traveling wave solutions of the model conduit
equation.

The conduit equation

$A_{t}+(A^{2})_{z}-(A^{2}(A^{-1}A_{t})_{z})_{z}=0,$ (1.1)

modeling the dynamics of the normalized cross-sectional area $A$ of a conduit
at vertical height $z$ and time $t$, is an accurate partial differential
equation model of conduit interfacial waves that has been studied since the
1980s (Scott et al. (1986); Olson & Christensen (1986); Lowman & Hoefer
(2013); Maiden & Hoefer (2016); Johnson & Perkins (2020); Simpson et al.
(2006); Harris et al. (2006)). Derivation of the conduit equation follows from
the Stokes equations for two viscous fluids under a long-wave, high-viscosity
contrast assumption (Lowman & Hoefer (2013)). The resulting nonlinear,
dispersive equation is an expression of mass conservation that realizes a
maximal balance between buoyancy and viscous stress effects.

The major objective of this work is the reliable generation of both linear and
nonlinear, cnoidal-type waves from a wavemaker. Laboratory measurement
constrained by boundary control is realized by injection of an interior fluid
of lower viscosity and lower density into a more viscous exterior fluid using
a piston pump. Limitations are investigated and the experiment is compared
with the theory of the boundary value problem for the conduit equation. We
find that, though the exact periodic traveling wave solutions to the conduit
equation exhibit quantitative agreement with our experimental observations, a
deviation in the linear dispersion relation occurs for sufficiently short
waves. In order to describe the linear waves in the shorter-wavelength regime,
we obtain the exact solution for two-Stokes flow linear interfacial wave
propagation. Previous work by Yih (1963, 1967) and Hickox (1971) studied the
linear instability of the flow in a core-annular flow for the Navier-Stokes
equations. Within a parameter regime consistent with our experimental
conditions, the instability has been shown to be convective (d’OLCE et al.
(2009); SELVAM et al. (2009)). However, none of these works obtained an
analytical solution to the problem. Herein, we incorporate a model for
recirculating core-annular flow and obtain the exact solution to the Stokes
equations for azimuthally symmetric, small-amplitude perturbations. Asymptotic
analysis is also conducted in the limit of high viscosity contrast, long
waves, and small Reynolds number.

Outline of this work

The paper is structured as follows. First we provide some background material
on periodic traveling waves. Then we derive the exact linear dispersion
relation for two-Stokes flow in section 2. Analytical results for an idealized
recirculating flow with azimuthally symmetric linear interfacial waves and
their asymptotic expansions are presented. Section 3 describes linear
modulation theory applied to the initial-boundary value problem for both the
conduit equation and two-Stokes flow. In section 4, experimental methods and
the analysis of periodic waves in a viscous fluid conduit are provided.
Finally, we conclude in section 5.

The wavemaker problem and periodic traveling wave solutions

The conduit equation (1.1) obeys the scaling invariance

$A^{*}=\frac{A}{\bar{A}},\quad z^{*}=\bar{A}^{-1/2}z,\quad
t^{*}=\bar{A}^{1/2}t,$ (1.2)

where $\bar{A}>0$. Owing to this scaling invariance, we can always assume the
periodic traveling wave has unit mean. In this work, we study the wavemaker
problem for a viscous fluid conduit so that a time-dependent periodic wave
with amplitude $a$ and injection frequency $\omega_{0}$ is created at the
boundary and propagates into a quiescent medium. The wavemaker initial-
boundary value problem considers

$\displaystyle A(z,0)$ $\displaystyle=1,\quad z\geq 0,$ (1.3a) $\displaystyle
A(0,t)$ $\displaystyle=1+f(t),\quad t\geq 0,\quad f(t+T_{0})=f(t),\quad
t>0,\quad T_{0}=\frac{2\pi}{\omega_{0}},$ (1.3b)

where $f(t)$ is the periodic profile with frequency $\omega_{0}>0$,
$f(0)=f^{\prime}(0)=0$, $\int_{0}^{T_{0}}f(t)dt=0$.

Conduit periodic traveling wave solutions satisfy $A(z,t)=\phi(\theta)$,
$\theta=kz-\omega t$, $\phi(\theta+2\pi)=\phi(\theta)$ for
$\theta\in\mathbb{R}$ and a nonlinear ODE (Olson & Christensen (1986))

$-\omega\phi^{\prime}+k(\phi^{2})^{\prime}+\omega
k^{2}(\phi^{2}(\phi^{-1}\phi^{\prime})^{\prime})^{\prime}=0.$ (1.4)

Integrating (1.4) twice results in

$(\phi^{\prime})^{2}=g(\phi)\equiv-\frac{2}{k^{2}}\phi-\frac{2}{\omega
k}\phi^{2}\ln\phi+c_{1}+c_{2}\phi^{2},$ (1.5)

where $c_{1,2}$ are integration constants.

Upon linearization of (1.1) as $A(z,t)-1\propto\cos(kz-\omega t)$, the conduit
equation exhibits the bounded dispersion relation

$\omega=\frac{2k}{1+k^{2}}$ (1.6)

on a unit mean background. For the wavemaker problem, it is helpful to invert
the linear dispersion relation (1.6) and solve for $k(\omega)$. In order to
select the correct branch, we apply the radiation condition (positive group
velocity) to obtain (Bühler (2014))

$k(\omega)=\begin{cases}\frac{1-\sqrt{1-\omega^{2}}}{\omega},&0<\omega\leqslant
1\\\ \frac{1+i\sqrt{\omega^{2}-1}}{\omega},&\omega>1\end{cases},$ (1.7)

where wavenumbers $k$ for frequencies $0<\omega\leqslant 1$ are real, implying
steady wave propagation. For $\omega>1$, $k(\omega)$ exhibits a positive
imaginary part so that these waves spatially decay. The frequency
$\omega_{cr}=1$ is therefore the critical frequency, above which waves will
not propagate. The corresponding critical wavenumber $k_{cr}=1$ satisfies
$k_{cr}=k(\omega_{cr})$.

In the weakly nonlinear regime, we can obtain approximate periodic traveling
wave solutions when the amplitude is small and the wavenumber is sufficiently
far away from 0. The weakly nonlinear approximation with three harmonics takes
the form

$\phi(\theta;k,a)=1+\dfrac{a}{2}\cos(\theta)+\dfrac{1+k^{2}}{48k^{2}}a^{2}\cos(2\theta)+\dfrac{1+k^{2}}{1536k^{4}}a^{3}\cos(3\theta)+O(a^{4}),$
(1.8)

where $\theta=kz-\omega t$, the amplitude $0<a\ll 1$, the wavenumber $k\gg
a^{2}$, and the frequency $\omega$ is

$\omega(k,a)=\frac{2k}{1+k^{2}}+a^{2}\frac{1-8k^{2}}{48k(1+k^{2})}+O(a^{4}).$
(1.9)

Alternatively, we can express the weakly nonlinear expansions in terms of
$\omega$ and get

$\displaystyle\phi(\theta,\omega,a)$
$\displaystyle=1+\frac{a}{2}\cos(\theta)+\frac{a^{2}}{24\left(1-\sqrt{1-\omega^{2}}\right)}\cos(2\theta)$
$\displaystyle\hphantom{=1}+\frac{\left(-1-\sqrt{1-\omega^{2}}\right)a^{3}}{768\left(-2+\omega^{2}+2\sqrt{1-\omega^{2}}\right)}\cos(3\theta)+\cdots,$
(1.10a) $\displaystyle k(\omega,a)$
$\displaystyle=\frac{1-\sqrt{1-\omega^{2}}}{\omega}+\frac{\left(-9+\frac{7}{\sqrt{1-\omega^{2}}}\right)}{96\omega}a^{2}+\cdots.$
(1.10b)

The nonlinear dispersion relation is explicitly expressed in terms of three
roots $0<u_{1}<u_{2}<u_{3}$ of $g(\phi)$ in (1.5) (Lowman & Hoefer (2013)) and
phase velocity $U$ so that

$\exp
U=u_{1}^{\frac{u_{1}^{2}(u_{2}+u_{3})}{(u_{3}-u_{1})(u_{2}-u_{1})}}u_{2}^{-\frac{u_{2}^{2}(u_{3}+u_{1})}{(u_{3}-u_{2})(u_{2}-u_{1})}}u_{3}^{\frac{u_{3}^{2}(u_{2}+u_{1})}{(u_{3}-u_{1})(u_{3}-u_{2})}}.$
(1.11)

The wave amplitude, wavenumber and mean are given by

$a=u_{3}-u_{2},\quad
k=\pi\left(\int_{u_{2}}^{u_{3}}\frac{du}{\sqrt{g(u)}}\right)^{-1},\quad\bar{\phi}=\frac{k}{\pi}\int_{u_{2}}^{u_{3}}\frac{udu}{\sqrt{g(u)}}.$
(1.12)

The numerically computed conduit nonlinear dispersion relation $\omega=kU$ is
presented in figure 1, where the red dashed line corresponds to the edge of
the existence region.

Figure 1: Contour plot of numerically computed nonlinear conduit dispersion
relation $k(\omega,a)$. Beyond the red dashed line, no real $k(\omega,a)$ is
obtained.

##  2 Two-Stokes flow

Since the conduit equation (1.1) and its linear dispersion relation (1.6),
(1.7) are derived under long wavelength assumptions, we now seek a more
accurate description of linear waves by analyzing the free interface between
two Stokes fluids.

###  2.1 Primary recirculating flow

Interfacial waves generated by the buoyant dynamics between two miscible,
Stokes fluids with small viscosity ratio $\epsilon=\mu^{(i)}/\mu^{(e)}$ are
described by continuity equations for mass conservation and Stokes equations
for momentum conservation

$\displaystyle\nabla\cdot\mathbf{\tilde{u}}^{(i,e)}$ $\displaystyle=0,$ (2.1a)
$\displaystyle-\nabla\tilde{p}^{(i,e)}+\nabla\cdot\tilde{\sigma}^{(i,e)}$
$\displaystyle=0,$ (2.1b)

for the interior and exterior fluids, denoted by a superscript $i$ or $e$,
respectively. The velocity fields
$\mathbf{\tilde{u}}^{(i,e)}=(\tilde{u}_{\tilde{r}}^{(i,e)},\tilde{u}_{\tilde{z}}^{(i,e)})$
are given in cylindrical coordinates assuming azimuthal symmetry. The pressure
deviation from hydrostatic is $\tilde{p}^{(i,e)}$ and
$\tilde{\sigma}^{(i,e)}=\mu^{(i,e)}[\nabla\mathbf{\tilde{u}}+(\nabla\mathbf{\tilde{u}})^{T}]$
is the deviatoric stress tensor, where $\mu^{(i,e)}$ denotes the fluid
visocity. The kinematic boundary condition between the two fluids at
$\tilde{r}=\tilde{R}(z,t)$ incorporates the time dynamics

$\tilde{u}_{\tilde{r}}^{(i)}=\frac{\partial\tilde{R}}{\partial\tilde{t}}+\tilde{u}_{\tilde{z}}^{(i)}\frac{\partial\tilde{R}}{\partial\tilde{z}},\quad\tilde{r}=\tilde{R}(\tilde{z},\tilde{t}).$
(2.2)

The remaining boundary conditions include the continuity of normal and
tangential velocities

$[\mathbf{\tilde{u}}\cdot\mathbf{\hat{\tilde{n}}}_{c}]_{j}=[\mathbf{\tilde{u}}\cdot\mathbf{\hat{\tilde{t}}}_{c}]_{j}=0,\quad\tilde{r}=\tilde{R}(\tilde{z},\tilde{t}),$
(2.3)

and continuity of the normal and tangential stresses

$[\mathbf{\hat{\tilde{n}}}_{c}\cdot\mathbf{\tilde{T}}\cdot\mathbf{\hat{\tilde{n}}}_{c}]_{j}=[\mathbf{\hat{\tilde{T}}}_{c}\cdot\mathbf{\tilde{T}}\cdot\mathbf{\hat{\tilde{n}}}_{c}]_{j}=0,\quad\tilde{r}=\tilde{R}(\tilde{z},\tilde{t}),$
(2.4)

where $[]_{j}$ represents the difference in the exterior and interior fluid
quantities, $\mathbf{\hat{\tilde{n}}}_{c}$ is the outward pointing unit normal
to the interface and $\mathbf{\hat{\tilde{t}}}_{c}$ is the unit tangent to the
interface.
$\mathbf{\tilde{T}}^{(i,e)}=\tilde{\sigma}^{(i,e)}-\mathbf{I}\tilde{p}^{(i,e)}$
is the stress tensor and $\mathbf{I}$ denotes the identity operator. No-slip
and no-penetration boundary conditions at the outer wall $\tilde{r}=\tilde{D}$
require $\tilde{u}_{\tilde{z}}^{(e)}(\tilde{D})=0$ and
$\tilde{u}_{\tilde{r}}^{(e)}(\tilde{D})=0$. Boundary conditions along the
symmetry axis are imposed as

$\dfrac{\partial\tilde{u}_{\tilde{z}}^{(i)}}{\partial\tilde{r}}=0,\quad\tilde{u}_{\tilde{r}}^{(i)}=0,\quad\dfrac{\partial\tilde{p}^{(i)}}{\partial\tilde{r}}=0,\quad\text{at
}\tilde{r}=0.$ (2.5)

A nondimensionalization that leads to a maximal balance between buoyant and
viscous stress effects taking vertical length scale $L$, velocity scale $U$
and time scale $T$ as $\epsilon\to 0$ is (Lowman & Hoefer (2013))

$\displaystyle r$ $\displaystyle=\epsilon^{-1/2}\tilde{r}/L,\quad
z=\tilde{z}/L,\quad t=\tilde{t}/T,$ (2.6a) $\displaystyle\mathbf{u}^{(i,e)}$
$\displaystyle=\mathbf{\tilde{u}}^{(i,e)}/U,$ (2.6b) $\displaystyle p^{(i,e)}$
$\displaystyle=\frac{\tilde{p}^{(i,e)}-\tilde{p}_{0}}{\Pi},\quad\Pi=\mu^{(i)}U/L,$
(2.6c) $\displaystyle\sigma^{(i,e)}$
$\displaystyle=(L/\mu^{(i,e)}U)\tilde{\sigma}^{(i,e)},$ (2.6d)

where $\tilde{p}_{0}$ is a constant reference pressure. The vertical variation
is assumed to be weaker than the radial length scale. With constant conduit
radius $R_{0}$, the scales are defined as

$L=R_{0}/\sqrt{8\epsilon},\quad U=\frac{gR_{0}^{2}\Delta}{8\mu^{(i)}},\quad
T=L/U,\quad\Omega=\frac{2\pi}{T},$ (2.7)

where $\Delta=\rho^{(e)}-\rho^{(i)}$ is the difference in the fluid densities.
Substituting the above scalings into the dimensional equations represented by
the tilde notation, the nondimensional continuity and linear momentum
equations for the interior and exterior fluid are obtained (see Appendix A).

Steady flow requires only that the axial velocity and pressure have initial
values different from zero and the vertical velocity is independent of $z$,
i.e. $u_{r,0}^{(i)}=u_{r,0}^{(e)}=0,\frac{\partial u_{z,0}^{(i)}}{\partial
z}=\frac{\partial u_{z,0}^{(e)}}{\partial z}=0$. Collecting the governing
equations as well as boundary conditions, we obtain the solution for the
primary flow

$\displaystyle p_{0}^{(i)}$ $\displaystyle=(\lambda-1)z+\Lambda,$ (2.8a)
$\displaystyle u_{z,0}^{(i)}$
$\displaystyle=\dfrac{(\lambda-1)}{4}(r^{2}-\eta^{2})+\dfrac{\epsilon\lambda}{4}(\eta^{2}-D^{2})+\dfrac{\epsilon\eta^{2}}{2}\ln\frac{D}{\eta},$
(2.8b) $\displaystyle p_{0}^{(e)}$ $\displaystyle=\lambda z+\Lambda,$ (2.8c)
$\displaystyle u_{z,0}^{(e)}$
$\displaystyle=\epsilon\left(\dfrac{\lambda}{4}(r^{2}-D^{2})+\dfrac{\eta^{2}}{2}\ln\frac{D}{r}\right),$
(2.8d)

with external fluid pressure gradient $\lambda$, internal fluid pressure
gradient $\lambda-1$ and interfacial radius $\eta$. Pressure is determined up
to an arbitrary constant, so we take $\Lambda=0$. The vertical volumetric flow
rates of the interior and exterior flow through a cross section normal to the
$z$ axis are

$\displaystyle Q^{(i)}$ $\displaystyle=2\pi\int^{\eta}_{0}u_{z,0}^{(i)}(r)rdr$
$\displaystyle=\frac{\pi}{8}\left(1-\lambda+2\epsilon\lambda-2\epsilon\lambda\frac{D^{2}}{\eta^{2}}+4\epsilon\ln\frac{D}{\eta}\right)\eta^{4},$
(2.9a) $\displaystyle Q^{(e)}$
$\displaystyle=2\pi\int^{D}_{\eta}u_{z,0}^{(e)}(r)rdr$
$\displaystyle=-\frac{\pi}{8}\epsilon\left(2+\lambda-2(1+\lambda)\frac{D^{2}}{\eta^{2}}+\lambda\frac{D^{4}}{\eta^{4}}+4\ln\frac{D}{\eta}\right)\eta^{4}.$
(2.9b)

In the case of a motionless external fluid so that $\epsilon=0$ and
$\lambda=0$, the interior flow rate is $Q^{(i)}=\frac{\pi}{8}\eta^{4}$, which,
upon dimensionalization (2.7) results in the Poiseuille, pipe flow scaling law
$R_{0}^{4}=\frac{8\mu^{(i)}}{\pi g\Delta}Q_{0}$, where $Q_{0}$ is the interior
fluid injection rate. Note that the externel pressure gradient $\lambda$ is a
parameter which is constrained by the direction of fluid flow. Previous
studies that derived the conduit equation assumed a high viscosity contrast
($\epsilon\ll 1$) and zero external pressure gradient. However, in this work,
$\epsilon$ and $\lambda$, as well as $D$ are adapted to model our experiments
by introducing an idealized recirculating flow in which the interior fluid
rises upward dragging some of the exterior fluid with it, while the rest of
the exterior fluid flows downward. This situation dictates that

$\left(\frac{\eta}{D}\right)^{2}<\lambda<\frac{2\eta^{2}\ln(D/\eta)}{D^{2}-\eta^{2}}.$
(2.10)

So far, we have only specified the vertical flow directions for the interior
and exterior fluids of an infinitely long conduit. In order to model flow
recirculation, we need to appeal to the observed fluid dynamics in our
experimental apparatus. The rising interior fluid is observed to pool at the
very top of the fluid column, creating a less dense, lower viscosity fluid
layer above the exterior fluid (see figure 5a for a schematic). On the other
hand, the exterior fluid entrained by the rising interior fluid is observed to
recirculate downward and, upon reaching the bottom of the column, recirculates
upward. We model these dynamics by invoking the additional requirement of mass
conservation for the exterior fluid (i.e. $Q^{(e)}=0$ in (2.9b)), resulting in
the external pressure gradient

$\lambda=\frac{2\eta^{2}(D^{2}-\eta^{2}-2\eta^{2}\ln(D/\eta))}{(D^{2}-\eta^{2})^{2}}.$
(2.11)

The validity of the requirement (2.11) will be discussed in section 4.4 using
experimental measurements.

###  2.2 Linear dispersion relation

The linear dispersion relation for waves on a straight conduit in two-Stokes
flow is investigated by linearizing about the primary flow. Following the
procedure as introduced in Batchelor & Gill (1962) and Hickox (1971), first
order velocities and pressures are written as

$\displaystyle u_{1}^{(i)}$
$\displaystyle=[g_{z}(r),i\epsilon^{1/2}g_{r}(r)]e^{i(kz-\omega t)},$ (2.12)
$\displaystyle u_{1}^{(e)}$ $\displaystyle=[\epsilon
G_{z}(r),i\epsilon^{1/2}G_{r}(r)]e^{i(kz-\omega t)},$ (2.13) $\displaystyle
p_{1}^{(i)}$ $\displaystyle=h(r)e^{i(kz-\omega t)},$ (2.14) $\displaystyle
p_{1}^{(e)}$ $\displaystyle=H(r)e^{i(kz-\omega t)}.$ (2.15)

The $\epsilon$ scaling of the exterior and interior velocity perturbations are
chosen to reflect tangential stress balance and the kinematic boundary
condition. We assume a constant conduit radius $\eta$ subject to a small
disturbance so that

$R(z,t)=\eta+ae^{i(kz-\omega t)},\quad|a|\ll 1.$ (2.16)

Solving the governing equations directly and applying the boundary conditions,
the eigenfunctions $g_{r},g_{z},h$ and $G_{r},G_{z},H$ are obtained in terms
of a linear algebraic system given in Appendix B. Analytically, inverting the
6 by 6 matrix leads to overly burdensome expressions, so we either evaluate
the linear system numerically or consider its asymptotics for small
$\epsilon$. The linear dispersion relation is given by the kinematic boundary
condition as

$\omega(k)=ku_{z,0}^{(i)}(\eta)-g_{r}(\eta,k).$ (2.17)

In regions where (2.17) is monotonic in $k$, we can obtain $k(\omega)$ for
consideration of the wavemaker problem. In the vicinity of $\omega_{k}=0$, we
use the radiation condition to select the appropriate branch. Selecting
$\eta=\sqrt{8}$ for convenience when comparing with the $\epsilon\to 0$
dispersion relation (1.7), example dispersion curves for waves on a unit-mean
background are shown in figure 2. The maximum of $\text{Re}(k(\omega))$ occurs
at the critical frequency $\omega=\omega_{cr}$, defined by the group velocity
$c_{g}(\omega_{cr})=\frac{d\omega}{dk}(\omega_{cr})=0$, above which
$k(\omega)$ incorporates a nonzero imaginary part and the wave is spatially
damped. The corresponding critical wavenumber is $k_{cr}=k(\omega_{cr})$. In
the experiments reported in section 4.4, $\epsilon>0.01$ and the critical
frequency is subject to a deviation from 1.

(a)  (b)

Figure 2: Two-Stokes recirculating flow dispersion relation with $D/\eta=10$,
$\lambda$ in (2.11) and different $\epsilon$. Black squares represent the
critical frequency $\omega_{cr}$. (a) Real part of the dispersion relation
$\text{Re}(k(\omega))$. (b) Imaginary part $\text{Im}(k(\omega))$.

Consider the limit $\epsilon\to 0$. In order for the leading order behavior of
the two-Stokes system to be well approximated by the conduit dynamics, we
require a maximal balance further incorporating the outer wall’s no-slip
boundary condition so that

$D=\epsilon^{-1/2}d,\quad d=O(1).$ (2.18)

In this case, the corresponding external pressure gradient (2.11) becomes

$\lambda=\frac{16\epsilon\left(d^{2}+8\epsilon(\ln(8\epsilon)-2\ln(d))-8\epsilon\right)}{\left(d^{2}-8\epsilon\right)^{2}}=\epsilon\frac{16}{d^{2}}+O\left(\epsilon^{2}\ln(\epsilon)\right).$
(2.19)

This results in higher order corrections to the conduit unit-mean linear
dispersion relation as

$\omega(k)=\frac{2k}{1+k^{2}}+4\epsilon\ln\left(1/\epsilon\right)k\frac{1+2k^{2}}{(1+k^{2})^{2}}+\epsilon\omega_{2}(k)+\cdots,$
(2.20)

where $\omega_{2}(k)$ involves a complicated expression that depends on $d$
which can be found in appendix B. Following the assumptions of high viscosity
contrast and long-wavelengths, the conduit linear dispersion relation is
obtained at leading order in (2.20).

We provide the 3-term asymptotic expansion (2.20) to the two-Stokes linear
dispersion relation because the 2-term expansion has limited applicability.
For example, the absolute error in the wave frequency for
$0<\omega<\omega_{cr}$ between (2.20) and the full dispersion relation in
figure 3 shows that the second term does not significantly improve upon the
leading order term, even for $\epsilon\approx 10^{-6}$ but the 3-term
expansion does. Nevertheless, the 2-term expansion predicts an upshift in the
critical frequency $\omega_{cr}=1+3\epsilon\ln(1/\epsilon)$ and the critical
wavenumber $k_{cr}=1+\epsilon\ln(1/\epsilon)$. While the former is consistent
with the numerical calculation of the full dispersion in our experimental
regime where $\epsilon\ll 1$ (see figure 2), the latter is not.

Figure 3: Maximum absolute error in the comparison of exact and asymptotic
(2.20) solutions to the two-Stokes linear dispersion relation at $d=5$.

##  3 Linear modulation theory

To describe slow modulations of the linear periodic wave’s amplitude and
wavenumber for the wavemaker problem, we use modulation theory. The modulation
equation is the conservation of waves (Whitham (1974))

$k_{t}+\omega_{z}=0,$ (3.1)

where $\omega$ is the angular frequency and $k(\omega)$ is the local linear
dispersion relation. Translating the initial-boundary value problem for
$u(z,t)$ to an initial-boundary value problem for $k(z,t)$, we have

$\displaystyle\omega(z,0)$ $\displaystyle=0,\quad z\geq 0,$ (3.2a)
$\displaystyle\omega(0,t)$ $\displaystyle=\omega_{0},\quad t>0.$ (3.2b)

We seek a self-similar solution such that $\omega=\omega(\xi)$, $\xi=z/t$,
implying $\omega^{\prime}(\xi)=\xi$ or $k^{\prime}(\omega)=1/\xi$. The
solution takes the form

$\omega(\xi)=\begin{cases}\omega_{0},&0\leq\xi\leq c_{g}(\omega_{0})\\\
g(\xi),&c_{g}(\omega_{0})<\xi\leq c_{g}(0)\\\ 0,&\xi>c_{g}(0)\end{cases}$
(3.3)

where $\omega_{0}$ is the input frequency and $c_{g}(\omega_{0})$ represents
the group velocity for the wavenumber $k_{0}=k(\omega_{0})$. Comparison of the
solution (3.3) for the conduit and two-Stokes dispersion $k(\omega)$ is
depicted in figure 4. In particular, the conduit equation admits a plane wave
solution with constant frequency $\omega_{0}$ and wavenumber
$k_{0}=\frac{1-\sqrt{1-\omega_{0}^{2}}}{\omega_{0}}$ for $0\leq\xi<\xi_{0}$,
and a decaying oscillatory solution with
$\omega(\xi)=\frac{\sqrt{1-2\xi+\sqrt{1+4\xi}}}{\sqrt{2}}$ and
$k(\xi)=\frac{\sqrt{-1-\xi+\sqrt{1+4\xi}}}{\sqrt{\xi}}$ for
$\xi_{0}\leq\xi\leq 2$. The two regions are separated by the line
$\xi_{0}=1-\omega_{0}^{2}+\sqrt{1-\omega_{0}^{2}}$. For the two-Stokes
recirculating flow modulation solution, we have a similar profile for
$\omega(\xi)$ but the group velocities $c_{g}(\omega_{0}),c_{g}(0)$ are faster
than the corresponding conduit group velocities. These calculations suggest
that for an input $\omega_{0}$ smaller than the critical frequency, we can
generate a plane wave that will eventually propagate and fill the entire
vertical conduit. We now realize this in experiments described in the next
section.

Figure 4: Solution to the linear modulation equation conservation of waves
with an input frequency $\omega_{0}=0.8$ for the conduit equation (black
solid) and two-Stokes flow with $\epsilon=0.04,D/\eta=10$ and $\lambda$ at
(2.11) (blue dashed).

##  4 Experiment

###  4.1 Setup

The experimental setup in figure 5a is identical to those used by Anderson et
al. (2019) and Maiden et al. (2020). A conduit is formed in a square acrylic
column with dimensions $5\text{ cm}\times 5\text{ cm}\times 200\text{ cm}$.
The exterior fluid is clean glycerin of high viscosity and the interior fluid
consists of a miscible mixture of glycerin ($\sim 70\%$), water ($\sim 20\%$)
and food coloring ($\sim 10\%$), which is less dense and less viscous than the
exterior fluid. The nominal parameter values used in one of the experiments
(dataset 3 in sections 4.4, 4.5) are presented in table 1. Miscibility of the
two fluids implies negligible surface tension at the interface and the free
interface is sustained over the timescale of experiment because mass diffusion
occurs much more slowly than momentum diffusion.

A computer-controlled piston pump is used as a wavemaker to inject the
interior fluid into the extrusive fluid at the nozzle that then rises
buoyantly. Variation in the flow rate results in interfacial waves. Periodic
waves are created by a suitable periodic injection rate $Q^{(i)}(t)$:

$Q^{(i)}(t)=q(t),\quad t\geq 0,\quad q(t+T_{0})=q(t),\quad
T_{0}=\frac{2\pi}{\omega_{0}}.$ (4.1)

Data acquisition is performed using high resolution cameras with one (camera 1
in figure 5a) equipped with macro lenses at the injection nozzle for precise
data measurement and one to image the far field (camera 2) and capture fully
developed periodic wavetrains. Spatial calibration is achieved with a ruler
inside the column within camera view. Example large-amplitude periodic waves
with different injection frequencies and the same wave amplitude in time are
shown in figure 5b. For a relatively small frequency, the periodic wave
maintains its fully developed structure, i.e. it propagates. However, by
increasing the input frequency, the spatial decay of the amplitude is
observed. The middle wave is weakly, spatially damped, and the rightmost wave
is arrested at a larger frequency.

(a)

(b)

Figure 5: (a) Schematic of the experimental apparatus. (b) Waves
$\omega_{0}=(0.85,0.95,1.00)/\Omega$ with fixed amplitude $a=0.8$ but
increasing frequencies from left to right. The spatial wave is observed to be
periodic at the leftmost, but get damped in the middle and arrested at the
rightmost. Conduit edges are outlined by the red curves. Measured experimental
parameters follows from table 1.

$\mu^{(i)}$ $0.045$ Pa s $\mu^{(e)}$ $1.11$ Pa s $\epsilon$ 0.041 $\rho^{(i)}$
1.21 g cm-3 $\rho^{(e)}$ 1.26 g cm-3 $Q_{0}$ 3.00 cm3 min-1 $2R_{0}$ 0.36 cm
$2D_{0}$ 5.00 cm $Re^{(i)}$ 0.47 $Re^{(e)}$ 0.020 Camera resolution 375.60 pix
cm-1

Table 1: Example fluid properties measured in experiments: viscosities
$\mu^{(i,e)}$, viscosity ratio $\epsilon$, densities $\rho^{(i,e)}$,
background flow rate $Q_{0}$, associated conduit diameter $2R_{0}$ computed by
Poiseuille’s law, outer wall diameter $2D_{0}$, and Reynolds numbers
$Re^{(i,e)}$ for interior and exterior fluids.

###  4.2 Method

Camera images are processed to extract the conduit edges. A low-pass filter is
applied to reduce noise due to pixelation of the photograph and any impurities
in the exterior fluid. The conduit diameter is obtained by calculating the
number of pixels between the two edges and normalizing by the conduit mean
diameter. In each experimental trial, images are taken after the wave has
fully developed in the camera view and we identify $t=0$ with the initiation
of imaging. A $3$ Hz sampling rate for the close camera and a $1$ Hz sampling
rate for the full camera were used.

We present all dimensional quantities in tilde notations. The circular cross-
sectional area $\tilde{A}$ of the conduit is then modeled by the dimensional
conduit equation

$\tilde{A}_{\tilde{t}}+\frac{g\Delta}{8\pi\mu^{(i)}}\left(\tilde{A}^{2}\right)_{\tilde{z}}-\frac{\mu^{(e)}}{8\pi\mu^{(i)}}\left(\tilde{A}^{2}\left(\tilde{A}^{-1}\tilde{A}_{\tilde{t}}\right)_{\tilde{z}}\right)_{\tilde{z}}=0,$
(4.2)

where $\tilde{z},\tilde{t}$ follow from (2.6) and $\tilde{A}=\pi R_{0}^{2}A$.
Profile measurements of linear waves are conducted by scaling the data to unit
mean and fitting at every temporal and spatial snapshot a sinusoidal plane
wave

$A(\tilde{z},\tilde{t})=1+\frac{a}{2}\cos(\tilde{k}\tilde{z}-\tilde{\omega}\tilde{t}+\psi),$
(4.3)

where $\psi$ represents the phase. The relative $L^{2}$ norm error in the fit
is below $15\%$ in space and $10\%$ in time. To extract the wavenumber for a
trial, we compute the mean of $\tilde{k}$ across a range of times $\tilde{t}$;
similarly, the angular frequency is obtained by averaging $\tilde{\omega}$
over a range of heights $\tilde{z}$. The standard error in parameter
measurements is roughly as low as $1\%$ for $\tilde{k}$ and $0.5\%$ for
$\tilde{\omega}$, indicating that the fit (4.3) to the data results in a
reliable and accurate measurement of linear periodic waves. However, this
method is not suitable for waves beyond the linear regime. Instead, larger-
amplitude waves with macroscopic edge changes are processed by directly
splitting the waves into periods at each snapshot and computing the averaged
wavelengths and wave periods. The standard error for
$\tilde{k},\tilde{\omega}$ measurements with this technique is typically below
$5\%$.

To compare the experimental results with theory, we require the vertical
length, speed and time scales $L,U,T=L/U$ in (2.7) for nondimensionlization.
Upon generation of linear periodic waves, we nondimensionalize the wavenumber
and angular frequency data by using two fitting methods: one is to determine
$L_{c}$ and $U_{c}$ by fitting the dimensional data $\tilde{\omega},\tilde{k}$
to the conduit linear dispersion relation

$f(\tilde{k})=\dfrac{2U_{c}\tilde{k}}{1+L_{c}^{2}\tilde{k}^{2}},$ (4.4)

and the other is by fitting $\tilde{\omega},\tilde{k}$ with the two-Stokes
linear dispersion relation using the spatial $L_{S}$ and velocity $U_{S}$
scales, as well as $\epsilon,\lambda$ and $D$. Treating $\epsilon,\lambda$ and
$D$ as effective fitting parameters, the two-Stokes linear dispersion relation
is determined and fitted to experiments. We note that the fitting method is
sensitive to the initial guess for the fitting parameters. These two
techniques, which construct a direct relationship between experimental data
and theory, return reliable approximations of the dimensional scaling
coefficients $L$ and $U$. Differences in the two methods will be discussed in
section 4.4. The scales in our experiments from either fitting approach are
typically in the range

$L:0.3-0.4\text{ cm},\quad U:0.3-0.4\text{ cm/s},$ (4.5)

with camera resolutions of $250-400$ pixel cm-1. Although $L$ scales like
$R_{0}/\sqrt{\epsilon}$ for modeling long-wavelength wave dynamics, in our
experimental regime, it is comparable with the conduit diameter $2R_{0}$.
Nevertheless, for wavenumbers $0<k\lessapprox 1$, the wavelengths are on the
centimeter scale while the conduit diameter is on the scale of millimeters.
Prior to each experiment described below, around 10 trials of linear periodic
waves were generated for calibrating $L$ and $U$ according to the method
described here.

###  4.3 Definition of linear, weakly nonlinear and fully nonlinear regimes

Figure 6: Fourier mode of typical periodic waves in experiments. Waves A,B
with dominant $a_{1}$ can be approximated by a single sinewave. To describe
larger amplitude waves C,D, higher harmonic terms are needed. Waves E,F with
the largest amplitudes acquire non-negligible Fourier modes even at the fourth
or fifth order.

To distinguish the traveling periodic wave structures obtained in experiments,
we perform a preliminary analysis by expressing the spatial waves as a Fourier
series at a fixed time $t=t_{0}$:

$A(z,t_{0})=a_{0}+\sum^{\infty}_{n=1}\frac{a_{n}}{2}\cos(nz+\psi),$ (4.6)

where $a_{0}=1$ is the wave mean. Dividing the periodic traveling wave data
into periods and applying the Fourier cosine transform, we obtain the Fourier
components $a_{n}$. Fourier components of some selected periodic waves are
plotted in figure 6. Waves with negligible second-order and higher modes
$|a_{n}|\lessapprox 0.05$, $n=2,3,\dots$ are identified as lying in the linear
regime and can be well described by the linear theory. For larger-amplitude
waves $(0.4\lessapprox|a_{1}|\lessapprox 1.6)$, a single harmonic is
insufficient and higher order harmonic terms are needed. These waves are
considered to be in the weakly nonlinear regime. A clear distinction between
weakly nonlinear and fully nonlinear waves $(|a_{1}|\gtrapprox 1.6)$ occurs
when the fourth-order coefficient is sufficiently large, in which case the
weakly nonlinear approximation of three harmonics fails. We now assess each
wave class beginning with linear waves.

###  4.4 Linear periodic traveling waves

(a)

(b)

Figure 7: (a) Example experimental data of linear periodic traveling waves in
the viscous fluid conduit. Measured wave parameters are $a=0.12$, $k=1.34\pm
0.04$ rad/cm and $\omega=0.860\pm 0.002$ rad/s. (b) Linear wave data in
spatial domain fitted with sinusoidal waveform.

An example linear periodic traveling wave launched in a viscous fluid conduit
by the sinusoidal boundary condition
$Q(\tilde{t})=Q_{0}(1+\frac{a}{2}\cos(\tilde{\omega}\tilde{t}+\pi/2))^{2}$ for
$\tilde{t}>0$ is depicted in figure 7a. The wave traveling in the positive
$\tilde{z}$ direction is periodic in both time and space. Extracting the
spatial wave at a fixed time, we fit the periodic wave with a cosine function
in figure 7b. This result confirms the prediction from modulation theory that
a fully developed periodic wave is generated for input frequency less than the
critical value (see equations (1.7) and (3.3)). Three independent experiments
with different fluid properties resulting in three datasets of measured
wavenumbers $\tilde{k}$ and wave frequencies $\tilde{\omega}$ were collected.
In figure 8, we compare the dimensionless experimental observations of
$k(\omega)$ with the conduit linear dispersion relation (1.7). Across all
three datasets, good agreement between the experimental data and the
theoretical prediction is achieved for sufficiently long wavelengths
$\tilde{\zeta}=\frac{2\pi}{k}L_{c}\gtrapprox 2.3$ cm ($k<0.8$). The linear
dispersion relation of the conduit equation is quantitatively verified in this
regime. However, the experimental critical frequency $\omega_{cr}$ and
critical wavenumber $k_{cr}$ corresponding to the values between the last
datapoint of each dataset in figure 8 and the first point where the wave is
damped are subject to a significant shift from 1, the prediction from the
conduit equation. The relative discrepancies in $\omega_{cr}$ and $k_{cr}$ are
$(\Delta\omega_{cr},\Delta k_{cr})=(5\%,27\%),(7\%,33\%),(6\%,29\%)$,
respectively. It is notable that the observed critical frequency and
wavenumber vary across the three datasets where the fluid properties differ.
The interior flow rate is calculated using Poiseuille’s law and (2.7) for a
motionless external fluid so that $Q^{(i)}=\pi UR_{0}^{2}$, where $U=U_{c}$ is
the fitted velocity scale and $R_{0}$ is the conduit radius measured from
experiments. Datasets 1-3 have measured radius $R_{0}=0.19,0.17,0.18$ cm. As
expected for the asymptotic validity of the conduit equation, the vertical
wavelength is much larger than the horizontal diameter, i.e. $\tilde{\zeta}\gg
2R_{0}$. We obtain $Q^{(i)}=2.7,1.7,2.5$ ml/min for each dataset compared to
the nominal input flow rates $Q_{0}=3.0\pm 0.1,2.0\pm 0.1,3.0\pm 0.1$ ml/min,
respectively.

(a)  (b)

(c)

Figure 8: Fit of the experimental measurements $k(\omega)$ (dots) to the
conduit linear dispersion relation (black curve). Error bars take account of
the errors in measurements using (4.3) and errors in nondimensionalization.

To further interpret the linear periodic waves obtained in experiments and
investigate the shorter wave regime in the neighborhood of the critical
frequency, we compare the experimental results with the linear dispersion
relation for two-Stokes interfacial waves (2.17) subject to the assumption
$Q^{(e)}=0$ in (2.9b) so that $\lambda$ follows from (2.11). We fit the
experimental data $\tilde{k}(\tilde{\omega})$ in the subcritical regime with
the two-Stokes dispersion relation using effective fitting paramters
$\epsilon,D$, as well as the scales $L_{S}$ and $U_{S}$ to give an exactly
matched critical frequency $\omega_{cr}$. Each dataset was collected
independently and therefore possesses distinct fluid parameters although they
are all similar. Figure 9 shows the fits and table 2 reports the corresponding
$\omega_{cr}$ and $k_{cr}$. It is demonstrated that this fitting procedure
results in an upshift in the critical frequency and a downshift in the
critical wavenumber. Two-Stokes linear dispersion curves well describe the
dispersion relation $k(\omega)$ and provide improved predictions for the
critical values. Errors in the critical wavenumber $k_{cr}$ are caused by the
step size of the input $\omega$ in experiments, the sensitivity of $k_{cr}$ on
fluid parameters, and limitations in the fitting method. As $\omega$
approaches $\omega_{cr}$, $k^{\prime}(\omega)\to\infty$ so that a small change
in $\omega$ will lead to a substantial difference in $k$, leading to
difficulties in approximating $k_{cr}$. The predicted interior flow rates from
the fits using (2.9a) for dataset 1 to 3 are $Q^{(i)}=3.1,1.9,2.7$ ml/min,
reasonably consistent with the input $Q_{0}=3.0,2.0,3.0$ ml/min, and the
predicted interior radius obtained from (2.7) are $R_{0}=0.20,0.17,0.18$ cm,
respectively. The two-Stokes dispersion relation follows from the physical
scenario that we observed in experiments where the interior fluid flows up and
pools at the top over time while the exterior flow recirculates with a zero
net flow rate.

We have shown that the two-Stokes recirculating flow that includes an external
fluid pressure gradient and exterior flow mass conservation is a slightly more
accurate model of propagating linear waves generated by a wavemaker than the
long-wave conduit equation for describing relatively short waves between two
viscous fluids in the linear regime. It is notable, however, that the conduit
equation accurately reproduces the subcritical data.

(a)  (b)  (c)

Figure 9: Experimental results of linear periodic waves $k(\omega)$ (dots)
fitted with the two-Stokes dispersion relation assuming (2.11) (black curve).
Subfigures report predicted recirculating vertical flow velocities.

Dataset Measured Two-Stokes Relative error $\omega_{cr}$ $k_{cr}$
$\omega_{cr}$ $k_{cr}$ $\Delta k_{cr}$ 01 $1.04\pm 0.01$ $0.75\pm 0.01$ 1.04
0.88 15% 02 $1.08\pm 0.07$ $0.71\pm 0.05$ 1.08 0.90 21% 03 $1.07\pm 0.05$
$0.73\pm 0.03$ 1.07 0.90 19%

Table 2: Comparison of $\omega_{cr}$ and $k_{cr}$ between experimental
measurements and the two-Stokes linear dispersion relation assuming (2.11).
The fitting method requires an exactly correct $\omega_{cr}$. $\Delta k_{cr}$
reports the relative difference of $k_{cr}$.

###  4.5 Linear wave damping and arresting

(a)

(b)

(c)

Figure 10: (a) Surface and (b) contour plots of an experimental small
amplitude, damped, non-propagating wave when the injection frequency at the
boundary exceeds the critical value. Dimensional wave parameters are
$\tilde{\omega}=1.12\pm 0.03$ rad/s and $\tilde{k}=2.32\pm 0.04$ rad/cm with
scales $L=0.35$ cm, $U=0.40$ cm/s. (c) Amplitude decay (dotted blue) fitted
with $a\exp(-bz)$, $a=0.26\pm 0.01$, $b=0.027\pm 0.001$ (solid red).

(a) The conduit linear dispersion relation.

(b) The two-Stokes dispersion relation.

Figure 11: Comparison between experimental measurements (red dots) of the
exponential spatial decay rate (left), spatial frequency (right) and theory
(black line) for waves in the supercritical regime $\omega>\omega_{cr}$.

(a)  (b)  (c)

Figure 12: Investigation of a linear damping wave profile during propagation.
(a) Decay of amplitude is fitted with an exponential function. (b) Frequency
is also found to slightly decrease ($\sim 3\%$) in space. (c) The phase slope
is subject to a decay of $\sim 20\%$. Phase is measured from the temporal wave
at fixed $z$ and smoothed.

Linear periodic waves in the short-wave, supercritical regime when the
injection frequency $\omega$ exceeds the critical frequency $\omega_{cr}$
exhibit spatial amplitude damping and do not propagate unhindered. Figure 10a
and 10b demonstrate an example of a damped wave observed in a viscous fluid
conduit. The wave generated by a time harmonic injection rate at the nozzle
quickly decays spatially while maintaining periodicity in time at each fixed
spatial location. By measuring the spatial dependence of the amplitude of the
wave in figure 10c, we find that it is well fitted by an exponential function
in space. This agrees with our hypothesis derived from the radiation condition
(1.7) so that for complex dispersion $k=k_{\text{Re}}+ik_{\text{Im}}$, the
plane wave solution takes the form

$A(z,t)=ae^{-k_{\text{Im}}z}\cos(k_{\text{Re}}z+\omega t+\psi),$ (4.7)

where $k_{\text{Im}}$ is the spatial damping rate and $k=k(\omega)$ is
obtained from the linear dispersion relation in the supercritical regime.

The exact linear dispersion relation of the conduit equation and the two-
Stokes flow system obtained in section 2 are used to compare with the
experimental wave profiles. Here we report the data points in figure 11 as
collected using the same experimental setup for those in figures 8c and 9c;
therefore, identical parameters and nondimensionalization scales are used. No
new fitting is applied. For a total of 18 experimental trials, $k_{\text{Im}}$
is obtained as the spatial amplitude damping rate by fitting an exponential to
the wave amplitude, $\omega$ is the averaged frequency over a spatial window
size of $z\in[0,2]\to[0,0.7\text{cm}]$ (250 data points) measured through the
cosine fitting method described in section 4.2, and $k_{\text{Re}}$ is the
slope of the temporal wave phase over the same spatial window. Figure 11
reports the comparison of the linear dispersion relation in the supercritical
regime from experiments and theory. It is shown that the experimental results
closely match the two-Stokes theoretical predictions and the first data point
of the damping wave is also well approximated by the two-Stokes predicted
$\omega_{cr},k_{cr}$. The conduit dispersion relation also performs well at
describing $k_{\text{Im}}$ but significantly overshoots $k_{\text{Re}}$.

Damped waves are found to exhibit a spatially dependent frequency and
wavenumber as shown in figure 12. Through temporal data fitting with the
function $\phi(t)=a\cos(-\omega t+\psi)$, we obtain the amplitude $a$,
frequency $\omega$ and phase $\psi$, whose slope is the local wavenumber $k$.
With an exponentially decreasing amplitude consistent with linear theory,
$\omega$ and $k$ are found to spatially decrease. We are unable to explain
this spatially dependent frequency, and more notably, the decaying wavenumber.
For example, nonlinear effects are insufficient to explain these observations.

Nevertheless, the conduit equation again provides a reasonable description of
part of the supercritical data. We now focus solely on the conduit equation to
model nonlinear waves.

###  4.6 Periodic traveling waves in the weakly nonlinear regime

Larger-amplitude, periodic traveling waves in the weakly nonlinear regime are
distinguished from their linear counterpart by no longer maintaining symmetry
about the wave mean and cannot be described by a single trigonometric
harmonic. The waves are more accurately modeled by the weakly nonlinear
approximation (1.10a, 1.10b) of the conduit equation obtained via the
amplitude expansions $\phi\sim 1+a\phi_{1}+a^{2}\phi_{2}+a^{3}\phi_{3}$ and
$k\sim k_{0}+a^{2}k_{1}$, where $0<a\ll 1$ is an amplitude scale. An
experimental weakly nonlinear wave is shown in figure 13. Length and velocity
scales for nondimenionalization are $(L_{c},U_{c})=(0.33\pm 0.01\text{
cm},0.38\pm 0.01\text{ cm/sec})$. The dimensional spatial wave in figure 13a
is fitted with the three harmonic expansion (1.9) at $\tilde{t}=\tilde{t}_{0}$

$\phi(\tilde{z})=1+\dfrac{a_{z}}{2}\cos(\tilde{k}\tilde{z}+\psi_{1})+\dfrac{1+\tilde{k}^{2}}{48\tilde{k}^{2}}a_{z}^{2}\cos(2(\tilde{k}\tilde{z}+\psi_{1}))+\dfrac{1+\tilde{k}^{2}}{1536\tilde{k}^{4}}a_{z}^{3}\cos(3(\tilde{k}\tilde{z}+\psi_{1})),$
(4.8)

where $a_{z}=1.10\pm 0.01$, $k=0.50\pm 0.01$ $(\tilde{k}=1.52\pm 0.01$ rad/cm)
and $\psi_{1}=2.05\pm 0.01$. The approximate frequency is derived from figure
13c, where the dimensional temporal wave is approximated by (1.10a) at
$\tilde{z}=\tilde{z}_{0}$

$\begin{split}\phi(\tilde{t})=1&+\frac{a_{t}}{2}\cos(\tilde{\omega}\tilde{t}+\psi_{2})+\frac{a_{t}^{2}}{24\left(1-\sqrt{1-\tilde{\omega}^{2}}\right)}\cos(2(\tilde{\omega}\tilde{t}+\psi_{2}))\\\
&+\frac{\left(-1-\sqrt{1-\tilde{\omega}^{2}}\right)a_{t}^{3}}{768\left(-2+\tilde{\omega}^{2}+2\sqrt{1-\tilde{\omega}^{2}}\right)}\cos(3(\tilde{\omega}\tilde{t}+\psi_{2})),\end{split}$
(4.9)

with $a_{t}=1.06\pm 0.02$, $\omega=0.78\pm 0.01$ ($\tilde{\omega}=0.90\pm
0.01$ rad/s) and $\psi_{2}=-2.97\pm 0.03$. Fits exhibit a maximum absolute
error of $0.11$ and have larger discrepancies at the wave peaks. Adding higher
order harmonic terms yield an even better approximation.

A total of 11 trials were measured in the weakly nonlinear regime with the
same input amplitude and various input frequencies. Wave profiles are measured
as follows. Fitting the spatial and temporal waves with the weakly nonlinear
expansions (4.8, 4.9), we obtain dimensionless
$k=\tilde{k}L,\omega=\tilde{\omega}\frac{L}{U}$ as well as
$a=\frac{a_{z}+a_{t}}{2}$ for those 11 trials as the red circles in figure 14.
Experimental results are then compared with the weakly nonlinear dispersion
relation (1.10b) and the linear dispersion relation (1.6) in figure 14.
Although weak nonlinearity does not lead to a significant deviation from the
linear dispersion relation, we have shown that the corresponding wave profiles
require amplitude-dependent higher harmonics to quantitatively model the
experimental observations. We have quantitatively validated the weakly
nonlinear approximation for measured conduit periodic traveling waves of
sufficiently large amplitude and demonstrated the necessity of including
higher order harmonic terms in this regime.

(a)

(b)

(c)

(d)

Figure 13: Measured periodic waves (blue dots) and the weakly nonlinear
approximation (red solid lines). Error bars are represented by gray dashed
lines and the horizontal black dashed line is the wave mean. (a) The spatial
wave at a fixed $t$ fitted with expansion (4.8). (b) Percentage relative error
in space. (c) Wave in the time domain at a fixed $z$ fitted with expansion
(4.9). (d) Percentage relative error in time. Figure 14: Nondimensional
experimental wave profiles $k(\omega)$ (circles) compared with the weakly
nonlinear dispersion relation (1.10b) at $a=1.04$ (black solid line) and the
linear dispersion relation (blue dashed line).

###  4.7 Cnoidal-like nonlinear periodic waves

(a)  (b)

(c)  (d)

Figure 15: Averaged experimental wave profiles over one wavelength or period
(blue dots) compared with linear (green dashed), weakly nonlinear (red dash-
dotted) approximations and cnoidal-like solutions (black solid) in both
spatial (a, c) and temporal (b, d) domains. (a-b) The extracted wave
parameters are $(a_{\text{avg}},\omega,k)=(2.20,0.44,0.24)$ (trial 1 in figure
16). (c-d) The extracted parameters are
$(a_{\text{avg}},\omega,k)=(1.71,0.50,0.28)$ (trial 3 in figure 16).

We now investigate periodic traveling wave solutions with large amplitudes in
the fully nonlinear regime. We refer to these waves as cnoidal-like, owing to
their similarities to cnoidal shallow water waves and solutions of the
Korteweg–de Vries equation (Korteweg & de Vries (1895)). First, we generate
linear periodic waves for calibrating the scales $L_{c}=0.26\pm 0.01$ cm and
$U_{c}=0.25\pm 0.01$ cm/sec. We then utilize a numerical database of pre-
computed conduit cnoidal-like solutions to experimentally generate and analyze
the observed large-amplitude periodic waves. Waves are generated by varying
the injection flow rate according to

$Q(\tilde{t})=Q_{0}\phi^{2}(-\tilde{\omega}\tilde{t}),$ (4.10)

where $\phi$ is a cnoidal-like solution with unit mean satisfying (1.5).
Measured waves are processed by dividing each space and time window into
spatial and temporal periods. The average of all processed periods and wave
profiles are then reported. Figure 15 displays two example wave profiles,
separately in space and time. By direct comparison, it is clear that the
observed large amplitude waves cannot be described by linear or weakly
nonlinear approximations. Instead, we compare the observed waves with cnoidal-
like solutions to the conduit equation. The amplitude is obtained as
$a_{z}=\max(\phi(z))-\min(\phi(z))$ for waves in the spatial domain, and
$a_{t}=\max(\phi(t))-\min(\phi(t))$ for waves in the temporal domain. $\omega$
and $k$ are measured by the relation
$\omega=\frac{2\pi}{\overline{T_{0}}},k=\frac{2\pi}{\overline{\zeta}}$ with
nondimensional averaged time period $\overline{T_{0}}$ and averaged wavelength
$\overline{\zeta}$. The set of unit-mean cnoidal-like solutions is a two-
parameter family. Taking $(a_{z},k)$ as the parameters for spatial waves and
$(a_{t},\omega)$ as the parameters for temporal waves, the corresponding
cnoidal-like solutions are compared with the experiments in figure 15. The
cnoidal-like wave profiles are barely distinguishable from the observed wave
profiles. We stress that we have not performed any fitting in this procedure,
other than the independent determination of the scales $L_{c}$ and $U_{c}$
from linear wave dispersion measurements. The cnoidal-like solution of the
conduit equation is uniquely selected from the two-parameter family using the
extracted $(a_{z},k)$ in the spatial case or $(a_{t},\omega)$ in the temporal
case. No profile fitting was performed. However, discrepancies in spatial and
temporal profile measurements lead to $(a_{z},k)$ and $(a_{t},\omega)$
belonging to two slightly different cnoidal-like solutions for the same trial.
We report the average of the two amplitudes
$a_{\text{avg}}=\frac{a_{z}+a_{t}}{2}$ in table 3. To see how close the
solutions are, measurements of wave parameters for 15 trials are reported in
table 3 and figure 16.

(a)  (b)

Figure 16: (a) Comparison of the conduit equation nonlinear dispersion
relation for cnoidal-like solutions. Data points are from the trials reported
in table 3. Utilizing measured $(a_{\text{avg}},\omega)$ as the parameters,
experimental wavenumbers (blue circles) are compared with the numerically
computed nonlinear dispersion relation $k(a_{\text{avg}},\omega)$ (black
squares). (b) Percentage relative error in $k$ for each trial.

Trial $a_{\text{avg}}$ $\omega$ (rad) $k$ (rad) $k(a_{\text{avg}},\omega)$
Relative error 01 $2.20\pm 0.03$ $0.44\pm 0.02$ $0.24\pm 0.01$ 0.22 $9\%$ 02
$2.13\pm 0.01$ $0.45\pm 0.02$ $0.24\pm 0.01$ 0.23 $8\%$ 03 $1.71\pm 0.01$
$0.50\pm 0.01$ $0.28\pm 0.01$ 0.26 $6\%$ 04 $1.42\pm 0.01$ $0.54\pm 0.01$
$0.30\pm 0.01$ 0.29 $4\%$ 05 $1.30\pm 0.01$ $0.62\pm 0.02$ $0.37\pm 0.01$ 0.36
$2\%$ 06 $1.27\pm 0.01$ $0.61\pm 0.03$ $0.36\pm 0.01$ 0.35 $4\%$ 07 $1.26\pm
0.01$ $0.61\pm 0.02$ $0.36\pm 0.01$ 0.34 $5\%$ 08 $1.25\pm 0.01$ $0.60\pm
0.01$ $0.35\pm 0.01$ 0.34 $3\%$ 09 $1.24\pm 0.01$ $0.60\pm 0.03$ $0.34\pm
0.01$ 0.34 $3\%$ 10 $1.23\pm 0.01$ $0.59\pm 0.01$ $0.33\pm 0.01$ 0.33 $2\%$ 11
$1.22\pm 0.01$ $0.57\pm 0.01$ $0.32\pm 0.01$ 0.31 $2\%$ 12 $1.18\pm 0.01$
$0.54\pm 0.02$ $0.30\pm 0.01$ 0.29 $5\%$ 13 $1.15\pm 0.01$ $0.52\pm 0.02$
$0.28\pm 0.01$ 0.27 $4\%$ 14 $1.09\pm 0.01$ $0.59\pm 0.01$ $0.33\pm 0.01$ 0.32
$2\%$ 15 $0.89\pm 0.01$ $0.60\pm 0.01$ $0.35\pm 0.01$ 0.34 $2\%$

Table 3: Nonlinear periodic wave measurements. Columns 2-4: measured wave
parameters. Column 5: numerically computed solution $k$ at measured
$(a_{\text{avg}},\omega)$. Column 6: relative error between experimental $k$
and cnoidal-like wave solutions. Trials 1 and 3 are depicted in figure 15. All
trials are shown in figure 16.

The nonlinear dispersion relation $k(\omega,a_{\text{avg}})$ is compared
between experiments and cnoidal-like solutions in figure 16. Agreement is
quantitatively achieved for all but three trials within $5\%$ relative error,
with at most $10\%$ error otherwise. These measurements are across a range of
wave amplitudes $a_{\text{avg}}\in(0.89,2.2)$ and frequencies
$\omega\in(0.44,0.62)$. The trials in figure 16 are numbered from highest
amplitude to lowest, as shown in table 3. Notably, the experimentally measured
wavenumbers exceed the exact values across all 15 trials. Although a different
class of waves, we remark that previous measurements of large-amplitude
soliton solutions in (Olson & Christensen (1986)) found that the soliton speed
was consistently underpredicted by conduit equation soliton solutions. We
hypothesize that inclusion of the weakly recirculating flow in the model,
could improve upon the already good predictions.

##  5 Discussion and conclusion

The main result of this study is the reliable generation of periodic traveling
waves in a viscous fluid conduit across a wide range of amplitudes and long to
short wavelengths that are well-modeled by the conduit equation. This provides
further evidence of the viscous fluid conduit as an accessible and flexible
laboratory system for the study of nonlinear wave dynamics with an accurate
partial differential equation model in the conduit equation. In order to
quantitatively model shorter small amplitude waves, we found it helpful to
derive the full two-Stokes linear dispersion relation that accounts for the
actual conditions in our experiment. Radial variation in the vertical velocity
allows for a non-negligible exterior pressure gradient and recirculating flow.
Together with an outer wall, these two main factors give rise to an upshifted
critical frequency for the transition from propagating to non-propagating,
spatially damped waves. This upshift from the conduit equation’s critical
frequency in the short-wave regime is observed and verified quantitatively.
The two-Stokes recirculating flow system with exterior flow mass conservation
is therefore considered to be a more precise model in which we successfully
extend the linear dispersive properties of the model to a shorter wave regime.
The two-Stokes linear dispersion relation is exactly obtained in terms of
Bessel functions. Nevertheless, the conduit equation remains an excellent
model of propagating linear waves.

Another key result of this experimental investigation is the existence and
stability of periodic traveling waves in the fluid conduit wavemaker problem.
Low frequency periodic waves generated at the boundary with a time-periodic
injection rate are observed to coherently propagate as traveling waves.
Arrested wave propagation for waves generated above the critical injection
frequency was observed and found to quantitatively match the two-Stokes
dispersion relation for real frequencies and complex wavenumbers. These
observations are explained in terms of the radiation condition and linear wave
modulation theory.

In addition to presenting a careful analysis of linear periodic waves, our
work also contributes to an understanding of conduit periodic traveling wave
solutions in the weakly nonlinear and strongly nonlinear regimes. Nonlinear
periodic waves are found to require more Fourier harmonics with strongly
nonlinear waves well-described by cnoidal-like solutions of the conduit
equation. Quantitative agreement is achieved between experimental results and
conduit cnoidal-like solutions for both the nonlinear dispersion relation and
the wave profiles. The characterization of fully nonlinear periodic traveling
waves and their reliable generation in a viscous fluid conduit experiment hold
promise for future studies of more general boundary value problems in
dispersive hydrodynamics (El & Hoefer (2016)), including the generalized
Riemann problem and the generation of breathers and breather trains.

Funding. This work was supported by the National Science Foundation (grant
number DMS-1816934).

Declaration of interests. The authors report no conflict of interest.

##  Appendix A Two-Stokes flow and nondimensionalization

The governing equations (2.1a-2.1b) for continuity and momentum conservation
will now be nondimensionalized according to (2.6). Equations for the interior
fluid take the nondimensional form

$\displaystyle\dfrac{1}{r}\dfrac{\partial}{\partial
r}(ru_{r}^{(i)})+\epsilon^{1/2}\dfrac{\partial u_{z}^{(i)}}{\partial z}$
$\displaystyle=0,$ (A.1a) $\displaystyle-\epsilon^{-3/2}\dfrac{\partial
p^{(i)}}{\partial
r}+\nabla^{2}u_{r}^{(i)}-\epsilon^{-1}\dfrac{u_{r}^{(i)}}{r^{2}}$
$\displaystyle=0,$ (A.1b) $\displaystyle-\epsilon^{-1}\dfrac{\partial
p^{(i)}}{\partial z}+\nabla^{2}u_{z}^{(i)}$ $\displaystyle=0,$ (A.1c)

where $\nabla^{2}=\epsilon^{-1}\frac{1}{r}\frac{\partial}{\partial
r}(r\frac{\partial}{\partial r})+\frac{\partial^{2}}{\partial z^{2}}.$ For the
exterior fluid, the equations are

$\displaystyle\dfrac{1}{r}\dfrac{\partial}{\partial
r}(ru_{r}^{(e)})+\epsilon^{1/2}\dfrac{\partial u_{z}^{(e)}}{\partial z}$
$\displaystyle=0,$ (A.2a) $\displaystyle-\epsilon^{-1/2}\dfrac{\partial
p^{(e)}}{\partial
r}+\nabla^{2}u_{r}^{(e)}-\epsilon^{-1}\dfrac{u_{r}^{(e)}}{r^{2}}$
$\displaystyle=0,$ (A.2b) $\displaystyle-\dfrac{\partial p^{(e)}}{\partial
z}+\nabla^{2}u_{z}^{(e)}$ $\displaystyle=0.$ (A.2c)

The kinematic boundary condition is

$u_{r}^{(i)}=\epsilon^{1/2}\left(\dfrac{\partial R}{\partial
t}+u_{z}^{(i)}\dfrac{\partial R}{\partial z}\right),\quad r=R(z,t).$ (A.3)

Equations for continuity of the velocity components are

$\displaystyle(u_{r}^{(i)}-u_{r}^{(e)})$
$\displaystyle=\epsilon^{1/2}\dfrac{\partial R}{\partial
z}(u_{z}^{(i)}-u_{z}^{(e)}),\quad r=R(z,t),$ (A.4a)
$\displaystyle(u_{z}^{(e)}-u_{z}^{(i)})$
$\displaystyle=\epsilon^{1/2}\dfrac{\partial R}{\partial
z}(u_{r}^{(i)}-u_{r}^{(e)}),\quad r=R(z,t).$ (A.4b)

The nondimensional stress balance conditions in the normal and tangential
directions are, respectively,

$\displaystyle\left[-\|\mathbf{n}_{c}\|^{2}P+\kappa\left(\sigma_{rr}-2\epsilon^{1/2}\dfrac{\partial
R}{\partial z}\sigma_{rz}+\epsilon\left(\dfrac{\partial R}{\partial
z}\right)^{2}\sigma_{zz}\right)\right]_{j}$ $\displaystyle=0,\quad r=R(z,t),$
(A.5a)
$\displaystyle\left[\kappa\left\\{-\left(1-\epsilon\left(\dfrac{\partial
R}{\partial z}\right)^{2}\right)\sigma_{rz}+\epsilon^{1/2}\dfrac{\partial
R}{\partial z}(\sigma_{zz}-\sigma_{rr})\right\\}\right]_{j}$
$\displaystyle=0,\quad r=R(z,t),$ (A.5b)

where $\kappa^{(e)}=\epsilon^{-1},\kappa^{(i)}=1$ are the fluid-specific
coefficients, $\mathbf{n}_{c}=\begin{bmatrix}1\\\
-\epsilon^{1/2}\frac{\partial R}{\partial z}\end{bmatrix}$ is the normal
vector for the conduit, and $P$ is the scaled, absolute pressure with
$P^{(i,e)}=p^{(i,e)}-\frac{\rho^{(i,e)}z}{\rho^{(e)}-\rho^{(i)}}$. The
normalized deviatoric stress tensor is

$\sigma^{(i,e)}=\begin{bmatrix}\sigma_{rr}^{(i,e)}&\sigma_{rz}^{(i,e)}\\\
\sigma_{zr}^{(i,e)}&\sigma_{zz}^{(i,e)}\end{bmatrix}=\epsilon^{1/2}\begin{bmatrix}2\dfrac{\partial
u_{r}^{(i,e)}}{\partial r}&\dfrac{\partial u_{z}^{(i,e)}}{\partial
r}+\epsilon^{1/2}\dfrac{\partial u_{r}^{(i,e)}}{\partial z}\\\ \dfrac{\partial
u_{z}^{(i,e)}}{\partial r}+\epsilon^{1/2}\dfrac{\partial
u_{r}^{(i,e)}}{\partial z}&2\epsilon^{1/2}\dfrac{\partial
u_{z}^{(i,e)}}{\partial z}\end{bmatrix}.$ (A.6)

The above nondimensional governing equations associated with the interfacial
boundary conditions, as well as boundary conditions along the symmetry axis
$u_{r}^{(i)}=\frac{\partial u_{z}^{(i)}}{\partial r}=\frac{\partial
p^{(i)}}{\partial r}=0$ at $r=0$ and the outer wall
$u_{r}^{(e)}=u_{z}^{(e)}=0$ at $r=D$ provide a complete description of the
dynamics of two-Stokes interfacial flow.

##  Appendix B Two-Stokes exact linear dispersion relation

In this section, we derive the solutions for the eigenfunctions in the two-
Stokes flow linearization (2.12-2.15) and compute the exact linear dispersion
relation. The differential equations for two-Stokes flow disturbances are

$\displaystyle\dfrac{dg_{r}(r)}{dr}+\frac{g_{r}(r)}{r}+kg_{z}(r)$
$\displaystyle=0,$ (B.1a) $\displaystyle\epsilon
g^{\prime\prime}_{r}+\frac{\epsilon
g^{\prime}_{r}}{r}-k^{2}\epsilon^{2}g_{r}-\frac{\epsilon
g_{r}}{r^{2}}+ih^{\prime}(r)$ $\displaystyle=0,$ (B.1b) $\displaystyle
g^{\prime\prime}_{z}+\frac{g^{\prime}_{z}}{r}-k^{2}\epsilon g_{z}-ikh(r)$
$\displaystyle=0$ (B.1c)

for the interior fluid on the interval $0\leq r\leq\eta$ and

$\displaystyle\dfrac{dG_{r}(r)}{dr}+\frac{G_{r}(r)}{r}+k\epsilon G_{z}(r)$
$\displaystyle=0,$ (B.2a) $\displaystyle
G_{r}^{\prime\prime}+\frac{G_{r}^{\prime}}{r}-k^{2}\epsilon
G_{r}-\frac{G_{r}}{r^{2}}+iH^{\prime}(r)$ $\displaystyle=0,$ (B.2b)
$\displaystyle G_{z}^{\prime\prime}+\frac{G_{z}^{\prime}}{r}-k^{2}\epsilon
G_{z}-ikH(r)$ $\displaystyle=0$ (B.2c)

for the exterior fluid on the interval $\eta\leq r\leq D$. Along with the
interfacial boundary conditions expanded about $r=R(z,t)=\eta+ae^{i(kz-\omega
t)}$:

$\displaystyle ku_{z,0}^{(i)}(\eta)-g_{r}(\eta)-\omega$ $\displaystyle=0,$
(B.3a) $\displaystyle\epsilon
ku_{z,0}^{(e)}(\eta)-ku_{z,0}^{(i)}(\eta)+g_{r}(\eta)-G_{r}(\eta)$
$\displaystyle=0,$ (B.3b) $\displaystyle\epsilon\dfrac{\partial
u_{z,0}^{(e)}}{\partial r}(\eta)-\dfrac{\partial u_{z,0}^{(i)}}{\partial
r}(\eta)-g_{z}(\eta)+\epsilon G_{z}(\eta)$ $\displaystyle=0,$ (B.3c)
$\displaystyle-2i\epsilon k\dfrac{\partial u_{z,0}^{(e)}}{\partial
r}(\eta)+2i\epsilon k\dfrac{\partial u_{z,0}^{(i)}}{\partial
r}(\eta)-2i\epsilon
g^{\prime}_{r}(\eta)+2iG^{\prime}_{r}(\eta)+h(\eta)-\epsilon H(\eta)$
$\displaystyle=0,$ (B.3d)
$\displaystyle-\dfrac{\partial^{2}u_{z,0}^{(e)}}{\partial
r^{2}}(\eta)+\dfrac{\partial^{2}u_{z,0}^{(i)}}{\partial r^{2}}(\eta)-\epsilon
kg_{r}(\eta)+kG_{r}(\eta)+g^{\prime}_{z}(\eta)-G^{\prime}_{z}(\eta)$
$\displaystyle=0,$ (B.3e)

and $g_{r}(0)=0,g_{z}^{\prime}(0)=0,G_{r}(D)=0,G_{z}(D)=0$, the equations can
be solved in terms of modified Bessel functions. The internal flow admits the
general solution:

$\displaystyle g_{r}(r)$
$\displaystyle=a_{1}K_{1}\left(\epsilon^{1/2}kr\right)+a_{2}I_{1}\left(\epsilon^{1/2}kr\right)+a_{3}rK_{0}\left(\epsilon^{1/2}kr\right)+a_{4}rI_{0}\left(\epsilon^{1/2}kr\right),$
(B.4a) $\displaystyle g_{z}(r)$
$\displaystyle=a_{1}\epsilon^{1/2}K_{0}\left(\epsilon^{1/2}kr\right)-a_{2}\epsilon^{1/2}I_{0}\left(\epsilon^{1/2}kr\right)-\frac{2a_{3}}{k}K_{0}\left(\epsilon^{1/2}kr\right)+a_{3}\epsilon^{1/2}rK_{1}\left(\epsilon^{1/2}kr\right)$
$\displaystyle\hphantom{=a_{1}}-\frac{2a_{4}}{k}I_{0}\left(\epsilon^{1/2}kr\right)-a_{4}\epsilon^{1/2}rI_{1}\left(\epsilon^{1/2}kr\right),$
(B.4b) $\displaystyle ih(r)$
$\displaystyle=-2\epsilon\left(a_{3}K_{0}\left(\epsilon^{1/2}kr\right)+a_{4}I_{0}\left(\epsilon^{1/2}kr\right)\right).$
(B.4c)

The solution set is linearly independent since the Wronskian
$W=-\frac{4k^{2}\epsilon}{r^{2}}$ is nonzero. Symmetry along $r=0$ requires
$a_{1}=a_{3}=0$. Similarly, we obtain the general solution for the extrusive
velocities and pressure as

$\displaystyle G_{r}(r)$
$\displaystyle=b_{1}K_{1}\left(\epsilon^{1/2}kr\right)+b_{2}I_{1}\left(\epsilon^{1/2}kr\right)+b_{3}rK_{0}\left(\epsilon^{1/2}kr\right)+b_{4}rI_{0}\left(\epsilon^{1/2}kr\right),$
(B.5a) $\displaystyle G_{z}(r)$
$\displaystyle=\frac{1}{\epsilon}\Bigg{[}b_{1}\epsilon^{1/2}K_{0}\left(\epsilon^{1/2}kr\right)-b_{2}\epsilon^{1/2}I_{0}\left(\epsilon^{1/2}kr\right)-\frac{2b_{3}}{k}K_{0}\left(\epsilon^{1/2}kr\right)+b_{3}\epsilon^{1/2}rK_{1}\left(\epsilon^{1/2}kr\right)$
$\displaystyle\hphantom{=a_{1}}-\frac{2b_{4}}{k}I_{0}\left(\epsilon^{1/2}kr\right)-b_{4}\epsilon^{1/2}rI_{1}\left(\epsilon^{1/2}kr\right)\Bigg{]},$
(B.5b) $\displaystyle iH(r)$
$\displaystyle=-2\left(b_{3}K_{0}\left(\epsilon^{1/2}kr\right)+b_{4}I_{0}\left(\epsilon^{1/2}kr\right)\right).$
(B.5c)

Applying the interfacial boundary conditions and no-slip, no-penetration outer
conditions, a linear system is derived for the coefficients $a_{2,4},b_{1-4}$:

$A\begin{pmatrix}a_{2}\\\ a_{4}\\\ b_{1}\\\ b_{2}\\\ b_{3}\\\
b_{4}\end{pmatrix}=\begin{pmatrix}0\\\
-\frac{\eta}{2}(1-\epsilon)(1-\lambda)\\\ 0\\\ 1\\\ 0\\\ 0\end{pmatrix},$
(B.6)

where

$A=\begin{pmatrix}I_{1}&I_{0}\eta&-K_{1}&-I_{1}&-K_{0}\eta&-I_{0}\eta\\\
A_{21}&A_{22}&A_{23}&A_{24}&A_{25}&A_{26}\\\
A_{31}&A_{32}&A_{33}&A_{34}&A_{35}&A_{36}\\\
A_{41}&A_{42}&A_{43}&A_{44}&A_{45},&A_{46}\\\
0&0&K_{1}(\epsilon^{1/2}Dk)&I_{1}(\epsilon^{1/2}Dk)&DK_{0}(\epsilon^{1/2}Dk)&DI_{0}(\epsilon^{1/2}Dk)\\\
0&0&K_{0}(\epsilon^{1/2}Dk)&-I_{0}(\epsilon^{1/2}Dk)&A_{65}&A_{66}\end{pmatrix}$
(B.7)

with

$\displaystyle\begin{aligned}
A_{21}&=\epsilon^{1/2}I_{0},\quad&&A_{22}=2\frac{I_{0}}{k}+\epsilon^{1/2}I_{1}\eta,\quad&&A_{23}=\epsilon^{1/2}K_{0},\\\
A_{24}&=-\epsilon^{1/2}I_{0},\quad&&A_{25}=-2\frac{K_{0}}{k}+\epsilon^{1/2}K_{1}\eta,\quad&&A_{26}=-2\frac{I_{0}}{k}-\epsilon^{1/2}I_{1}\eta,\\\
A_{31}&=\frac{2}{\eta}\left(-\epsilon
I_{1}+\epsilon^{3/2}I_{0}k\eta\right),\quad&&A_{32}=2\epsilon^{3/2}I_{1}k\eta,\quad&&A_{33}=\frac{2}{\eta}\left(K_{1}+\epsilon^{1/2}K_{0}k\eta\right),\\\
A_{34}&=\frac{2}{\eta}\left(I_{1}-\epsilon^{1/2}I_{0}k\eta\right),\quad&&A_{35}=2\epsilon^{1/2}K_{1}k\eta,\quad&&A_{36}=-2\epsilon^{1/2}I_{1}k\eta,\\\
A_{41}&=-2\epsilon I_{1}k,\quad&&A_{42}=-2\epsilon^{1/2}I_{1}-2\epsilon
I_{0}k\eta,\quad&&A_{43}=2K_{1}k,\\\
A_{44}&=2I_{1}k,\quad&&A_{45}=-\frac{2}{\epsilon^{1/2}}K_{1}+2K_{0}k\eta,\quad&&A_{46}=\frac{2I_{1}}{\epsilon^{1/2}}+2I_{0}k\eta,\end{aligned}$
(B.8) $\displaystyle
A_{65}=\frac{1}{\epsilon^{1/2}k}\left(-2K_{0}(\epsilon^{1/2}Dk)+\epsilon^{1/2}DkK_{1}(\epsilon^{1/2}Dk)\right),$
$\displaystyle
A_{66}=\frac{1}{\epsilon^{1/2}k}\left(-2I_{0}(\epsilon^{1/2}Dk)-\epsilon^{1/2}DkI_{1}(\epsilon^{1/2}Dk)\right).$

and $I_{0}=I_{0}(\epsilon^{1/2}\eta k)$, $I_{1}=I_{1}(\epsilon^{1/2}\eta k)$,
$K_{0}=K_{0}(\epsilon^{1/2}\eta k)$, $K_{1}=K_{1}(\epsilon^{1/2}\eta k)$.
Varying $\epsilon,D$ leads to changes in the dispersion curves from the
conduit dispersion relation as well as changes in the critical frequency
$\omega_{cr}$, critical wavenumber $k_{cr}$ and inflection point $k_{inf}$
(see figure 17).

(a)  (b)  (c)

Figure 17: Investigation of (a) the critical frequency $\omega_{cr}$, (b) the
critical wavenumber $k_{cr}$, and (c) the inflection point $k_{inf}$ at
selected parameters $\epsilon,D$ for recirculating flow (2.11).

Following the asymptotic behavior $D=\epsilon^{-1/2}d$ and
$\lambda=\epsilon\frac{16}{d^{2}}+O(\epsilon^{2}\ln(\epsilon))$, where
$d=O(1)$, we expand the modified Bessel functions in the linear system (B.6)
for $\epsilon\to 0$. We seek solutions to the system (B.6) with the asymptotic
expansions

$v=\begin{pmatrix}a_{2}\\\ a_{4}\\\ b_{1}\\\ b_{2}\\\ b_{3}\\\
b_{4}\end{pmatrix}=\begin{pmatrix}\epsilon^{-3/2}\left(a_{20}+\epsilon\ln(\epsilon)a_{21}+\epsilon
a_{22}+\epsilon^{2}\ln(\epsilon)a_{23}+\epsilon^{2}a_{24}+\dots\right)\\\
\epsilon^{-1}\left(a_{40}+\epsilon\ln(\epsilon)a_{41}+\epsilon
a_{42}+\epsilon^{2}\ln(\epsilon)a_{43}+\epsilon^{2}a_{44}+\dots\right)\\\
\epsilon^{1/2}\left(b_{10}+\epsilon\ln(\epsilon)b_{11}+\epsilon
b_{12}+\epsilon^{2}\ln(\epsilon)b_{13}+\epsilon^{2}b_{14}+\dots\right)\\\
\epsilon^{1/2}\left(b_{20}+\epsilon\ln(\epsilon)b_{21}+\epsilon
b_{22}+\epsilon^{2}\ln(\epsilon)b_{23}+\epsilon^{2}b_{24}+\dots\right)\\\
\epsilon\left(b_{30}+\epsilon\ln(\epsilon)b_{31}+\epsilon
b_{32}+\epsilon^{2}\ln(\epsilon)b_{33}+\epsilon^{2}b_{34}+\dots\right)\\\
\epsilon\left(b_{40}+\epsilon\ln(\epsilon)b_{41}+\epsilon
b_{42}+\epsilon^{2}\ln(\epsilon)b_{43}+\epsilon^{2}b_{44}+\dots\right)\end{pmatrix}.$
(B.9)

Collecting (B.6) with (B.9) at each order and applying the kinematic boundary
condition (B.3a), we achieve the asymptotic expression (2.20) for the linear
dispersion relation with

$\displaystyle\omega_{2}(k)$
$\displaystyle=4k\left(-1+\ln\left(\frac{d}{2\sqrt{2}}\right)\right)-\frac{\sqrt{2}k(1+2k^{2})}{1+k^{2}}b_{20}-2\sqrt{2}b_{40}$
$\displaystyle\hphantom{=}+\frac{2k}{3d^{2}(1+k^{2})^{2}}\Bigg{[}d^{2}\left(6\Gamma\left(k^{4}-2k^{2}-1\right)+k^{4}(\ln(8)-19)-6k^{2}(4+\ln(2))-6-\ln(8)\right)$
$\displaystyle\hphantom{=}+6d^{2}\left(k^{4}-2k^{2}-1\right)\ln(k)-48\left(k^{2}+1\right)\Bigg{]},$
(B 10) where $\Gamma$ is the Euler–Mascheroni constant (0.5772…) and the
coefficients are $\displaystyle b_{10}$
$\displaystyle=-\frac{4\sqrt{2}k^{2}}{1+k^{2}},$ (B.10a) $\displaystyle
b_{20}$
$\displaystyle=\frac{I_{0}(dk)\left(K_{1}(dk)\left(d^{2}kb_{30}+2b_{10}\right)+dkb_{10}K_{0}(dk)\right)+dkI_{1}(dk)(b_{10}K_{1}(dk)+db_{30}K_{0}(dk))}{dkI_{0}(dk){}^{2}-2I_{1}(dk)I_{0}(dk)-dkI_{1}(dk){}^{2}},$
(B.10b) $\displaystyle b_{30}$
$\displaystyle=-\frac{\sqrt{2}k(1+3k^{2})}{1+k^{2}},$ (B.10c) $\displaystyle
b_{40}$
$\displaystyle=-\frac{kI_{0}(dk)(b_{10}K_{1}(dk)+db_{30}K_{0}(dk))+I_{1}(dk)(K_{0}(dk)(kb_{10}-2b_{30})+dkb_{30}K_{1}(dk))}{dkI_{0}(dk){}^{2}-2I_{1}(dk)I_{0}(dk)-dkI_{1}(dk){}^{2}}.$
(B.10d)

## References

  * Anderson et al. (2019) Anderson, Dalton V, Maiden, Michelle D & Hoefer, Mark A 2019 Controlling dispersive hydrodynamic wavebreaking in a viscous fluid conduit. Physical Review Fluids 4 (7), 074804.
  * Batchelor & Gill (1962) Batchelor, GK & Gill, A Eo 1962 Analysis of the stability of axisymmetric jets. Journal of Fluid Mechanics 14 (4), 529–551.
  * Bühler (2014) Bühler, Oliver 2014 Waves and mean flows. Cambridge University Press.
  * Chabchoub et al. (2011) Chabchoub, A, Hoffmann, NP & Akhmediev, Nail 2011 Rogue wave observation in a water wave tank. Physical Review Letters 106 (20), 204502.
  * Chwang (1983a) Chwang, Allen T 1983a Nonlinear hydrodynamic pressure on an accelerating plate. The Physics of Fluids 26 (2), 383–387.
  * Chwang (1983b) Chwang, Allen T 1983b A porous-wavemaker theory. Journal of Fluid Mechanics 132, 395–406.
  * d’OLCE et al. (2009) d’OLCE, M., MARTIN, J., RAKOTOMALALA, N., SALIN, D. & TALON, L. 2009 Convective/absolute instability in miscible core-annular flow. part 1: Experiments. Journal of Fluid Mechanics 618, 305–322.
  * El & Hoefer (2016) El, GA & Hoefer, MA 2016 Dispersive shock waves and modulation theory. Physica D: Nonlinear Phenomena 333, 11–65.
  * Goring & Raichlen (1980) Goring, Derek & Raichlen, Fredric 1980 The generation of long waves in the laboratory. In Coastal Engineering 1980, pp. 763–783.
  * Hammack et al. (2005) Hammack, Joseph L, Henderson, Diane M & Segur, Harvey 2005 Progressive waves with persistent two-dimensional surface patterns in deep water. Journal of Fluid Mechanics 532, 1–52.
  * Harris et al. (2006) Harris, Shirley E, Clarkson, Peter A & others 2006 Painlevé analysis and similarity reductions for the magma equation. SIGMA. Symmetry, Integrability and Geometry: Methods and Applications 2, 068.
  * Havelock (1929) Havelock, Th H 1929 Lix. forced surface-waves on water. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 8 (51), 569–576.
  * Helfrich & Whitehead (1990) Helfrich, Karl R & Whitehead, John A 1990 Solitary waves on conduits of buoyant fluid in a more viscous fluid. Geophysical & Astrophysical Fluid Dynamics 51 (1-4), 35–52.
  * Henderson et al. (2006) Henderson, Diane M, Patterson, Matthew S & Segur, Harvey 2006 On the laboratory generation of two-dimensional, progressive, surface waves of nearly permanent form on deep water. Journal of Fluid Mechanics 559, 413–427.
  * Hickox (1971) Hickox, Charles E 1971 Instability due to viscosity and density stratification in axisymmetric pipe flow. The physics of Fluids 14 (2), 251–262.
  * Johnson & Perkins (2020) Johnson, Mathew A & Perkins, Wesley R 2020 Modulational instability of viscous fluid conduit periodic waves. SIAM Journal on Mathematical Analysis 52 (1), 277–305.
  * Kennard (1949) Kennard, EH 1949 Generation of surface waves by a moving partition. Quarterly of applied mathematics 7 (3), 303–312.
  * Kimmoun et al. (1999) Kimmoun, Olivier, Branger, Hubert & Kharif, Christian 1999 On short-crested waves: experimental and analytical investigations. European Journal of Mechanics-B/Fluids 18 (5), 889–930.
  * Korteweg & de Vries (1895) Korteweg, Diederik Johannes & de Vries, Gustav 1895 Xli. on the change of form of long waves advancing in a rectangular canal, and on a new type of long stationary waves. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 39 (240), 422–443.
  * Kumagai & Kurita (1998) Kumagai, Ichiro & Kurita, Kei 1998 Interaction of periodic wave trains in magma conduits. Physics of the earth and planetary interiors 107 (1-3), 131–141.
  * Lowman & Hoefer (2013) Lowman, Nicholas K & Hoefer, Mark A 2013 Dispersive hydrodynamics in viscous fluid conduits. Physical Review E 88 (2), 023016.
  * Lowman et al. (2014) Lowman, Nicholas K, Hoefer, Mark A & El, Gennady A 2014 Interactions of large amplitude solitary waves in viscous fluid conduits. Journal of fluid mechanics 750, 372–384.
  * Maiden et al. (2018) Maiden, Michelle D, Anderson, Dalton V, Franco, Nevil A, El, Gennady A & Hoefer, Mark A 2018 Solitonic dispersive hydrodynamics: theory and observation. Physical review letters 120 (14), 144101.
  * Maiden et al. (2020) Maiden, Michelle D, Franco, Nevil A, Webb, EG, El, GA & Hoefer, MA 2020 Solitary wave fission of a large disturbance in a viscous fluid conduit. Journal of Fluid Mechanics 883.
  * Maiden & Hoefer (2016) Maiden, Michelle D & Hoefer, MARK A 2016 Modulations of viscous fluid conduit periodic waves. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 472 (2196), 20160533.
  * Maiden et al. (2016) Maiden, Michelle D, Lowman, Nicholas K, Anderson, Dalton V, Schubert, Marika E & Hoefer, Mark A 2016 Observation of dispersive shock waves, solitons, and their interactions in viscous fluid conduits. Physical review letters 116 (17), 174501.
  * Olson & Christensen (1986) Olson, Peter & Christensen, Ulrich 1986 Solitary wave propagation in a fluid conduit within a viscous matrix. Journal of Geophysical Research: Solid Earth 91 (B6), 6367–6374.
  * Scott & Stevenson (1984) Scott, David R & Stevenson, David J 1984 Magma solitons. Geophysical Research Letters 11 (11), 1161–1164.
  * Scott et al. (1986) Scott, David R, Stevenson, David J & Whitehead, John A 1986 Observations of solitary waves in a viscously deformable pipe. Nature 319 (6056), 759–761.
  * SELVAM et al. (2009) SELVAM, B., TALON, L., LESSHAFFT, L. & MEIBURG, E. 2009 Convective/absolute instability in miscible core-annular flow. part 2. numerical simulations and nonlinear global modes. Journal of Fluid Mechanics 618, 323–348.
  * Simpson et al. (2006) Simpson, Gideon, Spiegelman, Marc & Weinstein, Michael I 2006 Degenerate dispersive equations arising in the study of magma dynamics. Nonlinearity 20 (1), 21.
  * Sommerfeld (1912) Sommerfeld, Arnold 1912 Die greensche funktion der schwingungsgleichung. J.-Ber. Deutsch Math.-Verein 21, 309–353.
  * Sommerfeld (1949) Sommerfeld, Arnold 1949 Partial differential equations in physics. Academic press.
  * Stubblefield et al. (2020) Stubblefield, Aaron G, Spiegelman, Marc & Creyts, Timothy T 2020 Solitary waves in power-law deformable conduits with laminar or turbulent fluid flow. Journal of Fluid Mechanics 886.
  * Whitham (1974) Whitham, GB 1974 Linear and nonlinear waves john wiley & sons. Inc., New York .
  * Yih (1963) Yih, Chia-Shun 1963 Stability of liquid flow down an inclined plane. The physics of Fluids 6 (3), 321–334.
  * Yih (1967) Yih, Chia-Shun 1967 Instability due to viscosity stratification. Journal of Fluid Mechanics 27 (2), 337–352.