## Abstract

Modulated periodic interfacial waves along a conduit of viscous liquid are explored using nonlinear wave modulation theory and numerical methods. Large-amplitude periodic-wave modulation (Whitham) theory does not require integrability of the underlying model equation, yet often either integrable equations are studied or the full extent of Whitham theory is not developed. Periodic wave solutions of the nonlinear, dispersive, non-integrable conduit equation are characterized by their wavenumber and amplitude. In the weakly nonlinear regime, both the defocusing and focusing variants of the nonlinear Schrödinger (NLS) equation are derived, depending on the carrier wavenumber. Dark and bright envelope solitons are found to persist in long-time numerical solutions of the conduit equation, providing numerical evidence for the existence of strongly nonlinear, large-amplitude envelope solitons. Due to non-convex dispersion, modulational instability for periodic waves above a critical wavenumber is predicted and observed. In the large-amplitude regime, structural properties of the Whitham modulation equations are computed, including strict hyperbolicity, genuine nonlinearity and linear degeneracy. Bifurcating from the NLS critical wavenumber at zero amplitude is an amplitude-dependent elliptic region for the Whitham equations within which a maximally unstable periodic wave is identified. The viscous fluid conduit system is a mathematically tractable, experimentally viable model system for wide-ranging nonlinear, dispersive wave dynamics.

## 1. Introduction

Nonlinear wave modulation is a major mathematical component of the description of dispersive hydrodynamic phenomena. Dispersive hydrodynamics encompass the study of fluid-like media where dissipative effects are weak compared to dispersion [1]. Solitary waves and dispersive shock waves (DSWs) are typical coherent structures. Model equations include the integrable Korteweg–de Vries (KdV) and nonlinear Schrödinger (NLS) equations as well as non-integrable counterparts that are important for applications to superfluids, geophysical fluids and laser light. Modulation theory assumes the existence of a multi-parameter family of nonlinear, periodic travelling wave solutions whose parameters change slowly relative to the wavelength and period of the periodic solution under perturbation. Such variation is described by modulation equations. In the weakly nonlinear regime, the NLS equation is a universal model for the slowly varying envelope and phase, incorporating both cubic nonlinearity and dispersion. In the large-amplitude regime, the Whitham equations [2] describe slow modulations of the wave's mean, amplitude and wavenumber. At leading order, they are a dispersionless system of quasi-linear equations.

In this paper, we investigate nonlinear wave modulations in both the weakly nonlinear and large-amplitude regimes for the dimensionless conduit equation
*A* at time *t* and vertical spatial coordinate *z*, separating a light, viscous fluid rising buoyantly through a heavier, more viscous, miscible fluid at small Reynolds numbers [3,4]. Our motivation for studying equation (1.1) is twofold. First, the conduit equation is not integrable [5], so there are mathematical challenges in analysing its rich variety of nonlinear wave features. Second, equation (1.1) is an accurate model of viscous fluid conduit interfacial waves where hallmark experiments have been performed on solitary waves [6–8], their interactions [3,9,10] and DSWs [11]. We believe the conduit system is an ideal model for the study of a broad range of dispersive hydrodynamic phenomena.

Indeed, in this work, we elucidate additional nonlinear wave phenomena predicted by equation (1.1) by analysing its weakly nonlinear, NLS reduction, the structural properties of the large-amplitude Whitham equations, and numerical simulations. An example of a non-modulated periodic wave is shown in figure 1 that defines the wave's mean *a*, wavenumber *k*, minimum *ϕ*_{0}, and phase velocity *c*_{p}=*ω*/*k* (*ω* is the wave's frequency). We identify persistent bright and dark envelope soliton solutions in both the weakly nonlinear, NLS regime and the large-amplitude, strongly nonlinear regime. The quasi-linear Whitham equations are analysed asymptotically and with numerical computation. Regions in parameter space of strict hyperbolicity, ellipticity and linear degeneracy are identified. The elliptic regime corresponds to modulationally unstable periodic waves and a maximally unstable wave is identified.

### (a) Background on viscous fluid conduits

Conduits generated by the low Reynolds number, buoyant dynamics of two miscible fluids with differing densities and viscosities were first studied in the context of geological and geophysical processes [12]. A system of equations describing the dynamics of melted rock within a solid rock matrix was derived by treating molten rock and its solid, porous surroundings as two fluids with a large density and viscosity difference [13]. Under appropriate assumptions, the family of magma equations
*φ* of molten rock, can be derived [14,15]. There are two constitutive model parameters (*n*,*m*) that relate the porosity of the rock matrix to its permeability and viscosity, respectively. Physical values of the parameters are 2≤*n*≤5, 0≤*m*≤1. The conduit equation (1.1) happens to coincide with the magma equation (1.2) when (*n*,*m*)=(2,1) [3]. For this reason, viscous fluid conduits were utilized as a laboratory model of magma dynamics. The conduit equation can be derived directly from the Stokes equations for two viscous fluids under a long-wave, high-viscosity contrast assumption [4]. Viscous fluid conduits, contrary to magma, are easily accessible in a laboratory setting, typically with a sugar solution or glycerine for the exterior fluid, and a dyed, diluted version of the same for the interior fluid [3,6,9,11].

Early experiments primarily explored the development of the conduit itself, which results in a diapir followed by a periodic wavetrain [12,6]. Solitary waves in an established conduit have also been extensively studied, including their amplitude–speed relation, interactions and fluid transport properties [3,6–8,10]. Experiments have also shown that interactions between solitons are nearly elastic, with a phase shift the primary quantifiable change [8,10]. Furthermore, soliton interaction geometries predicted by Lax for the KdV equation [16] were observed and agreed well with numerical simulations [10]. This is particularly notable because the Lax categories persisted into the large-amplitude regime, although unlike KdV, the short-wavelength behaviour of conduit dispersion is bounded, akin to the Benjamin–Bona–Mahoney (BBM) equation [10]. Recently, dispersive shock waves were observed, yielding good agreement with predictions from Whitham averaging theory [11]. The accompanying observations of soliton–DSW interaction suggest a high degree of coherence, i.e. the sustenance of dissipationless/dispersive hydrodynamics over long spatial and temporal time scales. It is for this reason that we further investigate modulations of periodic conduit waves.

### (b) Properties of the conduit equation

To fully describe the two-fluid interface of the conduit system, one can consider the full Navier–Stokes equations with boundary conditions along a moving, free interface. However, in the low Reynolds number, small interfacial slope, long-wave regime, a balance between the viscous stress force of the exterior fluid and the buoyancy force acting on the interior fluid leads to the asymptotically resolved conduit equation (1.1) with no amplitude assumption [6,4]. The force balance is achieved with small interfacial slopes on the order of the square root of the ratio of the interior to exterior fluid viscosities.

The conduit equation (1.1) has been studied since the 1980s and is known to have exactly two conservation laws [17,18]:
*k*, similar to the bounded dispersion of the BBM equation. This leads to the linear phase *c*_{p} and group *c*_{g} velocities
*c*_{g}<*c*_{p} for *k*>0. While the phase velocity is always positive, the group velocity is negative for *k*>1. Failure of the Painlevé test suggests that the conduit equation is not completely integrable [5]. The conduit equation is globally well-posed for initial data *A*(*z*,0) physically relevant data bounded away from zero in order to avoid the singularity [19].

Solitary waves have been studied numerically for the more general magma equation (1.2), where it has been found that they exhibit near-elastic interactions resulting in a phase shift and a physically negligible dispersive tail [10,20,21]. The asymptotic stability of solitary waves has also been proved [22]. General (unmodulated) periodic wave solutions have been found, and an implicit dispersion relation has been computed for these waves [6]. In the long-wavelength, small-amplitude regime, the conduit equation reduces to KdV [23]. However, the fact that solitary waves exhibit KdV-like interaction behaviour including almost elastic interactions, the three Lax interaction categories [10] and coherent interactions with DSWs [11], all for strongly nonlinear, large-amplitude solitary waves is notable.

There have been several works that apply Whitham modulation theory to the magma equations (1.2). Marchant & Smyth [24] considered equation (1.2) with (*n*,*m*)=(3,0), describing DSWs and some structural properties of the Whitham equations. Modulations of periodic travelling waves in the magma equation (1.2) and a generalization of it were investigated in the weakly nonlinear, KdV regime [25]. Modulated periodic waves in the form of DSWs were investigated for the entire family of magma equations (1.2) in [26]. This work differs from previous studies by concentrating on the case (*n*,*m*)=(2,1) for the conduit equation, identifying new coherent structures (envelope solitons) and determining structural properties of the associated Whitham equations (hyperbolicity, ellipticity, linear degeneracy).

### (c) Outline of this work

The paper is organized as follows. Periodic travelling wave solutions to the conduit equation are studied in §2 both numerically and asymptotically in the weakly nonlinear regime. In §3, we consider weakly nonlinear periodic wave modulations and include long-wave dispersion to derive the NLS equation. By an appropriate choice of the periodic travelling wave's wavenumber, both the focusing and defocusing variants of the NLS equation are possible. We numerically demonstrate the persistence of large-amplitude dark and bright envelope solitary wave solutions in the full conduit equation. In §4, we analyse modulated periodic waves of arbitrary amplitude via the conduit Whitham modulation equations. A weakly nonlinear analysis and direct numerical computation are used to determine structural properties of the Whitham equations including hyperbolicity or ellipticity and genuine nonlinearity or linear degeneracy. Consequences for the stability of periodic waves are examined. The manuscript is concluded with a discussion of the implications of this work in §5.

## 2. Periodic travelling wave solutions

We seek periodic travelling wave solutions to equation (1.1) in the form *A*(*z*,*t*)=*ϕ*(*θ*), *θ*=*kz*−*ωt*, *ϕ*(*θ*+2*π*)=*ϕ*(*θ*) for *A* and *B* are real integration constants.

The above equation exhibits at most three real roots [26]. When there are three distinct roots, a periodic solution oscillates between the largest two. The solution can therefore be parametrized by three independent variables. Defining the wave minimum *ϕ*_{0} according to *ϕ* is 2*π*-periodic is enforced through
*g*(*ϕ*_{0})=*g*(*ϕ*_{0}+*a*)=0 determine *A* and *B* from equation (2.2).

Owing to the scaling invariance equation (1.4), the wave mean can be scaled to unity. This implies that only *ω*(*k*,*a*,1) and *ϕ*_{0}(*k*,*a*,1) need be determined. Then the general cases follow according to

### (a) Stokes expansion

We can obtain approximate periodic travelling wave solutions in the weakly nonlinear regime via the Stokes wave expansion [27]:
*ε*≪1 is an amplitude scale. Inserting this ansatz into equation (2.1), equating like coefficients in *ε* and enforcing solvability conditions yield the approximate solution
*a*,*b* shows the dispersion and phase velocities for numerically computed periodic waves, and figure 3*c* compares the full, nonlinear dispersion

Of interest is that the approximate solution (2.10) can result in an unphysical, negative conduit cross-sectional area. The minimum of the approximate solution *θ*=*π*. Equating the minimum to zero, we find that physical, positive values for approximate

## 3. Weakly nonlinear, dispersive modulations

The aim of this section is to describe wave modulation in the weakly nonlinear, dispersive regime. The approximate NLS equation is found using multiple scales in appendix B by seeking a solution in the form [2,28]
*ε* is an amplitude scale. By introducing the standard, scaled coordinate system
*B*(*ζ*,*τ*)
*ζ*_{0}, *ψ*_{0} and *α*, where 0<*α*≤*π*/2 is the phase jump across the soliton. The grey soliton reduces to a black soliton when *α*=*π*/2.

For the focusing case, periodic waves are modulationally unstable [2,29]. Bright envelope solitons for the NLS equation exist, which have the form
*ζ*_{0} and *Θ*_{0} are arbitrary, real constants. To validate these approximate solutions, we numerically simulate the conduit equation (1.1) with envelope soliton initial conditions (3.1) and (3.4) or (3.1) and (3.5), depicted in figure 4*a*,*b*. In figure 4*a*, a black soliton is observed to coherently propagate, maintaining essentially the same shape. The NLS approximation is asymptotically valid up to times *a*, *ε*^{2}=400. The black envelope soliton shows no sign of instability over times up to *t*=1000. Figure 4*b* shows the long-time evolution of an envelope bright soliton. The envelope appears to steepen and become peaked by *t*=1000 but otherwise maintain its essential structure. The observed speeds of propagation of the black and bright solitons are 0.0002 and −0.1589 very close to the predicted group velocities 0 and −0.16, respectively.

We also numerically studied the large-amplitude regime with dark and bright envelope soliton initial conditions. Figure 5*a* shows the numerical evolution of dark envelope soliton initial data (3.1), (3.4) for *t*=750 over the truncated domain *z*∈[100,175] and use this as an initial condition for the conduit equation. Additional periods of the unmodulated wave were prepended to the profile in order to increase the spatial domain. The evolution is shown in figure 5*b*, displaying a remarkable coherence and persistence out to *t*=1000.

Figure 6*a* depicts the evolution of bright envelope soliton initial conditions with *t*=1000 over the truncated domain *z*∈[350,360] and superimposed on a unit background and then used as a new initial condition for the conduit equation. The result is shown in figure 6*b* that shows the persistence of a large-amplitude envelope structure accompanied by the emission of small-amplitude dispersive radiation.

These results reflect the fact that the NLS approximation models the envelope of a weakly nonlinear, dispersive carrier wave. Yet numerical evolution of large-amplitude initial data present intriguing coherent structures deserving of further study. We now turn to the Whitham equations for an asymptotic description of nonlinear wave modulations in the moderate- to large-amplitude regime. However, the trade-off for using these quasi-linear equations is their lack of dispersion at the first order of approximation; consequently, they cannot describe envelope solitons.

## 4. Whitham equations

To describe modulated, large-amplitude periodic waves, we appeal to the Whitham modulation equations. Whitham's original formulation invoked averaged conservation laws [30], later shown to be equivalent to a perturbative, multiple-scales reduction [31]. For completeness, we have implemented both approaches in appendix C. For this, we seek modulations of an arbitrary amplitude, 2*π*-periodic, travelling wave solution *ϕ* to equation (1.1) (see §2). As we will incorporate only the leading-order Whitham equations, the large parameter 1/*ε* here corresponds to the time scale of their validity. Note that the lack of an amplitude restriction in leading order Whitham modulation theory is at the expense of a shorter time scale of validity relative to the

It is convenient to express the Whitham equations in the conservative form
*ϕ*′^{2}=*g*(*ϕ*) (2.2) and use the notation

The scaling invariance (1.4) can be used so that the dependence on *i*=1,2,3. Therefore, computation of the averaging integrals is only required for

We will be interested in structural properties of the Whitham equations such as hyperbolicity (strict or non-strict), ellipticity and genuine nonlinearity. All of these criteria rely on the eigenvalues *c* and eigenvectors *c*_{1}≤*c*_{2}≤*c*_{3} when the Whitham equations are hyperbolic or, in the case of one real and two complex conjugate eigenvalues, the Whitham equations are elliptic. If the eigenvalues are all real and they are strictly ordered *c*_{1}<*c*_{2}<*c*_{3}, then the Whitham equations are strictly hyperbolic.

The coefficient matrix

Using the identities in (4.9), we find that the quantity
*μ*≠0, then the Whitham equations are genuinely nonlinear [2]. For those values of *μ*=0, the Whitham equations are linearly degenerate. The sign definiteness of *μ* corresponds to a monotonicity condition that is required for the existence of simple wave solutions to the Whitham equations, of particular importance for the study of DSWs [1].

### (a) Weakly nonlinear regime

Now consider equation (4.4) in the small *a* regime by inserting the Stokes expansion (2.10), (2.11), yielding

Next, we determine the approximate eigenvectors *μ*_{j}, *j*=1,2,3 (4.11). The expressions are cumbersome, so we do not report them here but there are two noteworthy findings. We find that *μ*_{1}=0 when

We also find linear degeneracy in the second characteristic field: *μ*_{2}=0 when

### (b) Large-amplitude regime

We now investigate modulations of large-amplitude, periodic waves by direct computation of the Whitham equations. For this, we examine the Whitham equations in the form (4.6), so that the dependence on *j*,*l*=1,2,…,*N*. We chose *N*=4000, Δ=0.001 so that

Using our direct computation of the coefficient matrix *a* the region in the

As noted earlier, ellipticity of the Whitham equations implies modulational instability (MI) of the periodic travelling wave [2]. In agreement with our weakly nonlinear analysis (4.13), we find that, in the elliptic region, *t*=500, whichever was longer. Some waves, especially those in the small-amplitude regime, were evolved for even longer time periods. The modulational (in)stability of several of these runs are shown in figure 7*a*. We find excellent agreement with the MI predictions from Whitham theory. The long-time evolution of two particular waves are shown in figure 8, showing both a stable and an unstable case. The unstable case in figure 8*b* appears to show the formation of large-amplitude, bright envelope coherent structures. This is additional numerical evidence for the existence of bright envelope soliton solutions of the conduit equation.

A periodic travelling wave solution of the conduit equation (1.1) corresponds to a constant solution *κ*. Because the Whitham equations are quasi-linear, first-order equations, any wavenumber is permissible (determined by the initial data), suggesting that the physical growth rate (4.15) is unbounded. In practice, there is a dominant wavenumber associated with the instability that is determined by higher order effects, which in this case would be higher order dispersion in the Whitham equations. The NLS equation (3.3) resolves this feature in the weakly nonlinear regime, but we are interested in large-amplitude modulations. We therefore identify the imaginary part of the characteristic velocity *b* that there is a maximum of *κ* in (4.15).

Next, we compute the quantities *a* where the curves correspond to the largest value of *μ* changes sign. The curve where *μ*_{1} changes sign bifurcates from the edge of the elliptic region at the point *μ*_{2,3} change sign apparently bifurcates from (0,0) and occurs for small

To understand the small *b* a zoom-in of this region. The accurate determination of the loss of genuine nonlinearity in this region is numerically challenging because the characteristic velocities *b* where, for each *μ*_{1}=0 in figure 9*a* occurs in a strictly hyperbolic region.

## 5. Discussion/conclusion

Our study of the structural properties of the conduit Whitham equations sheds some light on recent theoretical and experimental studies of dispersive shock waves. The DSW fitting method allows one to determine a dispersive shock's harmonic and soliton edge speeds, even for non-integrable systems [36]. However, the method is known to break down when the Whitham equations lose genuine nonlinearity in the second characteristic field [26,1]. It was observed in [26] that the fitting method failed to accurately predict conduit DSW soliton edge speeds for sufficiently large jump heights. Our results here suggest that this could be due to the loss of strict hyperbolicity and/or genuine nonlinearity in the small wavenumber (soliton train) regime.

In addition to the hyperbolic modulation regime where DSWs and dark envelope solitons are prominent coherent structures, we have found an elliptic regime where the periodic wave breaks up into coherent wave-packets or bright envelope solitons. The accessibility of both hyperbolic and elliptic modulation regimes in one system motivates further study of each and the transition between the two. One potential future, novel direction is to explore the possibility of creating a soliton gas [37].

It remains to generate a periodic wave from an initially uniform conduit and explore its properties experimentally. Accurate control of wave-breaking via a dispersionless simple wave (a rarefaction wave) has been achieved by slow modulation of the conduit area from a boundary [11]. One generalization of this is to use simple wave solutions of the Whitham equations to efficiently and smoothly transition between a constant conduit

We have presented strong numerical evidence for the existence of large-amplitude dark and bright envelope solitary waves in viscous fluid conduits, bifurcating from weakly nonlinear NLS solutions. Dark envelope solitons can have either positive or negative velocities. All bright envelope solitons for NLS have negative velocities. It remains to be determined if this is the case in the large-amplitude case. Existing laboratory studies of viscous fluid conduits implement control of the conduit interface by varying the injected flow rate through a nozzle at the bottom of a fluid column. This allows for the creation of waves with positive (upward) propagation velocities such as dark envelope solitons. If bright envelope solitons only have negative velocities, then an alternative experimental approach will be required to create them.

We have shown that the non-convexity of the conduit linear dispersion relation leads to the existence of elliptic Whitham equations and modulational instability. This is just one possible implication of non-convex dispersion in dispersive hydrodynamics. We note that non-convex dispersion in other, higher order equations has also been found to give rise to a resonance between the DSW soliton edge and linear waves, leading to the generation of radiating DSWs [38–40].

This study has identified and categorized modulations of both small- and large-amplitude periodic travelling waves for the conduit equation (1.1). These findings, along with previous theoretical and experimental studies of solitons and DSWs in the strongly nonlinear regime imply that the viscous fluid conduit system is an accessible environment in which to investigate rich and diverse nonlinear wave phenomena.

## Data accessibility

This paper contains no experimental data. All computational results are reproducible.

## Authors' contributions

Both authors contributed equally to problem formulation, solution and writing of the paper.

## Competing interests

We have no competing interests.

## Funding

This work was partially supported by NSF CAREER DMS-1255422 (M.A.H.) and NSF GRFP (M.D.M.).

## Appendix A. Numerical methods

**(a) Periodic solutions**

We compute unit-mean conduit periodic travelling wave solutions *A*, which we do not require because of our imposition of unit mean. Projection of equation (A1) onto *n*=1,…,*N* yields *N* equations for the *N*+1 unknowns *k*+1, i.e. by ^{−13} in the 2-norm of the residual and choosing *N* so that |*a*_{n}| is below 5⋅10^{−12} for *n*>3*N*/4. The number of coefficients required strongly depends on the wavenumber *N*=2^{6} provides sufficient accuracy, whereas for *N*=2^{12}.

With the cosine series coefficients of *j*=1,2,3 in equations (4.5), (4.2) using the spectrally accurate trapezoidal rule. We then use sixth-order finite differencing to compute derivatives of

**(b) Time-stepping**

For the direct numerical simulation of the conduit equation (1.1), it is convenient to write it in the form of two coupled equations:
*A*(0,*t*) and *A*(*L*,*t*) so that the first equation in (A2) yields the boundary conditions for

## Appendix B. Nonlinear Schrödinger equation derivation

Here, we derive an approximation of wave modulations in the small-amplitude, weakly nonlinear regime. Consider the ansatz
*A*_{i}=*A*_{i}(*θ*,*Z*,*T*) for *i*=0,1,…, *Z*=*εz* and *T*=*εt*, where *A*_{0}: *A*_{0}=*ψ*(*Z*,*T*) *e*^{iθ}+*c*.*c*., where ‘c.c.’ denotes the complex conjugate of the previous terms. At *g*_{1}. Solving for *A*_{1}, we include second harmonic and mean terms *M* to be determined at the next order. Solvability with respect to constants at *g*_{1}, which, upon entering the moving reference frame and scaling to long time

## Appendix C. Derivation of the Whitham equations

For completeness, we supply a synopsis of the multiple scales asymptotic derivation of the Whitham modulation equations. For the formal derivation here, we introduce slow space and time scales *Z*=*εz*, *T*=*εt* and consider the ansatz *A*(*z*,*t*)=*A*_{0}(*θ*,*Z*,*T*)+*εA*_{1}(*θ*,*Z*,*T*)+⋯ , 0<*ε*≪1, where *θ*_{z}=*k* and *θ*_{t}=−*ω*. Continuity of mixed partials *θ*_{zt}=*θ*_{tz} implies the conservation of waves *k*_{T}+*ω*_{Z}=0, one of the Whitham equations. We insert the asymptotic ansatz into the conduit equation (1.1) and equate like orders in *ε*. The *Z*,*T*). Note that in order to remove secularity at this order, the period of the solution must be scaled to a constant [31,2], which we choose to be 2*π* without loss of generality.

At the next order,

There are two solvability conditions in the form *π*-periodic boundary conditions. Note that there is a third, linearly independent function that is annihilated by *π*-periodic. Applying the two solvability conditions *ε*=1, i.e. considering the Whitham equations as the long time *t*≫1 asymptotic, we obtain equations (4.1).

Averaging of the conservation laws (1.3) is achieved by inserting the ansatz *A*(*z*,*t*)=*ϕ*(*θ*) and averaging the densities and fluxes over a period:

- Received July 1, 2016.
- Accepted November 14, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.