Stability of Viscous St. Venant Roll Waves: From Onset to Infinite Froude Number Limit

Abstract

We study the spectral stability of roll wave solutions of the viscous St. Venant equations modeling inclined shallow water flow, both at onset in the small Froude number or “weakly unstable” limit $F\rightarrow 2^+$ and for general values of the Froude number F, including the limit $F\rightarrow +\infty$. In the former, $F\rightarrow 2^+$, limit, the shallow water equations are formally approximated by a Korteweg-de Vries/Kuramoto–Sivashinsky (KdV–KS) equation that is a singular perturbation of the standard Korteweg-de Vries (KdV) equation modeling horizontal shallow water flow. Our main analytical result is to rigorously validate this formal limit, showing that stability as $F\rightarrow 2^+$ is equivalent to stability of the corresponding KdV–KS waves in the KdV limit. Together with recent results obtained for KdV–KS by Johnson–Noble–Rodrigues–Zumbrun and Barker, this gives not only the first rigorous verification of stability for any single viscous St. Venant roll wave, but a complete classification of stability in the weakly unstable limit. In the remainder of the paper, we investigate numerically and analytically the evolution of the stability diagram as Froude number increases to infinity. Notably, we find transition at around $F=2.3$ from weakly unstable to different, large-F behavior, with stability determined by simple power-law relations. The latter stability criteria are potentially useful in hydraulic engineering applications, for which typically $2.5\le F\le 6.0$.

Appendices

Appendix 1: Notation for High-Frequency Bounds

The proof of Lemma 2.2 in Sect. 2.2.1 relied heavily on computations previously carried out in detail in Section 4.1 of Barker et al. (2011). There the authors were concerned with the stability analysis of traveling solitary wave solutions of the viscous St. Venant equation (1.11), i.e., those traveling wave solutions that decay exponentially fast to zero as $x-\bar{c}t\rightarrow \pm \infty$. By a straightforward adaptation of this analysis to the periodic case,²⁴ it follows that, for each $k\in \mathcal {P}$, there exists an $X_\delta$-periodic change in variables $W(\cdot )=P(\cdot ;\lambda ,\delta )Z(\cdot )$ that transforms the spectral problem (2.13) into an equivalent system of the form $W'(x)=\left( D(x,\lambda )+\Theta (x,\lambda )\right) W(x)$, supplemented with the boundary conditions $W(X_\delta )=e^{i\xi }W(0)$ for some $\xi \in [-\pi /X_\delta ,\pi /X_\delta )$, where the $3\times 3$ matrix-valued D is defined via $D(\cdot ,\lambda )=\mathrm{diag }\left( \frac{\lambda }{\bar{c}}+\theta _0+\frac{\theta _1}{\lambda },~~\bar{\tau }\sqrt{\frac{\lambda }{\nu }},~~-\bar{\tau }\sqrt{\frac{\lambda }{\nu }}\right) ,$ where $\theta _0=\frac{\bar{\alpha } \bar{\tau }^2}{\bar{c}\nu }$ and $\theta _1=-\frac{\bar{\tau }^2(q-\bar{c}\bar{\tau })^2}{\nu }-\frac{\bar{c}\bar{\alpha }}{\sqrt{\nu }}\frac{\bar{\tau }^3}{\nu ^{3/2}}$ with $\bar{ \alpha }:=\bar{\tau }^{-3}(F^{-2}+2\bar{c}\nu \bar{\tau }')$. Moreover, the $3\times 3$ matrix $\Theta (x,\lambda )$ has the block structure $\Theta =\left( \begin{array}{cc} \Theta _{++} &{} \Theta _{+-}\\ \Theta _{-+} &{} \Theta _{--} \end{array}\right) ,$ where $\Theta _{++}$ is a $2\times 2$ matrix. The individual blocks of the matrix $\Theta$ can be expanded as cubic polynomials in $\lambda ^{-1/2}$ with matrix-valued coefficients. More precisely, they expand as

$$\begin{aligned} \left\{ \begin{aligned} \Theta _{++}(\cdot ,\lambda )&=\Theta _{++}^0 + \lambda ^{-1/2}\Theta _{++}^1 + \lambda ^{-1}\Theta _{++}^2 +\lambda ^{-3/2}\Theta _{++}^3\\ \Theta _{+-}(\cdot ,\lambda )&=\Theta _{+-}^0 + \lambda ^{-1/2}\Theta _{+-}^1 + \lambda ^{-1}\Theta _{+-}^2 +\lambda ^{-3/2}\Theta _{+-}^3\\ \Theta _{-+}(\cdot ,\lambda )&=\Theta _{-+}^0 + \lambda ^{-1/2}\Theta _{-+}^1 + \lambda ^{-1}\Theta _{-+}^2 +\lambda ^{-3/2}\Theta _{-+}^3\\ \Theta _{--}(\cdot ,\lambda )&=\Theta _{--}^0 + \lambda ^{-1/2}\Theta _{--}^1 + \lambda ^{-1}\Theta _{--}^2 +\lambda ^{-3/2}\Theta _{--}^3, \end{aligned}\right. \end{aligned}$$

with

$$\begin{aligned} \Theta ^0_{++}&= \begin{pmatrix} 0 &{}-\frac{\bar{\tau }^2}{\nu }\\ -\frac{\bar{\alpha }\bar{\tau }^2}{2\nu }&{}\frac{\bar{\tau }'}{\bar{\tau }}-\frac{\bar{c}\bar{\tau }^2}{2\nu }\end{pmatrix}, \qquad \Theta ^1_{+-}= \begin{pmatrix} -\frac{2\bar{\tau }'}{\sqrt{\nu }}+\frac{\bar{\tau }^3}{\nu ^{3/2}}(\frac{\bar{\alpha }}{\bar{c}}+\bar{c})\\ -\frac{\bar{\tau }^{2}\left( q-\bar{c}\bar{\tau }\right) }{\sqrt{\nu }}-\frac{\bar{\alpha }}{2}\frac{\bar{\tau }^3}{\nu ^{3/2}} \end{pmatrix},\\ \Theta ^1_{++}&= \begin{pmatrix} 0&{} -\frac{2\bar{\tau }'}{\sqrt{\nu }}+\frac{\bar{\tau }^3}{\nu ^{3/2}}(\frac{\bar{\alpha }}{\bar{c}}+\bar{c})\\ \frac{\bar{\tau }}{2\sqrt{\nu }}\left( \left( q-\bar{c}\bar{\tau }\right) ^2+c\bar{\alpha }\frac{\bar{\tau }^2}{\nu }\right) &{} \frac{\bar{\tau }^{2}\left( q-\bar{c}\bar{\tau }\right) }{\sqrt{\nu }}+\frac{\bar{\alpha }}{2}\frac{\bar{\tau }^3}{\nu ^{3/2}} \end{pmatrix},\\ \Theta ^2_{++}&=\begin{pmatrix} 0 &{} -\frac{2\bar{\tau }^{3}\left( q-\bar{c}\bar{\tau }\right) }{\nu }\\ 0&{}\frac{\bar{\tau }^{2}\left( q-\bar{c}\bar{\tau }\right) ^2}{2\nu }+\frac{\bar{c}\bar{\alpha }}{2\sqrt{\nu }}\frac{\bar{\tau }^3}{\nu ^{3/2}} \end{pmatrix},\\ \Theta ^3_{++}&=\begin{pmatrix} 0 &{} -\frac{\bar{\tau }^{3}\left( q-\bar{c}\bar{\tau }\right) ^2}{\nu ^{3/2}}-\frac{\bar{c}\bar{\alpha }}{\sqrt{\nu }}\frac{\bar{\tau }^4}{\nu ^2}\\ 0&{}0\end{pmatrix}, \qquad \Theta ^0_{+-}= \begin{pmatrix} \frac{\bar{\tau }^2}{\nu }\\ \frac{\bar{\tau }'}{2\bar{\tau }}-\frac{\bar{c}}{2}\frac{\bar{\tau }^2}{\nu }\end{pmatrix},\\ \Theta ^2_{+-}&= \begin{pmatrix} \frac{2\bar{\tau }^{3}\left( q-\bar{c}\bar{\tau }\right) }{\nu }\\ \frac{\bar{\tau }^{2}\left( q-\bar{c}\bar{\tau }\right) ^2}{2\nu }+\frac{\bar{c}\bar{\alpha }}{2\sqrt{\nu }}\frac{\bar{\tau }^3}{\nu ^{3/2}}\end{pmatrix}, \qquad \Theta ^3_{+-}= \begin{pmatrix} -\frac{\bar{\tau }^{3}\left( q-\bar{c}\bar{\tau }\right) ^2}{\nu ^{3/2}}-\frac{\bar{c}\bar{\alpha }}{\sqrt{\nu }}\frac{\bar{\tau }^4}{\nu ^2}\\ 0\end{pmatrix},\\ \Theta ^0_{-+}&= \begin{pmatrix} \frac{\bar{\alpha }\bar{\tau }^2}{2\nu }&\frac{\bar{\tau }'}{\bar{\tau }}-\frac{\bar{c}\bar{\tau }^2}{2\nu } \end{pmatrix}, \qquad \Theta ^2_{-+}= \begin{pmatrix} 0&\frac{\bar{\tau }^{2}\left( q-\bar{c}\bar{\tau }\right) ^2}{2\nu }+\frac{\bar{c}\bar{\alpha }}{2\sqrt{\nu }}\frac{\bar{\tau }^3}{\nu ^{3/2}} \end{pmatrix},\\ \Theta ^1_{-+}&= \begin{pmatrix} \frac{\bar{\tau }}{2\sqrt{\nu }}\left( \left( q-\bar{c}\bar{\tau }\right) ^2+\bar{c}\bar{\alpha }\frac{\bar{\tau }^2}{\nu }\right)&\frac{\bar{\tau }^{2}\left( q-\bar{c}\bar{\tau }\right) }{\sqrt{\nu }}+\frac{\bar{\alpha }}{2}\frac{\bar{\tau }^3}{\nu ^{3/2}}\end{pmatrix}, \qquad \Theta ^0_{--}=\begin{pmatrix}\frac{\bar{\tau }'}{\bar{\tau }}-\frac{\bar{c}\bar{\tau }^2}{2\nu } \end{pmatrix},\\ \Theta ^1_{--}&= \begin{pmatrix}-\frac{\bar{\tau }^{2}\left( q-\bar{c}\bar{\tau }\right) ^2}{2\nu }+\frac{\bar{c}\bar{\alpha }}{2\sqrt{\nu }}\frac{\bar{\tau }^3}{\nu ^{3/2}}\end{pmatrix}, \qquad \Theta ^2_{--}=\begin{pmatrix}\frac{\bar{\tau }^{2}\left( q-\bar{c}\bar{\tau }\right) ^2}{2\nu }+\frac{\bar{c}\bar{\alpha }}{2\sqrt{\nu }}\frac{\bar{\tau }^3}{\nu ^{3/2}}\end{pmatrix}. \end{aligned}$$

From this, the matrices $N_{D,D},N_{H,D},N_{D,H}$ in (2.17) are obtained through the identification

$$\begin{aligned} \left( \begin{array}{cc} \Theta _{++}&{}\Theta _{+-}\\ \Theta _{-+}&{}\Theta _{--}\\ \end{array} \right) \ =\ \left( \begin{array}{cc} 0 &{}N_{H,D}\\ N_{D,H}&{}N_{D,D}\\ \end{array} \right) , \end{aligned}$$

where $N_{DD}$ is a $2\times 2$ matrix.

Appendix 2: Computational Methods

For completeness, we very briefly describe the computational methods utilized in our investigations reported in Sects. 3.2 and 3.3 above. For more details, the interested reader is referred to Barker et al. (2013) where an analogous numerical study is performed on the generalized Kuramoto–Sivashinsky equation.

1.1 Hill’s Method

To determine the global picture of spectrum of a linear X-periodic operator L, we use Hill’s method. The linear operator L takes the form $L_{j,k}= \sum _{q=1}^{m_{jk}}f_{j,k,q}(x)\frac{\partial ^q}{\partial x^q}$ where the $f_{j,k,q}(x)$ are X periodic. Following Deconinck et al. (2007), we represent the coefficient functions $f_{j,k,q}(x)$ as a Fourier series $f_{j,k,q}(x)=\sum _{j=-\infty }^{\infty } \hat{\phi }_{j,k,q}e^{i2\pi jx/X}$. We use MATLAB’s fast Fourier transform to determine the coefficients $\hat{\phi }_{j,k,q}$. The generalized eigenfunctions are represented as²⁵ $v(x) = e^{i\xi x}\sum _{j=-\infty }^{\infty } \hat{v}_j e^{i\pi jx/X}$, where $\xi \in [-\pi /2X,\pi /2X)$ is the Floquet exponent. Upon substituting the Fourier expansions into the eigenvalue problem, fixing $\xi$, and equating the coefficients of the resulting Fourier series, we arrive at the eigenvalue problem $\hat{L}^{\xi } \hat{v} = \lambda \hat{v}$ where $\hat{L}^{\xi }$ is an infinite-dimensional matrix. The spectrum of the operator L is given by $\sigma (L) = \bigcup _{\xi \in [-\pi /2X,\pi /2X)} \sigma (L_{\xi })$. Truncating the Fourier series at N terms leads to a finite-dimensional eigenvalue problem $L_N^{\xi }\hat{v} = \lambda \hat{v}$. The matrix $L_N^{\xi }$ is of the form $M_2^{-1}M_1$ where $M_1\hat{v} = \lambda M_2\hat{v}$ is the original eigenvalue problem. Typically $M_2$ is the identity, but in (3.14), $M_2$ is diagonal with jth diagonal entry $i(j+\xi )$; hence, we avoid $\xi = 0$ in that case so that $M_2$ is invertible. We compute $\sigma (L_N^{\xi })$ on a mesh to approximate the spectral curves of L. For our numerical studies, we used the implementation of Hill’s method built into STABLAB [BHZ]. For discussion of Hill’s method and its convergence, see Curtis and Deconinck (2010), Deconinck and Nathan Kutz (2006), Johnson and Zumbrun (2012).

1.2 Evans Function

Our results for Hill’s method are augmented by use of the Evans function. To this end, note all the spectral problems we study, such as the one given in (3.14), may be written as a first-order system of the form $W'(x)=\mathbb {A}(x;\lambda )W(x)$ and that, further, $\lambda \in \mathbb {C}$ belongs to the essential spectrum of the associated linearized operator L if and only if this first-order system admits a nontrivial solution satisfying

$$\begin{aligned} Y(x+X;\lambda )=e^{i\xi X}Y(x;\lambda ),\quad \forall x\in \mathbb {R}\end{aligned}$$

where here X denotes the period of the coefficients of L. Following Gardner (1993), the Evans function is defined as $D(\lambda ,\xi ):= \det \left( \Psi (X,\lambda )-e^{i\xi X} \right)$ where the matrix $\Psi (x,\lambda )$ satisfies $\partial _x\Psi (x,\lambda ) = \mathbb {A}\Psi (x,\lambda )$ and $\Psi (0,\lambda ) = \mathrm{Id }$. By construction then, the roots of $D(\cdot ;\xi )$ agree in location and algebraic multiplicity with the eigenvalues of the associated $\xi$-dependent spectral problem. Unfortunately, the Evans function as described here is poorly conditioned for numerical computation. To remedy this, as in Barker et al. (2013), we use the observation of Gardner (1993) that

$$\begin{aligned} D(\lambda ,\xi ):=\det (\Psi (X)-e^{i\xi X}\mathrm{Id })=\det \begin{pmatrix} \Psi (X)&{} e^{i\xi X}\mathrm{Id }\\ \mathrm{Id }&{} \mathrm{Id }\end{pmatrix}, \end{aligned}$$

to express the Evans function as an exterior product of solutions of

$$\begin{aligned} \begin{pmatrix} Y\\ \alpha \end{pmatrix}'= \begin{pmatrix} \mathbb {A}(\cdot ,\lambda )Y \\ 0\end{pmatrix}, \end{aligned}$$

with data $( \mathrm{Id },\mathrm{Id })^T$ at $x=0$ and $( e^{i\xi X}\mathrm{Id },\mathrm{Id })^T$ at $x=X$; for details see Barker et al. (2013). We then use the polar coordinate method of Humpherys and Zumbrun (2006) to evolve the solutions. This algorithm is numerically well conditioned Zumbrun (2009). All computations were carried out using STABLAB [BHZ].

As mentioned above, Hill’s method is ideal for obtaining a global picture of the spectrum and the Evans function can be evaluated on contours and the winding number evaluated to determine the presence of zeros. However, neither method determines definitively whether unstable spectra of arbitrarily small size exist, due to numerical error. In particular, such methods cannot be used to determine the spectrum of the associated linearized operators in a sufficiently small neighborhood of the origin, i.e., they cannot resolve the modulational instability problem. A strategy for rigorously computing stability, which we utilize in our intermediate F numerical studies reported in 3.3 above, involves a Taylor expansion of the Evans function as we now briefly describe; see Barker et al. (2013) for details.

Due to the presence of a conservation law in the governing system (see Serre (2005), Johnson et al. (2011), Johnson et al. (2014), Rodrigues (2013)) the Evans function has a double root at the origin when $\xi = 0$. As such, the Taylor expansion of the Evans function about the origin $(\lambda ,\xi )=(0,0)$ takes the form $D(\lambda ,\xi )=c_{2,0}\lambda ^2+c_{1,1}\lambda \xi +c_{0,2}\xi ^2+c_{3,0}\lambda ^3+c_{2,1}\lambda ^2\xi +c_{1,2}\lambda \xi ^2+c_{0,3}\xi ^3+O(|\lambda |^4+|\xi |^4)$ where the coefficients $c_{k,j}$ may be determined via Cauchy’s integral formula,

$$\begin{aligned} c_{k,j}=-\frac{1}{4\pi ^2}\oint _{\partial B(0,r)}\oint _{\partial B(0,r)}D(\lambda ,\xi )\lambda ^{-k-1}\xi ^{-j-1}d\lambda ~d\xi \end{aligned}$$

(3.16)

with $r>0$ sufficiently small. Setting $\alpha _j=\frac{-c_{1,1}+(-1)^{j+1}\sqrt{c_{1,1}^2-4c_{2,0}c_{0,2}}}{2c_{2,0}}$, $\beta _j=-\frac{c_{3,0}\alpha _j^3+c_{2,1}\alpha _j^2 +c_{1,2}\alpha _j+c_{0,3}}{2c_{2,0}\alpha _j+c_{1,1}}$, one readily checks that the roots of the Evans function near $(\lambda ,\xi )=(0,0)$ may be continued for $|\xi |\ll 1$ as

$$\begin{aligned} \lambda _j(\xi )= \alpha _j \xi +\beta _j\xi ^2+\frac{\xi ^3}{2}\int _0^1(1-s)^2\lambda _j'''(s\xi )ds. \end{aligned}$$

The spectral curves at the origin are thus approximated by $\alpha _j \xi + \beta _j \xi ^2$ with spectral stability corresponding to the case $\alpha _j\in \mathbb {R}i$ and $\mathfrak {R}(\beta _j)<0$; see Barker et al. (2013) for details.

In practice, to compute the Taylor expansion coefficients, rather than to compute the Evans function on the contour integral in the variable $\lambda$ for fixed $\xi$ given in (3.16), we interpolate with $\sum _{k=0}^K e^{i k\xi x}$ ($K =3$ is the largest power of $e^{i\xi x}$ that appears) and then use the Taylor expansion $e^{ik\xi x} = 1+ (ik\xi x)+(ik\xi x)^2/2 + \cdot$ yielding $D(\lambda ,\xi ) = \sum _{k=0}^{\infty } c_k \xi ^k$, from which the contour integral can be determined simply by reading off the corresponding coefficient. Calling the quantity just determined $\tilde{D}$, we see

$$\begin{aligned} \frac{1}{2\pi i} \oint _{|\lambda |= R_1} \frac{\tilde{D}(\lambda )}{\lambda ^{r+1}} d\lambda = \frac{1}{2\pi } \int _{-\pi }^{\pi } \frac{\tilde{D}(Re^{i\theta })}{R^re^{ir\theta }}d\theta = \frac{1}{2R^r}\int _{-1}^{1} \tilde{D}(Re^{i\pi \theta }) e^{-ir\pi \theta } d\theta , \end{aligned}$$

which we compute by approximating the integrand with Chebyshev interpolation and integrating.

Computational Effort

1.1 Computational Environment

In carrying out our numerical investigations, we used a MacBook laptop with 2GB memory and a duo core Intel processor with 2GHz processing speed, a 2009 Mac Pro with 16GB memory and two quad-core Intel processors with 2.26 GHz processing speed, and Quarry, a supercomputer at Indiana University consisting of 140 IBM HS21 Blade servers with two 2GH quad-core Intel Zeon 5335 processors per node and delivering 8.96 teraflops processing speed. All computations were done using MATLAB and the MATLAB-based stability package STABLAB.

1.2 Computational Time

We begin by providing computational statistics for the representative parameter set $\alpha = -2$, $q_0 = 0.4$, $F = 10$, and $X = 50$. We compute the Evans function, $D(\lambda ,\xi )$, on a semicircular contour, $\varOmega$, of radius $R = 0.2$ with 42 evenly spaced Floquet parameters $\xi \in [-\pi /X,-\pi /(10X)]\cup [\pi /(10X),\pi /X]$. We require the relative error between contour points of the image contour, $Y_{\xi }$, $D(\cdot ,\xi ): \varOmega \rightarrow Y_{\xi }$, not exceed 0.2 so that Rouché’s theorem implies the winding number of $Y_{\xi }$ corresponds to the number of roots of $D(\cdot ,\xi )$ in $\varOmega$. We use 277 points in $\varOmega$, chosen adaptively, to achieve 0.2 relative error in each $Y_{\xi }$ at a computational cost of 56.0 s using 8 MATLAB workers on Quarry to determine the winding number is zero. Computing the Taylor expansion of the Evans function at the origin requires 61.7 s on Quarry, while computing the spectrum via Hill’s method using 603 Fourier modes and 21 Floquet parameters comes at a computational cost of 143 s.

As the period X increases or as F increases, the number of Fourier modes needed in Hill’s method increases. Using 600 Fourier modes typically takes around 3 min on the Mac Pro, while using 1600 Fourier modes takes about 77 min, and using 3000 Fourier modes requires approximately 8 h.

In creating Fig. 4a, it took 4.36 days of computation time to evaluate the Taylor coefficients and 34.5 days to compute the spectrum using Hill’s method, while the Evans function required an estimated 58 h. A typical profile requires only a few seconds to compute, but we must use continuation whereby we use a nearby profile as an initial guess in the boundary value solver, so that computing the profiles also required a great computational effort. Overall, taking into account the use of parallel computing and all values of $\alpha$ investigated, we estimate that total computational time for the project exceeds a year.