Solitary wave fission of a large disturbance in a viscous fluid conduit

TL;DR

This study combines theory, numerical simulations, and experiments to analyze how large disturbances in viscous fluid conduits break into sequences of solitary waves, extending soliton resolution concepts to non-integrable systems.

Contribution

It introduces the solitary wave resolution method for non-integrable systems and validates its predictions through experiments and simulations in viscous fluid conduits.

Findings

Predicted number of solitary waves matches experimental observations within 1-2 waves.

Amplitude distribution of waves shows strong agreement with theoretical predictions.

Universal properties of solitary wave fission are identified across fluid dynamics problems.

Abstract

This paper presents a theoretical and experimental study of the long-standing fluid mechanics problem involving the temporal resolution of a large, localised initial disturbance into a sequence of solitary waves. This problem is of fundamental importance in a range of applications including tsunami and internal ocean wave modelling. This study is performed in the context of the viscous fluid conduit system-the driven, cylindrical, free interface between two miscible Stokes fluids with high viscosity contrast. Due to buoyancy induced nonlinear self-steepening balanced by stress induced interfacial dispersion, the disturbance evolves into a slowly modulated wavetrain and further, into a sequence of solitary waves. An extension of Whitham modulation theory, termed the solitary wave resolution method, is used to resolve the fission of an initial disturbance into solitary waves. The developed theory predicts the relationship between the initial disturbance's profile, the number of emergent solitary waves, and their amplitude distribution, quantifying an extension of the well-known soliton resolution conjecture from integrable systems to non-integrable systems that often provide a more accurate modelling of physical systems. The theoretical predictions for the fluid conduit system are confirmed both numerically and experimentally. The number of observed solitary waves is consistently within 1-2 waves of the prediction, and the amplitude distribution shows remarkable agreement. Universal properties of solitary wave fission in other fluid dynamics problems are identified.

Tables

Table 1. Table 1: Densities, viscosities, viscosity ratio, background flow rate Q 0 subscript 𝑄 0 Q_{0} , and associated background conduit diameter 2 ​ R 0 2 subscript 𝑅 0 2R_{0} , mean flow rate U 0 subscript 𝑈 0 U_{0} according to equations ( 8 ), ( 9 ), respectively for the reported experiments (except figure  2 (b)).

${\mu }^{\left(i\right)}$	$3.66\text{×}{10}^{-2}\mathrm{Pa}s$
${\mu }^{\left(e\right)}$	$1.296\mathrm{Pa}s$
${\rho }^{\left(i\right)}$	$1.198g{\mathrm{cm}}^{-3}$
${\rho }^{\left(e\right)}$	$1.260g{\mathrm{cm}}^{-3}$
$\epsilon$	0.0283
$Q_0$	$0.50{\mathrm{cm}}^3{\min }^{-1}$
$2R_0$	$2.1\mathrm{mm}$
$U_0$	$2.3\mathrm{mm}{s}^{-1}$

Table 2. Table 2: Integrals of conserved and nonconserved quantities for the solution depicted in figure  8 at t = 350 𝑡 350 t=350 over three spatial intervals corresponding to the whole domain [ 0 , L ] 0 𝐿 [0,L] , the subinterval containing only small amplitude, dispersive waves [ 0 , z ∗ ] 0 subscript 𝑧 [0,z_{*}] , and the solitary wavetrain subinterval [ z ∗ , L ] subscript 𝑧 𝐿 [z_{*},L] ( z ∗ = 1000 subscript 𝑧 1000 z_{*}=1000 , L = 1500 𝐿 1500 L=1500 ).

$I=$	$[0,L]$	$[0,z_{\ast }]$	$[z_{\ast },L]$
${\int }_I\left(a-1\right)dz$	84.48	-0.34	84.82
${\int }_I\left(1-\frac{1}{a}-\frac{a_z^2}{a^2}\right)dz$	46.36	-0.35	46.71
${\int }_I{\left(a-1\right)}^{2}dz$	81.39	0.01	81.38

Equations

$$a_t+(a^2{)}_z-{\left(a^2{\left(a^{-1}{a}_t\right)}_z\right)}_z=0.$$

$$\begin{array}[]{cc}a(z,0)=1+a_{0}(z),&\lim_{\left|z\right|\to\infty}a_{0}(z)=0,\end{array}$$

$$a_m=\underset{z\in \mathbb{R}}{\max }{a}_0(z).$$

$$u_{\tau }+u{u}_x+u_{xxx}=0$$

$$\begin{array}[]{ccc}\tau=\delta^{3/2}t,&x=\delta^{1/2}2^{-1/3}\left(z-2t\right),&u=2^{2/3}\delta(a-1),\quad 0<\delta\ll 1,\end{array}$$

$$\begin{array}{rl}\mathcal{N}& =\frac{1}{\pi \sqrt{6}}{\int }_{-\infty }^{\infty }\sqrt{u_0(x)}\mathrm{d}x,\\ f(\mathcal{A})& =\frac{1}{4\pi \sqrt{6}}{\int }_{x_1}^{x_2}\frac{\mathrm{d}x}{\sqrt{u_0(x)-\mathcal{A}/2}},0\le \mathcal{A}\le 2u_m.\end{array}$$

$$\begin{array}{rl}\mathrm{for}{u}_0(x)& =\{\begin{array}{ll}{u}_m& -w<x<0\\ 0& \mathrm{else}\end{array},\mathcal{N}=\frac{w\sqrt{u_m}}{\pi \sqrt{6}},\\ \mathcal{F}(\mathcal{A})& =\frac{1}{\mathcal{N}}{\int }_0^{\mathcal{A}}f({\mathcal{A}}^{\prime })\mathrm{d}{\mathcal{A}}^{\prime }=1-\sqrt{1-\frac{\mathcal{A}}{2u_m}},0\le \mathcal{A}\le 2u_m,\end{array}$$

$$2R_0={\left(\frac{2^7{\mu }^{(i)}{Q}_0}{\pi g({\rho }^{(e)}-{\rho }^{(i)})}\right)}^{1/4},$$

$$U_0=\frac{g{R}_0^2({\rho }^{(e)}-{\rho }^{(i)})}{8{\mu }^{(i)}}.$$

$$c(a_s)=\frac{2a_s^2\log a_s-a_s^2+1}{(a_s-1{)}^2},$$

$$\mathcal{F}(\mathcal{A})=\frac{\text{number of solitary waves with amplitudes}{a}_s\text{satisfying}{a}_s\le 1+\mathcal{A}}{\text{total number of solitary waves}}.$$

$$A_T+\frac{g\Delta }{8\pi {\mu }^{(i)}}{\left(A^2\right)}_Z-\frac{{\mu }^{(e)}}{8\pi {\mu }^{(i)}}{\left(A^2{\left(A^{-1}{A}_T\right)}_Z\right)}_Z=0,$$

$$a=\frac{1}{\pi R_0^2}A,z=\frac{\sqrt{8\epsilon }}{R_0}Z,t=\frac{\sqrt{8\epsilon }{U}_0}{R_0}T.$$

$$c_s(a_s,\overline{\varphi })=\frac{\overline{\varphi }\left(2a_s^2\left(\log a_s-\log \overline{\varphi }\right)-a_s^2+{\overline{\varphi }}^2\right)}{(a_s-\overline{\varphi }{)}^2}.$$

$$\begin{array}[]{ccc}a(z,t)=\phi(\theta),&\theta=kz-\omega t,&\phi(\theta+2\pi)=\phi(\theta).\end{array}$$

$$({\varphi }^{\prime }{)}^2=g(\varphi )=-\frac{2}{k^2}\varphi -\frac{2}{\omega k}{\varphi }^2\log \varphi +C_0+C_1{\varphi }^2,$$

$${\omega }_0(k,\overline{\varphi })=\frac{2k\overline{\varphi }}{1+k^2\overline{\varphi }}.$$

$$\{\begin{array}{l}{a}_t+(a^2-a^2(a^{-1}{a}_t{)}_z{)}_z=0,\\ {\left(\frac{1}{a}+\frac{a_z^2}{a^2}\right)}_t+{\left(\frac{a_{tz}}{a}-\frac{a_z{a}_t}{a^2}-2\ln a\right)}_z=0\end{array}$$

$$\overline{F}=\frac{1}{2\pi }{\int }_0^{2\pi }F(\theta )d\theta ,$$

$$\{\begin{array}{l}{\overline{\varphi }}_t+{\left(\overline{{\varphi }^2}-2k\omega \overline{{\varphi }_{\theta }^2}\right)}_z=0\\ {\left(\overline{\frac{1}{\varphi }}+k^2\overline{\frac{{\varphi }_{\theta }^2}{{\varphi }^2}}\right)}_t-2{\left(\overline{\ln \varphi }\right)}_z=0\\ k_t+{\omega }_z=0\end{array},$$

$$s_{+}=\sqrt{1+8{\overline{\varphi }}_{-}}-1,s_{-}=3+3{\overline{\varphi }}_{-}-3\sqrt{{\overline{\varphi }}_{-}(8+{\overline{\varphi }}_{-})}.$$

$$c_s(a_{+},1)=\sqrt{1+8{\overline{\varphi }}_{-}}-1.$$

$$k_{-}^2=\frac{1}{4}\left(1-\frac{4}{{\overline{\varphi }}_{-}}+\sqrt{\frac{1}{{\overline{\varphi }}_{-}}(8+{\overline{\varphi }}_{-})}\right).$$

$$\mathcal{N}=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }k(z,t)\mathrm{d}z.$$

$${\beta }_t+2\beta {\beta }_z=0.$$

$$\begin{array}[]{lll}z=z_{-}(t):&\mathcal{A}=0,&\overline{\phi}=\beta_{-}(t),\\ z=z_{+}(t):&k=0,&\overline{\phi}=\beta_{+}(t)=1,\end{array}$$

$$\mathcal{A}=0:$$

$${\overline{\varphi }}_t+2\overline{\varphi }{\overline{\varphi }}_z=0,$$

$$k_t+{\left({\omega }_0(k,\overline{\varphi })\right)}_z=0.$$

$$\frac{dz}{dt}=2\overline{\varphi }\mathrm{and}\frac{dz}{dt}={\omega }_{0,k}.$$

$$\frac{d{k}_{-}}{d\overline{\varphi }}=\frac{{\omega }_{0,\overline{\varphi }}}{2\overline{\varphi }-{\omega }_{0,k_{-}}}.$$

Full text

1 Introduction

A fundamental problem in fluid dynamics is the long-time resolution of a large, localised disturbance. In inviscid fluids, a prominent feature of this resolution is the emergence of a solitary wavetrain. This process is generally referred to as soliton fission and has been observed in a variety of fluid contexts. For example, while intense earthquakes can lead to the vertical displacement of the ocean surface by several meters, its horizontal extent can reach 10–100 kilometers Geist et al. (2007), which, under appropriate shallowness conditions, can evolve into a large number of surface solitary waves (Matsuyama et al., 2007; Arcas & Segur, 2012). Another important example is the generation of large amplitude internal ocean solitary waves with two identified soliton fission mechanisms: 1) an initial, broad displacement of internal temperature and salinity (Osborne & Burch, 1980) and 2) the propagation of a large internal solitary wave onto a shelf (Farmer & Armi, 1999; Vlasenko et al., 2014). In both scenarios, the result is the same—the generation of a large number of rank-ordered solitary waves. In fact, the well-known soliton fission law by Djordjevic & Redekopp (1978) for scenario 2 was obtained by modeling it with an initial, broad disturbance to the constant coefficient Korteweg-de Vries (KdV) equation, a weakly nonlinear, long wave model. More generally, the disintegration of a broad disturbance into solitary waves is the inevitable result of boundary or topography interaction with an undular bore or dispersive shock wave (DSW) that results from a sharp gradient due to a variety of reasons El & Hoefer (2016).

Despite the prevalence of soliton fission in fluid dynamics, its theoretical description has primarily been limited to completely integrable partial differential equations (PDEs) such as the KdV equation. First attempts to understand this problem began with the celebrated Zabusky-Kruskal numerical experiment of an initial cosine profile for KdV (Zabusky & Kruskal, 1965). Asymptotics of KdV conservation laws (Karpman, 1967; Johnson, 1973) and the inverse scattering transform (Segur, 1973; Deng et al., 2016) yield a prediction for the number of solitons based on eigenvalue counting and an estimate for the amplitudes of fissioned solitons from an initial profile.

Because of its ubiquity, we seek a deeper understanding of soliton fission that results from a broad initial condition, hereafter referred to as the box problem due to the initial profile’s wide shape. A new method based on Whitham averaging theory (Whitham, 1974) that does not require integrability was first proposed and applied to the Serre/Su-Gardner/Green-Naghdi equations for fully nonlinear shallow water waves in El et al. (2008) and, partially, to the defocusing nonlinear Schrödinger equation with saturable nonlinearity in El et al. (2007). The method draws upon principles first developed to describe DSWs that result from step initial data (El, 2005). The long-time evolution into a solitary wavetrain is one component of the soliton resolution conjecture, which proposes that localised initial conditions to nonlinear dispersive wave equations generically evolve into a soliton wavetrain and small amplitude dispersive radiation, originally formulated within the context of integrable PDEs such as the KdV equation (Segur, 1973; Schuur, 1986; Deift et al., 1994). Because the method presented here is not reliant on integrability of the underlying PDE and yields concrete predictions for the number of solitary waves and their amplitude distribution that result from broad initial disturbances, we refer to this approach as the solitary wave resolution method. An important feature of the solitary wave resolution method is that it bypasses an analysis of the full Whitham modulation equations—which are generally difficult to analyse—in favor of the exact zero amplitude and zero wavenumber reductions of the full Whitham equations that admit a general structure and form that is amenable to further analysis.

We note that the solitary wave resolution method does not resolve the second component of the soliton resolution conjecture—the small amplitude dispersive radiation. For sufficiently broad boxes, this component of the conjecture is negligible, as is well-known for the KdV equation (see, e.g., Karpman (1974); Whitham (1974)). We quantify the contributions of the radiation and solitary wave components numerically for a specific initial disturbance in the conduit equation.

Our basic hypothesis in this work is that a large initial disturbance in certain nonintegrable equations (e.g., the conduit equation described below) results most prominently in the fission of solitary waves, which enables us to apply Whitham averaging theory (Whitham, 1974). This hypothesis is motivated by rigorous semiclassical analysis of the KdV equation (Lax & Levermore, 1979) and is confirmed by numerical simulations for the conduit equation. Moreover, this hypothesis has been successfully applied to the nonintegrable Serre equations (El et al., 2008).

The solitary wave fission problem has been studied experimentally, primarily in water wave tanks modeled by the KdV equation (Hammack & Segur, 1974, 1978). More recent water wave experiments physically recreated the Zabusky-Kruskal numerical experiment, observing recurrence as well as soliton fission (Trillo et al., 2016). These experiments exhibit excellent agreement with WKB theory applied to the inverse scattering transform (Deng et al., 2016). The WKB approach has also been applied to the defocusing nonlinear Schrödinger equation, yielding the number and amplitudes of emergent solitons (Deng et al., 2017). However, none of these quantitative methods are applicable to non-integrable equations.

This paper presents solitary wave fission experiments and modulation theory for the interfacial dynamics between two high-viscosity, miscible fluids, one rising buoyantly within another. Original experiments demonstrated that solitary waves preserve their shape and form despite long distance propagation and interaction with other solitary waves (Olson & Christensen, 1986; Scott et al., 1986). In fact, in both of these experimental papers, solitary waves were generated by the fission of a large, initial disturbance. We apply modulation theory for solitary wave fission introduced in (El et al., 2008) to the box problem for a PDE model of this fluid context known as the conduit equation Lowman & Hoefer (2013a)

[TABLE]

In the derivation of the conduit equation, no restriction is placed on the magnitude of the nondimensional, circular cross-sectional area $a(z,t)$, where $z$, $t$ are the scaled height and time, respectively, assumed to be much larger than the conduit diameter and a characteristic advective time scale. This equation is also an asymptotic, long-wave model of magma flows rising through the Earth’s mantle (Barcilon & Richter, 1986; Whitehead & Helfrich, 1988; Helfrich & Whitehead, 1990) and to a comparatively simple laboratory experiment (Olson & Christensen, 1986; Scott et al., 1986; Whitehead & Helfrich, 1988; Helfrich & Whitehead, 1990; Lowman & Hoefer, 2013a; Maiden et al., 2016, 2018; Anderson et al., 2019). Equation (1) fails the so-called Painlevé test for integrability (Harris & Clarkson, 2006) and has at least two conservation laws (Harris, 1996) therefore is an excellent candidate to test the more broadly applicable solitary wave resolution method for the initial value problem consisting of (1) and

[TABLE]

where $a_{0}(z)$ is a broad, localised disturbance with exactly one critical point at the maximum

[TABLE]

Note that $a=a_{m}+1$ at the maximum, i.e., $a_{m}$ measures the amplitude of the disturbance exceeding the unit background area ratio $a=1$. We will quantify the profile’s broadness more precisely later on but for $a_{0}(z)$ in the shape of a box, then a sufficiently wide box will do. Formally, $a_{0}\in\mathcal{C}^{\infty}(\mathbb{R})$ as well, but the relaxation of this assumption still aligns with the theoretical results. We shall assume that the support of $a_{0}(z)$ is $[-w,0]$ where $w>0$ is the box width. The solitary wave resolution method utilizes the characteristics of the Whitham modulation equations to estimate the number of solitary waves and the solitary wave amplitude distribution resulting from a large-scale initial condition. An example initial condition and its numerically evolved state according to the conduit equation (1) are shown in figure 1. The theoretically predicted solitary wave number (12, derived in section 4) is correct and the predicted amplitudes fall well within the ranges determined from the quantization of the continuous amplitude distribution (derived in section 5).

The conduit equation (1) can be approximated by the KdV equation

[TABLE]

in the small-amplitude, long-wavelength regime with the scaling (Whitehead & Helfrich, 1986)

[TABLE]

where $\delta$ is a characteristic disturbance amplitude deviation from unit background. The formulae for the expected number of solitons $\mathcal{N}$ and the amplitude ($\mathcal{A}$) density function $f(\mathcal{A})$ for the initial profile $u(x,0)=u_{0}(x)$ to the KdV equation (4) are (El et al., 2008)

[TABLE]

Here, $x_{1}$ and $x_{2}$ are the intersections of the initial condition $u_{0}$ with the value $\mathcal{A}/2$. The initial data is assumed to be on a zero background with maximum $u_{m}=\max u_{0}(x)$. For initial data consisting of a box of width $w$ and height $u_{m}$, equation (6) becomes

[TABLE]

where $\mathcal{F}(\mathcal{A})$ is the cumulative distribution function of the soliton amplitudes normalised by the total number of solitons. The number of solitons agrees with that obtained by IST-related approaches (Karpman, 1967; Ablowitz et al., 2009) in its asymptotic regime of validity $\mathcal{N}\sim w\sqrt{u_{m}}\gg 1$. This modulation theory approach to solitary wave fission can be applied to any dispersive nonlinear wave equation that admits a Whitham modulation description (Whitham, 1974; El & Hoefer, 2016). We identify certain universal properties of solitary wave fission, including the independence of the normalised cumulative distribution function on box width (e.g., $\mathcal{F}(\mathcal{A})$ is independent of $w$). We also predict the linear dependence on box width $w$ of the solitary wave fission number $\mathcal{N}$.

This paper continues with section 2 where we present viscous fluid conduit fission experiments. Section 3 includes relevant background information on the conduit equation (1). In sections 4 and 5, we develop the solitary wave resolution method to estimate the number of solitary waves and their amplitude distribution. In light of the developed modulation theory, we return to the experiments in section 6. We wrap up with concluding remarks in section 7.

2 Observation of solitary wave fission

We motivate our analysis by first presenting viscous fluid conduit experiments on solitary wave fission.

2.1 Experimental setup

The experimental setup is nearly identical to that used by Anderson et al. (2019) and consists of a square acrylic column with dimensions 4\text{\,}\mathrm{cm}$\times$4\text{\,}\mathrm{cm}$\times$200\text{\,}\mathrm{cm}, filled with glycerine, as shown in figure 2(a). The interior fluid (identified by the superscript (i)) consists of certain ratio of glycerine, water, and black food coloring, which is injected through a nozzle installed at the column’s base. The ratio is chosen so that the interior fluid has both lower density, $\rho^{(i)}<\rho^{(e)}$, and significantly lower viscosity, $\mu^{(i)}\ll\mu^{(e)}$, than the exterior fluid denoted by the superscript (e). Miscibility of the two fluids implies that surface tension effects are negligible. The nominal parameter values used in the experiments presented here are those in table 1.

A high precision, computer-controlled piston pump is used to inject the interior fluid with a pre-determined temporal flow profile. Buoyancy and steady injection at a fixed volumetric flow rate ($Q_{0}$ in table 1) leads to a vertically uniform fluid conduit, which is referred to as the background conduit, and is verified to be well-approximated by the pipe (Poiseuille) flow relation (see, e.g., the supplementary material in Maiden et al. (2016))

[TABLE]

where $2R_{0}$ is the conduit diameter and $g$ is the acceleration due to gravity. The mean vertical advective velocity within the conduit according to pipe flow is

[TABLE]

Data acquisition is performed using high resolution digital cameras equipped with macro lenses, one to capture the initial box profile and one near the top of the apparatus (the far-field) to capture the solitary wavetrain. A ruler is positioned beside the column within camera view for calibration purposes in order to quantify the observed box width. All amplitudes (solitary wave and box) are reported as cross-sectional areas that are normalised to the observed mean background cross-sectional area. This facilitates future comparison with theoretical results for the nondimensional conduit equation (1).

2.2 Methods

We use the characteristic control method described in (Anderson et al., 2019) to generate a volumetric flow rate profile that results in a box-like structure in the lower part of the column with a pre-specified width $w$ and nondimensional cross-sectional area $a_{m}+1$. Figure 2(a) displays a schematic of the experiment and figure 2(b) depicts the experimental time development of a box-like profile. The lower “box camera” takes several images before, during, and after the predicted box development time. After the leading edge of the box forms, the pump rate is quickly reduced to the background rate $Q_{0}$, and the box evolves into oscillations that rise up the conduit. Once the leading oscillation reaches the upper “wave camera” imaging window, images are taken at $0.2\text{\,}\mathrm{Hz}$ for several minutes, to ensure that all waves originating from the box have had sufficient time to propagate through the viewing window. For the experiment reported only in figure 2(b), an additional camera (not shown in figure 2(a)) is employed to image a $1.25\text{\,}\mathrm{m}$ section of the column. The large aspect ratio of the full columnar dynamics imply the relatively low image resolution of $43\text{\,}\mathrm{pixel}\text{\,}{\mathrm{cm}}^{-1}$. These dynamics will be directly compared with the evolution predicted by the conduit equation in Section 6.

The approximately white background and the opaque, black conduit yield sufficient contrast for edge detection by identifying the two midpoints between the maximum and minimum of the spline interpolated horizontal image intensity. The edge data is then processed with a low-pass filter to reduce pixelation noise and the effects of impurities in the exterior fluid. The number of pixels between the two edges is identified as the conduit diameter, which is squared and normalised by the squared observed background conduit diameter to obtain the dimensionless cross-sectional area $a$. Our imaging setup at both the “box camera” and “wave camera” in figure 2(a) admit resolutions of 300 pixel cm*-1* and between 132–228 pixel cm*-1*, respectively (generally higher resolution for smaller boxes).

We use the lower camera to determine the box shape. Note that near the point of breaking, dispersion is no longer negligible; as a result, a pure box is difficult to realize in the conduit system. We use the characteristic control method presented in (Anderson et al., 2019) to extract the time of box profile formation, and use the nondimensionalized version of that profile as the initial condition in further analyses of the conduit equation. An example experimental box profile is shown in the $t=0$ panel of figure 2(b).

For the upper camera, a wave-tracking algorithm is utilised to follow all wave peaks across the imaging window. Each candidate peak’s amplitude and position are validated against the conduit equation’s solitary wave speed-amplitude relation Olson & Christensen (1986)

[TABLE]

during the temporal window that the peak is in view. The elevation solitary wave amplitude $a_{s}>1$ is measured from zero area, hence must be larger than the background area $a=1$. Since solitary waves exhibit the speed lower bound $c(a_{s})>2$, any observed wave peak with a slower speed was discarded as small amplitude, dispersive wave phenomena.

2.3 Results

A total of 30 experimental trials were executed according to the protocol described in this Section with the experimental parameter values identified in table 1. We generated 15 box geometries and carried out 2 trials per geometry with the nominal nondimensional box heights $a_{m}\in\{1,2,3\}$ and nominal box widths $\{20,25,30,35,40\}$ cm—corresponding to nondimensional widths $w\in\{90,112,134,156,178\}$ in the conduit equation (1) (the nondimensionalisation will be provided in Section 3). The results are shown in figures 3(a-d). Figure 3(a) depicts the observed number of solitary waves as functions of box width and nondimensional height $a_{m}$. Across all trials, between 8 and 20 solitary waves were generated, placing these fission experiments in the large number of solitary waves regime, as expected for broad initial conditions. For fixed box height, the data exhibit an approximately linear increase with width, as shown by the linear fits for fixed $a_{m}$ in figure 3(a). Figures 3(b-d) report the normalised cumulative distribution functions (cdfs) of solitary wave amplitudes. Each panel (b-d) includes the normalised cdfs for trials with a common box height value $a_{m}$. The normalised cdf $\mathcal{F}(\mathcal{A})$ depends parametrically on box width $w$ and box height $a_{m}$ and is defined as

[TABLE]

While the staircase cdfs plotted in figures 3(b-d) have different shapes for different box heights, the box width dependence of the normalised cdfs for fixed box height—i.e., for each fixed panel—shows little variation.

The rest of this paper is concerned with developing a modulation theory description for solitary wave fission in the conduit equation box problem (equations (1) and (2)). Our analysis results in explicit predictions for the number of solitary waves—linearly dependent on box width—and the normalised amplitude cdf—independent of box width. Following our analysis, we will reconsider the experiments presented here.

3 Conduit equation background

The conduit equation (1) describes the dynamics of the free interface between two viscous fluids: a highly dense, highly viscous exterior fluid, and a less dense, less viscous interior fluid. As the interior fluid is pumped steadily through the exterior fluid, the interface resembles a deformable pipe whose walls are the two-fluid boundary. The circular cross-sectional area $A$ of this pipe can be modeled as a function of time $T$ and vertical distance $Z$ by the dimensional conduit equation (Olson & Christensen, 1986; Lowman & Hoefer, 2013a)

[TABLE]

when $\varepsilon=\frac{\mu^{(i)}}{\mu^{(e)}}$, the interior to exterior dynamic viscosity ratio, is small, $\Delta=\rho^{(e)}-\rho^{(i)}$ is the difference in exterior to interior fluid densities, and $g$ is gravitational acceleration. Equation (12) results from the interplay between interior fluid buoyancy and continuity of both the fluid velocity and interfacial stress at the two fluid boundary. This nonlinear dispersive partial differential equation is a long-wave, slowly varying asymptotic reduction of the Navier-Stokes equations for two fluids. Restrictions include sufficiently small Reynolds number and small interfacial steepness but there is no restriction on the conduit amplitude, hence the dispersive term is nonlinear (Lowman & Hoefer, 2013a). The nondimensional form of (12) is (1), obtained via the scalings (c.f. (8), (9))

[TABLE]

This transformation rescales the background conduit area of radius $R_{0}$ to unity. The conduit equation has been shown to admit a variety of multiscale coherent wave solutions (Maiden & Hoefer, 2016).

Previous experimental comparisons to dynamics predicted by the conduit equation include solitary waves (Olson & Christensen, 1986), their interactions with each other (Helfrich & Whitehead, 1990; Lowman et al., 2014), and interactions with a dynamically changing mean flow (Maiden et al., 2016, 2018). Solitary wave solutions can be obtained from the ordinary differential equation that results from the travelling wave ansatz $a(z,t)=f(z-ct)$, where the solitary wave speed $c$ is related to its total amplitude $a_{s}$ (measured from $a=0$) on the background $\overline{\phi}$ by the speed-amplitude relation (Olson & Christensen, 1986)

[TABLE]

Dispersive shock waves have also been studied theoretically and experimentally in the viscous fluid conduit system (Lowman & Hoefer, 2013b; Maiden et al., 2016, 2018). Dispersive shock waves are the result of a sustained, large increase in background conduit area from 1 to $\overline{\phi}_{-}>1$ and can be characterised by a modulated periodic travelling wave solution of the conduit equation, i.e., a solution of the form

[TABLE]

Inserting this ansatz into equation (1) and integrating twice results in (Olson & Christensen, 1986)

[TABLE]

where $C_{0}$ and $C_{1}$ are real constants of integration. The right side of the equation can have up to three roots, $\phi_{1}\leq\phi_{2}\leq\phi_{3}$, which parameterise the solution.

A physically relevant parameterisation of the periodic wave $\phi(\theta)$ is given by three constants: the wavenumber $k$, the wave amplitude $\mathcal{A}$ (defined as the difference between the wave’s maximum and minimum), and the wave mean $\overline{\phi}$, which can be written in terms of $C_{0}$, $C_{1}$, and $k$, or equivalently, in terms of $\phi_{j}$, $j=1,2,3$. The wave frequency is determined by the $2\pi$ periodicity of $\phi(\theta)$ as $\omega=\omega(k,\overline{\phi},\mathcal{A})$. The modulation theory description of a DSW is achieved by allowing the periodic wave’s parameters to vary slowly relative to the wavelength $2\pi/k$ and period $2\pi/\omega$ while introducing the generalised wavenumber $k=\theta_{x}$ and frequency $\omega=-\theta_{t}$ (Lowman & Hoefer, 2013b). Then a DSW can be viewed as connecting two distinguished limits of these modulated wave parameters: the zero amplitude limit as $\mathcal{A}\to 0$ and the zero wavenumber limit $k\to 0$. When $\mathcal{A}\to 0$, the DSW solution limits to small amplitude harmonic waves with the linear dispersion relation

[TABLE]

When $k\to 0$, the DSW solution limits to a solitary wave that satisfies the speed-amplitude relation (14).

Allowing for slow modulations of $\overline{\phi}$, $k$, and $\mathcal{A}$ in space and time results in the conduit-Whitham equations. The conduit-Whitham equations consist of the conservation of waves $k_{t}+\omega_{x}=0$, resulting from $\theta_{tx}=\theta_{xt}$, and the averaging of the conduit equation’s two conservation laws (Barcilon & Richter, 1986)

[TABLE]

over the periodic wave family. Using the following notation for averaging over a wave period

[TABLE]

the Whitham equations are (Maiden & Hoefer, 2016)

[TABLE]

where $\omega=\omega(k,\overline{\phi},\mathcal{A})$ is the nonlinear wave frequency. That the averaging operator (19) approximately commutes with partial differentiation is a result of scale separation between the modulation—which is large and slow—and the periodic wave’s much shorter and faster spatial wavelength and temporal period, respectively (Whitham, 1965, 1974).

We remark that a rigorous, necessary condition for the stability of conduit periodic waves is the hyperbolicity of the conduit-Whitham equations (Johnson & Perkins, 2019). The conduit-Whitham equations are known to be hyperbolic in an amplitude/wavenumber dependent regime of phase space (Maiden & Hoefer, 2016) and we will operate within this regime.

If a certain self-similar, simple wave solution to the conduit-Whitham equations exists (a 2-wave (El & Hoefer, 2016)), we can obtain expressions for the leading (solitary wave) and trailing (harmonic) edge speeds in terms of the DSW jump parameter $\overline{\phi}_{-}$, labeled $s_{+}$ and $s_{-}$, respectively (Lowman & Hoefer, 2013b)

[TABLE]

The solitary wave amplitude $a_{+}$ is implicitly determined by equating $s_{+}$ with the solitary wave speed-amplitude relation (14)

[TABLE]

The trailing edge small amplitude wavepacket is characterized by the wavenumber $k_{-}$, explicitly determined by equating the linear group velocity $\partial_{k}\omega_{0}$ to $s_{-}$

[TABLE]

The group velocity of the harmonic edge is always less than the speed of the solitary wave edge. Thus, a DSW in the conduit system is always led by a solitary wave, with a trailing, continually expanding, oscillating wavetrain that can exhibit backflow and instabilities for sufficiently large jumps $\overline{\phi}_{-}$ (Lowman & Hoefer, 2013b; Maiden & Hoefer, 2016).

We now return to the initial value problem (2) for the conduit equation (1) and our development of the solitary wave resolution method. Initially, the edges of the wide box (2) can be treated as two well-separated, step-like (Riemann) initial value problems. As such, the rightmost edge will evolve similar to a DSW, and the leftmost edge similar to a rarefaction wave (RW). However, finite box extent necessarily implies the eventual interaction of the DSW and RW (El & Grimshaw, 2002; Ablowitz et al., 2009). Ultimately, a finite solitary wavetrain emerges from this interaction process. We now use a modification of conduit DSW theory (Lowman & Hoefer, 2013b) to determine the properties of this solitary wavetrain by applying the solitary wave resolution method originally developed in (El et al., 2008).

4 Number of solitary waves

In what follows, we make the assumption that the box initial value problem (2) for equation (1) will result in a slowly modulated wavetrain that can be described by the Whitham modulation equations (20). This assumption is the cornerstone of the solitary wave resolution method (El et al., 2008) and will later be verified by numerical simulations.

Allowed to evolve long enough, the individual wave crests resulting from the box initial conditions will separate with minimal overlap, i.e., will result in a non-interacting solitary wavetrain. To count these waves, note that they are separated by exactly their wavelength, defined in terms of the wavenumber as $2\pi/k$. Consequently, $k/2\pi$ is a wave crest density and we determine the total number of waves $\mathcal{N}$ in a wavetrain at time $t$ by

[TABLE]

This integral is finite at $t=0$ because the initial disturbance $a_{0}(z)$ has compact support, implying $k\to 0$ as $\left|z\right|\to\infty$ sufficiently fast. The conservation of waves equation in the conduit-Whitham modulation equations (20) implies that $\mathcal{N}$ is independent of time. Then the total number of fissioned solitary waves that emerge in long time can be determined by the wavenumber function $k(z,0)$ associated with the initial condition. The challenge is to determine $k(z,0)$ when the waves are initially so densely packed that there are no visible oscillations, i.e., the wave amplitude $\mathcal{A}=0$ and there is only the nonzero mean $\overline{\phi}(z,0)$.

Whitham modulation theory can now be utilised to find a relationship between the initial condition—the non-oscillatory data (2) equates to the initial mean $\overline{\phi}(z,0)=1+a_{0}(z)$ in modulation theory—and the wavenumber $k$. For this, we note that the conduit-Whitham equations (20) are supplemented by conditions that ensure continuity of the modulation solution at the trailing and leading edges of the oscillatory wavetrain for all $t$. There exist only two ways for the modulation solution to continuously match to the solution of the dispersionless conduit equation

[TABLE]

Either $k\to 0$ or $\mathcal{A}\to 0$. The case $k\to 0$ is the solitary wave limit and $\mathcal{A}\to 0$ is the small amplitude, harmonic wave limit. These limits are important for the modulation solution of a DSW, with $\mathcal{A}\to 0$ at the leftmost, trailing edge and $k\to 0$ at the rightmost, leading edge. Because early to intermediate time evolution leads to the generation of a DSW, we identify these edges as $z_{-}(t)$ and $z_{+}(t)$, respectively, and require

[TABLE]

where the wave mean $\overline{\phi}$ matches to the solution $\beta(z,t)$ of the dispersionless conduit equation (25) subject to the initial condition $\beta(z,0)=1+a_{0}(z)$ (cf. (2)). Then, $\beta_{\pm}(t)=\beta(z_{\pm}(t),t)$. Consequently, equation (25) is valid outside the oscillatory region influenced by the disturbance, i.e., for $z\in(-\infty,z_{-}(t))\cup(z_{+}(t),\infty)$. We note that the dispersionless conduit equation (25) (buoyancy driven flow with negligible curvature induced interfacial stress) has been experimentally shown to be a good approximation to the physical conduit system when there are no oscillations, i.e., when the interface is slowly varying (Anderson et al., 2019). The edges $z_{\pm}(t)$ in the boundary matching problem (20), (26) with $\beta_{-}\neq\hbox{const}$ can be determined by a recent extension of the DSW fitting method developed by Kamchatnov (2019). This determination won’t be necessary in our construction in which we seek the long-time solitary wave resolution.

When $\mathcal{A}\to 0$, the vanishing oscillations do not contribute to the averaging (19), so $\overline{F(\phi)}=F(\overline{\phi})$, for any differential or algebraic operator $F$ (El & Hoefer, 2016). Thus all $\theta$ derivatives of $\phi$ average to zero. In this case, the first and second conduit-Whitham equations (20) limit to the dispersionless conduit equation (25) but the conservation of waves modulation equation remains and the wave frequency is the linear dispersion relation (17) so that the modulation system reduces to

[TABLE]

Since the disturbance is initially non-oscillatory, we have $\overline{\phi}(z,0)=1+a_{0}(z)$, $z\in\mathbb{R}$ (c.f. 2). However, because there are no initial oscillations, the initial wavenumber is not well-defined. We must appeal to properties of the disturbance’s evolution in order to uniquely define $k(z,0)$. We do so by identifying a simple wave relationship $k=k_{-}(\overline{\phi})$ between the wavenumber and mean so that $k(z,0)=k_{-}(\overline{\phi}(z,0))$. The rationale for the use of the simple wave relation is detailed in (El et al., 2008) and is based on the fact that the DSW trailing edge is a characteristic. Equations (27) and (28) have two characteristic families

[TABLE]

The first family corresponds to the decoupled evolution of the mean flow equation (27) and coincides with the slowly varying evolution of the disturbance, e.g., the initial RW. The second family coincides with the vanishingly small amplitude oscillations emerging from the edge of the evolving disturbance with an envelope that moves with the group velocity. It is the second characteristic family that captures the evolution of the emergent solitary wavetrain. In order to obtain the relationship between $k$ and $\overline{\phi}$ along the second characteristic family, we make the simple wave ansatz $k=k_{-}(\overline{\phi})$ along $z_{-}(t)$ where $dz/dt=\omega_{0,k}$, which, when combined with the modulation equations (27), results in the ODE

[TABLE]

Substituting the linear dispersion relation (17) into this equation and integrating yields an expression for $k_{-}$ in terms of the wave mean $\overline{\phi}$ and an integration constant $\lambda$

[TABLE]

Matching to equation (25) at the disturbance’s initial termini $z\in\{-w,0\}$ via the first of equation (26), $k_{-}(\overline{\phi}=1;\lambda)=0$. Then $\lambda=1/2$. This choice of integration parameter results in the same expression for $k_{-}$ as found at the DSW’s harmonic edge from DSW fitting theory; see equation (23). It will be useful in the next section to use the fact that the simple wave curve $k=k_{-}(\overline{\phi};1/2)$ corresponds to the level curve $\lambda=1/2$ of the surface

[TABLE]

obtained by inverting the relationship in equation (31).

The simple wave relationship $k=k_{-}(\overline{\phi};1/2)$ in equation (31) provides the needed translation between the initial condition for the mean $\overline{\phi}(z,0)=1+a_{0}(z)$ and the initial condition for the wavenumber $k(z,0)=k_{-}(1+a_{0}(z);1/2)$. Then the number of solitary waves is obtained from equation (24) as

[TABLE]

For the case when $a_{0}(z)$ is a box of width $w$ and height $a_{m}$ above a background of $1$, equation (33) can be integrated exactly

[TABLE]

The small $a_{m}$ expansion of equation (34) is

[TABLE]

which agrees to leading order with the small-amplitude KdV result in equation (7) when we identify $u_{m}=a_{m}$. The large $a_{m}$ approximation, on the other hand, is independent of box height to a good approximation

[TABLE]

To compute the number of solitary waves for more general initial profiles $a_{0}(z)$, equation (33) can be integrated. Of course, the number of solitary waves should be an integer whereas $\mathcal{N}$ continuously depends on the initial profile $a_{0}(z)$. The result is asymptotic, i.e., equation (33) is asymptotic to the number of solitary waves due to solitary wave fission if $\mathcal{N}\gg 1$. Then, computing the ceiling, floor, or rounding $\mathcal{N}$ to the nearest integer are all asymptotically equivalent. If we approximate the initial disturbance by a box of width $w$ and height $a_{m}$, equation (34) gives an explicit determination of when modulation theory for solitary wave fission is valid, i.e., when the initial disturbance is sufficiently wide.

5 Distribution of solitary wave amplitudes

Next, we seek an estimate for the amplitudes of the fissioned solitary wavetrain. Because the conduit solitary wave speed-amplitude relation (10) is monotonically increasing with amplitude, sufficiently long evolution is expected to lead to the waves separating into an amplitude-ordered train that are well isolated from one another. We will treat them as a non-interacting solitary wavetrain, a concept that was recently exploited in (Maiden et al., 2018) to describe solitary wave interaction with a mean flow. Here, we analyse both the $\mathcal{A}\to 0$ (harmonic) and $k\to 0$ (solitary wave) limits and identify a relationship between them. This enables a mapping of the initial profile to the long-time solitary wave amplitude distribution.

5.1 Harmonic limit

In the harmonic limit, $\mathcal{A}\to 0$ and Equations (27) and (28) hold. To compute the total number of solitary waves in the previous section, we identified the edges of the initial disturbance’s support and set $k=0$ at the edges. This calculation resulted in the simple wave relationship determined by the level curve $\lambda(k,\overline{\phi})=1/2$ (cf. equation (32)). We now extend this to the interior of the initial disturbance’s support and study other level curves, $\lambda(k,\overline{\phi})=$ constant, to identify the number of solitary waves contained within a portion of the initial disturbance. By solving for $\lambda$ when $k=k_{-}(1+a_{0}(z);\lambda)=0$, we therefore consider the level curves $\lambda(k,\overline{\phi})=[2(1+a_{0}(z))]^{-1}\in[1/(2(1+a_{m})),1/2]$.

We can now extend the calculation of the total number of solitary waves $\mathcal{N}$ to the number of solitary waves that emerge from the section of the initial profile of total amplitude of at least $\overline{\phi}_{\mathrm{min}}$. We modify equation (33) for the total number of fissioned solitary waves to integrate only over the initial profile section in which $1+a_{0}(z)\geq\overline{\phi}_{\mathrm{min}}$ (see figure 4(a)) and consider the $\lambda$-level curve $\lambda(k,\overline{\phi})=1/(2\overline{\phi}_{\mathrm{min}})$, determined by the zero wavenumber condition $k_{-}(\overline{\phi}_{\mathrm{min}};\lambda)=0$. Then the number of solitary waves for this truncated portion of the initial profile is

[TABLE]

The justification for this calculation comes from hyperbolicity of the modulation system (20) in the requisite domain of dependent variables (Maiden & Hoefer, 2016) and the fact that asymptotically, as $t\to\infty$, the region of influence of the support of the $\lambda$ section of the initial profile is confined by the modulation characteristics emanating from $z_{2}$ and the maximum point $z_{m}$: $a_{0}(z_{m})=a_{m}$—see El et al. (2008).

The goal now is to relate $\overline{\phi}_{\mathrm{min}}$ to the solitary wave amplitude $\mathcal{A}$. We will use the intermediate variable $\lambda$ to relate the two. With a slight abuse of notation, we define $z_{1}(\lambda)$, $z_{2}(\lambda)$ as the $z$-values at which $\overline{\phi}=1/(2\lambda)$ and $G(\lambda)$ as the number of solitary waves that emerge from the $\lambda$-section of the initial profile of total amplitude at least $1/(2\lambda)$

[TABLE]

Since $G(\lambda)$ is an increasing function of $\lambda$ and its maximum is the total number of solitary waves $\mathcal{N}=G(1/2)$ (cf. equation (33)), we define the normalised cumulative density function (cdf) $\mathcal{G}(\lambda)$ as

[TABLE]

Reconsidering the smoothed box and its numerical evolution from figure 1 in this way, profile sections for different values of $\lambda$ are shown in figure 4(a) and the expected solitary waves from each truncation are shown in figure 4(b).

The integral endpoints $z_{1}$ and $z_{2}$ for the initial condition shown in figure 4(a) are shown as a function of $\lambda$ in figure 5. The endpoints depend monotonically on $\lambda$.

5.2 Solitary wave limit

So far, we have been focused on the number of fissioned solitary waves emerging from a $\lambda$ section of the initial profile. We need to relate $\lambda$ to the amplitudes of the fissioned solitary waves. For this, we now perform an analysis of the solitary wave limit $k\to 0$ of the conduit-Whitham equations (20), which describe the modulations of a non-interacting solitary wavetrain. By use of a clever change of variables (El, 2005; El & Hoefer, 2016), this limit can be put in a form that is analogous to the harmonic limit analysis of equation (27). In particular, we will determine a relationship between the wave amplitude and mean, $\mathcal{A}=\mathcal{A}(\overline{\phi})$, that is valid along the characteristic family associated with the propagation of non-interacting solitary waves. This relationship will be a simple wave curve.

We now consider the conduit-Whitham equations (20) in the solitary wave limit $k\to 0$. Here, the wavelength $2\pi/k$ tends to infinity, so again the contribution of oscillations is negligible and averaging commutes $\overline{F(\phi)}=F(\overline{\phi})$ (El & Hoefer, 2016). Then the modulation equations reduce to the dispersionless mean flow equation and an equation for the solitary wave amplitude $\mathcal{A}$ (Maiden et al., 2018)

[TABLE]

where $c_{s}$ is the solitary wave speed-amplitude relation (14) and $g$ is a coupling function that we will not need to explicitly determine. We now introduce the convenient change of modulation variables $(\overline{\phi},\mathcal{A},k)\to(\overline{\phi},\tilde{k},\Lambda)$ (El, 2005)

[TABLE]

where $\phi_{1,2}$ are the two smaller roots of the right side of the periodic wave ODE (16).

This change of variable is based on the idea of a conjugate conduit equation, where $\tilde{a}(\tilde{z},\tilde{t})=a(i\tilde{z},i\tilde{t})$ is substituted into the conduit equation (1) so that it becomes

[TABLE]

The parameter $\tilde{k}$ is the wavenumber of the conjugate travelling wave that satisfies the ODE

[TABLE]

with the conjugate linear dispersion relation

[TABLE]

We require that periodic solutions $\phi(\theta)$ and $\tilde{\phi}(\tilde{\theta})$ to equations (16) and (44), respectively, have identical phase velocities $c_{p}=\omega/k=\tilde{\omega}/\tilde{k}$, thus $\omega=\Lambda\tilde{\omega}$. The benefit of this formulation is that the solitary wave limit of the conduit equation periodic wave is the harmonic limit of the conjugate conduit equation periodic wave, and can be leveraged as such. It allows for a formulation of the solitary wave limit that is symmetric to the harmonic limit. By substituting $k=\Lambda\tilde{k}$, $\omega=\Lambda\tilde{\omega}$ into the equation for conservation of waves $k_{t}+\omega_{z}=0$, we obtain.

[TABLE]

In the solitary wave limit, $k\to 0$ and therefore $\Lambda\to 0$, but this limit is a singular one in that $|k_{x}|,|k_{t}|\to\infty$ for a DSW and therefore $|\Lambda_{t}|,|\Lambda_{x}|\to\infty$ (El, 2005). We therefore consider equation (46) when $|\Lambda|\ll|\Lambda_{t}|,|\Lambda_{x}|$ to obtain the leading order equation

[TABLE]

This equation admits the characteristics

[TABLE]

The specific characteristic in which $\Lambda=0$ corresponds to $k=0$ and is the solitary wave edge of the wavetrain. Along $\Lambda=0$, the characteristic speed $c_{p}$ is equal to the solitary wave speed-amplitude relation $c_{s}$ (14). This relation $c_{p}=\tilde{\omega}(\tilde{k},\overline{\phi})/\tilde{k}=c_{s}(\overline{\phi}+\mathcal{A},\overline{\phi})$ determines the change of variable $(\overline{\phi},\mathcal{A})\to(\overline{\phi},\tilde{k})$ when $\Lambda=0$ for a non-interacting solitary wavetrain

[TABLE]

One can see using (14) that $\mathcal{A}\to 0$ implies $\tilde{k}\to 0$ and vice versa so $\tilde{k}$ is an amplitude type variable (El, 2005; Lowman & Hoefer, 2013a).

The next order equation when $|\Lambda|\ll|\Lambda_{t}|,|\Lambda_{x}|$ is

[TABLE]

Similar to $k_{-}$ for harmonic waves, the simple wave assumption $\tilde{k}=\tilde{k}_{+}(\overline{\phi})$ results in the ordinary differential equation

[TABLE]

whose integration results in

[TABLE]

where $\tilde{\lambda}$ is the integration constant.

5.3 Combined solitary wave and harmonic limits

Combining the simple wave results for both the harmonic wave limit $\mathcal{A}\to 0$ and the solitary wave limit $k\to 0$, we have the following characteristic integrals

[TABLE]

Compatibility between the harmonic and the solitary wave regimes within a single structure—the DSW—implies a relation between the integration constants $\lambda$ and $\tilde{\lambda}$. Indeed, if $\overline{\phi}$ is such that $k_{-}(\overline{\phi};\lambda)=0$ in $I_{H}$ then, simultaneously, $\tilde{k}_{+}(\overline{\phi};\tilde{\lambda})=0$ in $I_{S}$ (see (El et al., 2008) for details). By eliminating $\overline{\phi}$, we obtain

[TABLE]

We remark that this same result—equivalence of the integral curve parameters $\lambda=\tilde{\lambda}$ for the harmonic and solitary wave reductions—was obtained for both the KdV and Serre equations in El et al. (2008).

In long time, the solitary waves are travelling on a unit background, so inserting $\overline{\phi}=1$ and $\tilde{\lambda}=\lambda$ into equation (52) relates $\lambda$ to $\tilde{k}$

[TABLE]

Then equations (49) and (56) together identify the desired relationship between $\lambda$ and $\mathcal{A}$, the solitary wave amplitude measured from unit background $\overline{\phi}=1$

[TABLE]

Since $\lambda\in[1/(2(1+a_{m})),1/2]$, $\mathcal{A}$ is limited to the values $[0,\mathcal{A}_{\rm max}]$, where $\lambda(0)=1/2$ and $\mathcal{A}_{\rm max}$ is defined such that

[TABLE]

thus $\lambda$ is a decreasing function of $\mathcal{A}$. Using (57), we obtain the implicit expression for $\mathcal{A}_{\rm max}$

[TABLE]

This equation for the total amplitude $1+\mathcal{A}_{\rm max}$ is the same expression one obtains for the DSW’s leading edge solitary wave amplitude $a_{+}$ in equation (22) that results from an initial jump of height $a_{m}$. This concurs with our interpretation of the initial box evolution as the generation of a DSW on the right and a RW on the left. Moreover, being entirely determined by the box height, the lead solitary wave’s amplitude is predicted to be independent of box width.

Then $G(\lambda)$ from equation (38) can be written in terms of $\mathcal{A}$

[TABLE]

Because $\lambda(\mathcal{A})$ is a decreasing function of $\mathcal{A}$ and $G(\lambda)$ is an increasing function of $\lambda$, the normalised cdf of the fissioned solitary wave amplitude distribution is

[TABLE]

Since this distribution is continuous and we have a fixed number of solitary waves, we will use the quantiled discretization of this distribution for comparison with experiment and numerics.

We now attempt to explain what is seen in figure 6(b), namely that initial conditions of differing widths but otherwise the same height have approximately the same normalised cdf. To do so, we approximate the initial condition with a box of width $w$ and height $a_{m}$. Thus the normalised cdf in $\lambda$ is

[TABLE]

Since there is no variation in $z$, the numerator can be trivially integrated

[TABLE]

Then, inserting $\mathcal{N}$ from equation (34) leads to no $w$-dependence in the normalised cdf

[TABLE]

This approximation is valid as long as the edges of the disturbance transition over a small $z$ relative to $w$.

Then the normalised cdf of the amplitude distribution is obtained from (64) by substituting the functional relationship $\lambda(\mathcal{A})$ from equation (57) and noting the reflection (61)

[TABLE]

An asymptotic expansion of $\mathcal{F}$ in (65) for small $\mathcal{A}$ and $a_{m}$ yields

[TABLE]

which agrees with the weakly nonlinear KdV result (7).

5.4 Summary of the solitary wave resolution method

The above derivation is readily generalised. Consider the initial value problem for a general dispersive hydrodynamic equation

[TABLE]

with integro-differential operator $D$ yielding the real valued, linear dispersion relation $\omega_{0}(k,\overline{u})$ with negative dispersion $\omega_{0,kk}<0$. Let equation (67) support solitary wave solutions propagating on the background $\overline{u}$ and characterised by the speed-amplitude relation $c_{s}(\overline{u}+\mathcal{A},\overline{u})$, where $\mathcal{A}$ is the soliton amplitude measured from background. Now we introduce $k_{-}(\overline{u};\lambda)$ as the solution of the ODE

[TABLE]

with $\lambda$ a constant of integration. The number of solitary waves resulting from the temporal evolution of $u_{0}(x)$ can then be calculated as

[TABLE]

where $\lambda_{\infty}$ is obtained from the boundary condition $k_{-}(u_{\infty};\lambda_{\infty})=0$. For $u_{0}(x)$ in the form of a box of width $w$ and height $u_{m}$,

[TABLE]

For the solitary wave amplitudes, we have the generic formula in terms of the relationship between the integration constant $\lambda$ and the cutoff mean $\overline{u}$

[TABLE]

where $\lambda_{m}$ is defined by $k_{-}(u_{m};\lambda_{m})=0$. Here, we are assuming that $k_{-}(\overline{u};\lambda)$, hence $\mathcal{G}(\lambda)$, is an increasing function of $\lambda$. Then to obtain $\lambda=\lambda(\mathcal{A})$, one first solves the ODE

[TABLE]

where $\tilde{\omega}_{0}(\tilde{k},\overline{u})=-i\omega_{0}(i\tilde{k},\overline{u})$. The solution of equation (72) is $\tilde{k}_{+}(\overline{u};\tilde{\lambda})$, where $\tilde{\lambda}$ is a constant of integration. Setting $k(\overline{u};\lambda)=\tilde{k}(\overline{u};\tilde{\lambda})=0$ gives the relationship between $\lambda$ and $\tilde{\lambda}$. Substituting $\tilde{k}=\tilde{k}(\overline{u};\tilde{\lambda}(\lambda))$ into $\tilde{\omega}/\tilde{k}=c_{s}(\overline{u}+\mathcal{A},\overline{u})$ yields the desired $\lambda=\lambda(\mathcal{A})$ and $\mathcal{F}(\mathcal{A})=1-\mathcal{G}(\lambda(\mathcal{A}))$ for $\mathcal{A}\in[0,\mathcal{A}_{\rm max}]$ where $\mathcal{A}_{\rm max}$ satisfies

[TABLE]

Here, we are assuming that $\lambda$ is a decreasing function of $\mathcal{A}$.

Two of the main results of this paper do not depend on the system under study so long as the necessary pre-requisites of the solitary wave resolution method are satisfied. Equation (70) for a pure box initial condition is always linear in the box width. Also, the maximum solitary wave amplitude $\mathcal{A}_{\rm max}$ and the normalised cumulative amplitude distribution function for a box, $\mathcal{F}(\mathcal{A})$, are independent of box width. These results can be used, for example, to identify the initial box height that yields a desired lead solitary wave with amplitude $\mathcal{A}_{*}$ by solving $\lambda_{*}=\lambda(\mathcal{A}_{*})$ where $k_{-}(u_{*};\lambda_{*})=0$ for the box height $u_{*}$ and the box width $w_{*}=2\pi\mathcal{N}_{*}/k_{-}(u_{*};\lambda_{\infty})$ that results in the desired number of solitary waves $\mathcal{N}_{*}$.

5.5 Numerical methods

Direct numerical simulations of the conduit equation were undertaken following the method described in (Maiden & Hoefer, 2016). Equation (1) is rewritten as two coupled equations, the spatial discretisation utilises fourth-order finite differences with periodic boundary conditions, and the temporal evolution is via a medium-order Runge-Kutta method. Numerical results presented show how the long-time box evolution is altered by width in figure 6 and by height in figure 7. We observe that the number of solitary waves produced approximately changes linearly with the width but does not change significantly with box height past a certain height. We observe that the amplitude distributions change with box height but not with width.

We also numerically investigate our basic hypothesis that an initial, broad disturbance for the conduit equation results primarily in the fission of solitary waves. In order to quantify this, we consider one simulation that represents an edge case in which a smoothed box with amplitude $a_{m}=0.88$ and width $w=96$ results in a relatively small number (9) of solitary waves and agrees with the predicted number from equation (33). The initial and final (at $t=350$) profiles are shown in figure 8. Solitary waves travel faster than the long wave speed $c_{s}(a_{s},1)>2$, whereas dispersive waves propagate with the group velocity that is slower and exhibits a minimum $-1/4\leq\partial_{k}\omega_{0}(k,1)\leq 2$. We identify the dividing location in figure 8 between small amplitude dispersive waves and the solitary wavetrain as the first $z$ value, $z_{*}=1000$ here, in which the solitary wavetrain departs from unity. Over the entire domain $[0,L]$ ($L=1500$ here) and each of the two subintervals $[0,z_{*}]$, $[z_{*},L]$, we compute integrals of the conserved densities $a-1$, $1-1/a-a_{z}^{2}/a^{2}$ and the nonnegative density $(a-1)^{2}$ at the final time. The results are reported in table 2. In all cases, the small amplitude dispersive wave contributions are less than 1% of the total. Consequently, the solitary wavetrain dominates these integral quantities and our basic hypothesis is confirmed.

6 Comparison of modulation theory with experiment

Theory predictions for the number of solitary waves in physical experiments are reported in figure 9. We calculate the prediction $\mathcal{N}$ from a smoothed version of the viscous fluid conduit’s observed profile at the time of wavebreaking. Filled dots and the vertical axis in figure 9(a) report the number of observed solitary waves for each trial. We observe excellent agreement between experiment and theory, with the observed number of solitary waves being at most two away from the predicted value. Consequently, there is a decrease in the relative percent error as the total number of solitary waves increases, as shown in figure 9(b).

Figure 10 details comparisons between asymptotic predictions and physical observations of the solitary wave amplitude cumulative distribution functions. The prediction from KdV analysis is included. For the amplitude distribution $\mathcal{F}(\mathcal{A})$, any relative minimum in the initial profile results in unphysical predictions. Therefore, instead of using the true profile generated from boundary control, we use an averaged version, similar to that used in numerical experiments (see the inset of figure 1). We fit the box amplitude $a_{m}$ and width $w$ by extracting these values from the experimental time, location, and mean height of breaking as determined by our previously introduced inflection point criterion (Anderson et al., 2019). The box profile is approximated by the following formula

[TABLE]

where the tanh nondimensional width 2.5 was identified as a good fit across all reported trials to both the leading edge transition and the final amplitude distributions. Although our analysis is explicit for sharp box profiles, we find that smoothing the box does slightly influence the smaller amplitude soliton distribution as depicted in figure 11 (see dashed versus dash-dotted curves).

We also observe a change in conduit diameter of roughly $15\%$ from the bottom of the apparatus to the location of solitary wave data-taking. While this does not affect $\mathcal{N}$, this is observed to have a profound affect on $\mathcal{F}(\mathcal{A})$. To compensate for this, we use the amplitude prediction from (Maiden et al., 2018) for a solitary wave on a changing background for the conduit equation. We measure error via the $\infty$-norm, as this relates to the Kolmogorov-Smirnov test for comparing cdfs in statistics. Across all trials, the conduit prediction has roughly half the error as the prediction from the rescaled KdV prediction.

We also compare our results to the explicit formula for a box (64). We observe in figure 11 that, as the initial condition’s width increases—corresponding to a larger number of solitary waves therefore improving the asymptotic approximation—the observed distribution approaches the expected distribution that is independent of width.

Our final comparison between experiment and theory involves the spatio-temporal data set reported in figure 2(b). Utilising the nominal, measured experimental parameter values reported in the caption of that figure, we determine the length and time scalings for the conduit equation in (13) to be $R_{0}/\sqrt{8\epsilon}=$1.6\text{,}\mathrm{mm} and $R_{0}/(U_{0}\sqrt{8\epsilon})=$1.17\text{\,}\mathrm{s}, respectively. These scalings and the measured parameters determine the initial, nondimensional box width $w=156$ and height $a_{m}=1.6$. A numerical simulation of the conduit equation with this smoothed box initial data (c.f. equation (74)) and these scalings is shown in figure 12(a). We report the nondimensional diameter $\sqrt{a}$ in the figure (black curves) in order to directly compare the numerics with experiment. The conduit equation simulation domain was taken to be larger ($160\text{\,}\mathrm{cm}$ in 12(a) and $180\text{\,}\mathrm{cm}$ in 12(b)) than the view shown so that a portion of the initial box is outside the displayed viewing window.

The experimental diameter profiles $D(Z,T)$ reported in figure 12 have been extracted from the images in figure 2(b) per the description in Section 2.2. Because of the large aspect ratio inherent in these long wave dynamics, low image resolution in the transverse direction implies that 8\text{\,}\mathrm{pixel}$\leq D\leq$20\text{\,}\mathrm{pixel}. However, from the box and wave cameras, we have much more accurate measurements of the conduit diameter near the bottom and top of the apparatus for both the equilibrium conduit diameter ($2\text{\,}\mathrm{mm}$) and the diameter of the box ($3.2\text{\,}\mathrm{mm}$). We use these measurements to normalise the pixel data via the linear transformation $D^{\prime}(Z,T)=\alpha D(Z,T)+\beta$. The coefficients $\alpha=$0.11\text{,}{\mathrm{pixel}}^{-1}$$ and $\beta=0.30$ are chosen so that the mean equilibrium conduit is $\overline{D^{\prime}}=1$ and the box diameter satisfies $D^{\prime}=3.2/2=1.6$. This normalisation effectively registers the low resolution data in figure 2 with the high resolution measurements from the other cameras. The profiles $D^{\prime}(Z,T)$ are reported in figure 12 with the lighter gray curves.

While the experiment gives rise to 14 solitary waves, the numerics produce 16. But the lead solitary wave diameter is only 4.5% larger than experiment. The numerical evolutionary timescale is somewhat slower than the experimental one, which is consistent with previous measurements of large amplitude solitary waves that were found to propagate faster than the conduit equation’s solitary wave speed-amplitude relation (10) predicts (Olson & Christensen, 1986; Maiden et al., 2016).

In figure 12(b), we utilise the same measured parameters as in (a) except we fit the interior $\mu^{(i)}=$6.88\text{\times}{10}^{-2}\text{,}\mathrm{Pa}\text{,}\mathrm{s} and exterior $\mu^{(e)}=$0.9\text{\,}\mathrm{Pa}\text{\,}\mathrm{s} viscosities by increasing the nondimensional length scale by the factor $16/14$ and reducing the nondimensional time scale by the factor $0.97$. This particular increase in length scale derives from the predicted linear scaling of the number of solitary waves by the initial box width. Indeed, figure 12(b) exhibits precisely 14 solitary waves from both numerics (nondimensional initial box width $w=137$) and experiment. The increased length scale and slightly reduced time scale lead to significantly improved solitary wavetrain evolution when compared with experiment. Remarkably, at the final reported time $t=$177\text{,}\mathrm{s}$$, the numerical and experimental normalised diameter profiles are almost indistinguishable for the 11 largest solitary waves in the solitary wavetrain.

For the fit, the exterior viscosity is reduced by 10%. The high viscosity glycerine utilised for the exterior fluid is extremely sensitive to even small amounts of interior fluid mass diffusion from the conduit. A 10% reduction from its nominal value is certainly possible. The interior viscosity’s fitted value is approximately 39% larger than its measured value, which is a bit more than expected. However, we have not accounted for uncertainty in the volumetric flow rate or fluid densities. Moreover, the conduit equation is a long-wave approximation of the full two-fluid, free-boundary dynamics that is valid in the small viscosity ratio $\epsilon=\mu^{(i)}/\mu^{(e)}$ regime (Lowman & Hoefer, 2013a). For these experiments, the measured value of this ratio is $\epsilon=0.049$ and the fitted value is $\epsilon=0.076$. Despite these reasonably small nondimensional values, a number of higher order effects, e.g., inertia and the finite-sized boundary (Lowman & Hoefer, 2013a), could be influencing the dynamics on the long timescales considered—the nondimensional final time is approximately 150 in both figures 12(a) and 12(b). For these reasons, we find the comparison between experiment and the conduit equation reported in figure 12(b) to be credible, strong evidence for both the conduit equation as an accurate model of viscous fluid conduit dynamics and the efficacy of the solitary wave resolution method.

7 Conclusion

We have derived explicit formulae accurately predicting the asymptotic number and amplitude distribution of solitary waves that emerge in long time from a localised, slowly varying initial disturbance for the conduit equation. Our analytical approach to the solitary wave fission problem is based upon Whitham modulation theory. The predictions have been quantitatively verified with experiments on the interfacial dynamics of two high-contrast, viscous fluids. While the solitary wave resolution method utilised here was previously developed by El et al. (2008) for the Serre-Green-Naghdi equations modeling large amplitude, shallow water waves, we have identified two new, universal predictions for box-shaped initial disturbances: i) the asymptotic number of solitary waves is linearly proportional to box width, ii) the asymptotic, normalised cumulative distribution function for the solitary wave amplitudes is independent of box width. Universal here means that these predictions apply to a broad class of dispersive hydrodynamic equations (67), not just the conduit equation.

Our experiments are the first that validate the solitary wave resolution method and we find that the physical evolution of viscous conduit profiles is well captured by the approximation, particularly for large width disturbances. All observed solitary wave counts are within 1–2 waves of their expected values and within 10% relative error for disturbances producing at least twelve waves. Amplitude distribution predictions agree quantitatively with experiment and demonstrate the necessity of going beyond the standard weakly nonlinear KdV model as it applies to the viscous fluid conduit system. The conduit’s observed, full spatio-temporal evolution exhibits remarkable fidelity to numerical simulations of the conduit equation when the two fluids’ viscosities are appropriately fitted to effectively account for a variety of uncertainties and higher order effects in the experiment. When this work is considered in conjunction with the large variety of previous experiments on viscous fluid conduits involving solitary waves (Olson & Christensen, 1986; Scott et al., 1986; Whitehead & Helfrich, 1988; Helfrich & Whitehead, 1990; Lowman et al., 2014), wavebreaking (Anderson et al., 2019), rarefaction waves, and dispersive shock waves (Maiden et al., 2016, 2018), the analytically tractable conduit equation, and the effectiveness of Whitham modulation theory further bolster the claim that the viscous fluid conduit system provides both an ideal laboratory environment and mathematical modelling framework in which to examine dispersive hydrodynamics and nonlinear dispersive wave dynamics more generally (Lowman & Hoefer, 2013a; Maiden & Hoefer, 2016).

Largely due to the paucity of analytical tools for studying strongly nonlinear wave dynamics, researchers have focused primarily on integrable models such as the KdV and modified KdV (mKdV) equations to obtain physical predictions for the soliton fission problem. A case in point is the field of internal ocean waves where strongly nonlinear solitary waves are known to be prevalent yet can only be explained with fully nonlinear models (Helfrich & Melville, 2006). Because the solitary wave resolution method is agnostic to integrable structure, a promising application direction is the fully nonlinear Miyata-Camassa-Choi (MMC) equations for two stratified fluid layers (Miyata, 1985; Choi & Camassa, 1999). The MMC equations have all the necessary ingredients to apply the solitary wave resolution method (Esler & Pearce, 2011).

These results quantify an extension of the soliton resolution conjecture from integrable systems to non-integrable systems that, often times, more closely model physical systems. The conjecture states that localised initial conditions to nonlinear, dispersive wave equations generically evolve in long time toward a rank-ordered train of solitary waves separated from small amplitude, dispersive waves (Segur, 1973; Schuur, 1986; Deift et al., 1994). Here, we have formally derived quantitative measures of the solitary wave component, which typically dominates the long-time outcome in problems with large-scale, localised initial data. Predicting properties of the accompanying small amplitude dispersive waves is the next step toward quantifying an extension of the full conjecture to non-integrable systems.

These promising results for solitary wave fission in the model conduit system provide inspiration for future studies on this problem in other complex systems.

This work was supported by NSF DMS-1255422 and DMS-1517291 (M. A. H.), the NSF GRFP (M. D. M.), NSF EXTREEMS-QED DMS-1407340 (E. G. W.), and EPSRC grant EP/R00515X/2 (G. A. E.).