Linear stability of plane Poiseuille flow in an infinite elastic medium and volcanic tremors

Abstract

The linear stability of a plane compressible laminar (Poiseuille) flow sandwiched between two semi-infinite elastic media was investigated with the aim of explaining the excitation of volcanic tremors. Our results show that there are several regimes of instability, and the nature of stability significantly depends on the symmetry of oscillatory fluid and solid motion. It has been shown that long-wave symmetric instability occurs at a very small value of the Reynolds number, but it is unlikely that this is the cause of volcanic tremors. We show that antisymmetric (flexural) instability also occurs, involving two parallel Rayleigh waves traveling against the Poiseuille flow, but the critical flow speed is faster than that of symmetric instability. However, if the basic flow profile is nonparabolic because of a nonuniform driving force or nonuniform viscosity, the critical flow speed of antisymmetric instability can be considerably slower than that of symmetric instability. Based on numerical calculations and analytical consideration, we conclude that this anomalous antisymmetric instability is possibly produced by a basaltic magma flow of a few meters per second through a dike with thickness of 1 m and extending for several kilometers; this origin can explain some of the characteristics of volcanic tremors.

Background

Volcanic tremors are characterized as long-period and long-lasting seismic events that are uniquely observed at active volcanoes (for reviews, see Chouet1996; Konstantinou and Schlindwein2002; McNutt2005). The frequency spectra of tremor oscillations have one or several distinct peaks, which in some cases of harmonic tremors are recognized as a main oscillation and its overtones. The characteristic periods often change over time with change in volcanic activity, but mostly lie in the range of 0.2 to 2 s, despite the fact that the source depth, style of eruption, and geological structure vary between volcanoes. Some tremors are known to last for several days, exhibiting a striking contrast to long-period volcanic earthquakes that share spectral content that is similar to those of tremor events but monotonically decay within a few minutes. However, it is often regarded that these two phenomena involve a common source process, at least in part, which is very different from the generation mechanisms of tectonic earthquakes.

Theoretical considerations and observational evidence have led to a general belief that the existence of an underground fluid plays an important role in the generation of volcanic tremors, but the detailed mechanisms underlying tremors have long been in debate (e.g., Aki et al.1977; Ferrick et al.1982; Chouet1988; Iwamura and Kaneshima2005; Jellinek and Bercovici2011). Peaked spectra are suggestive of the characteristic oscillations of certain structures buried under volcanoes, which may also explain the long-period nature if the slow sound speed of volcanic fluids has some relevance to the oscillations. A mutual relation between the occurrence of volcanic tremors and eruptions offers some grounds for considering that the process of magma transport plays a key role. In this paper, we discuss a wave generation process in which a fluid flows in an underground conduit and oscillations are self-excited because of fluid dynamical instability aided by elastic deformation of the surrounding rock. This mechanism is analogous to flow-induced vibration such as the aeroelastic flutter of airfoils and the vibration of reeds in some musical instruments, but the model we present differs from these practical applications in that it involves oscillations of an infinitely extended homogeneous elastic body. In active volcanoes, magmas, waters, and volcanic gases are ubiquitous, and we hypothesize that the transport of fluid through cracks or tube-like conduits is common. It is possible that the flow-induced vibration may account for the long-lasting feature of volcanic tremors, because the oscillations last as long as the flow continues.

Julian (1994) undertook a pioneering study on the possibility of flow-induced volcanic tremors, wherein a nonlinear lumped-parameter model was considered, in which a fluid rapidly flowed between two rigid walls, and a spring and a dashpot connected each rigid wall to a stationary background basis. This model produces an oscillatory instability with a peaked frequency spectrum after the flow speed reaches a certain threshold. The spectrum subsequently experiences period-doubling bifurcations with an increase in the flow speed and finally becomes chaotic. Further fundamental physics relevant to the tremor process can also be shown by analyzing a fluid flow in an infinite elastic continuum. Balmforth et al. (2005) considered instability in a plane viscous laminar (Poiseuille) flow embedded between two semi-infinite elastic media. In their linear stability analysis, they showed that the critical flow speed at which the laminar flow is destabilized approaches zero as the wavelength of the elastic deformation lengthens considerably. The instability accompanies a wave traveling at a speed comparable to the main laminar flow. This wave can be identified as a slow surface wave concentrated around the fluid layer (Rust et al.2008), known as a Krauklis wave or a crack wave in volcano seismology (Krauklis1962; Chouet1986; Korneev2011). However, in possible volcanic settings, this flow-induced instability appears to be unrealistic because it still requires a too-rapid fluid flow, and the wave amplification may not be efficient (Rust et al.2008). In Balmforth et al. (2005), the fluid is assumed to be incompressible and both the fluid and solid motions are assumed to have a wavelength that is considerably longer than the fluid layer thickness. Therefore, the full equation of the fluid motion is not treated, but either a boundary-layer model is applied to solve the Navier-Stokes equation or a more simplified slot-averaging model is used, in which the flow perturbation is uniform in the direction perpendicular to the main flow. Dunham and Ogden (2012) recently performed a similar analysis assuming that the fluid is compressible, but this study used only a slot-averaging model.

Here, we reconsider the linear stability of a plane Poiseuille flow surrounded by an infinite elastic medium. Our motivation for doing so is that the physical nature of this system as a wave generator is not fully understood because of the fact that previous studies have made considerable assumptions. A plane-layer model is not only physically simple and fundamental, but such a sheet-like structure, known geologically as a dike, is considered to be commonly distributed under volcanoes. Magma transport through a dike is one of the most fundamental processes in magmatism (e.g., Best2003). For example, a basaltic magma that has relatively small viscosity tends to ascend from a magma reservoir creating a thin planar crack in a largely vertical direction, driven by an increase in pressure in the reservoir and the buoyancy acting on the fluid magma, and magma rises to the surface in the unique case of fissure eruption. The relationship between magma-filled dikes and tremors (or long-period earthquakes) has been pointed out in previous studies (e.g., Aki et al.1977; Chouet1988; Furumoto et al.1990; Neuberg et al.2006), although these studies were not necessarily concerned with the aspect of flow-induced vibration.

If the conduit wall is rigid and the fluid incompressible, the Orr-Sommerfeld equation describes the linear stability and the resulting flow structure of the unstable flow (e.g., Drazin and Reid2004). We extend this equation so that it is applicable to a compressible fluid flowing in an infinite elastic continuum. We make no assumption that the wavelength is considerably longer than the conduit thickness or that the flow profile is uniform in the direction perpendicular to the main flow. Compressibility might have only minor effects, but we retain this because the dispersion relation of the Krauklis wave generally depends on the ratio between the elastic wave speeds of a fluid and a solid. In general, because of the symmetry of the system, the fluid and solid motion can be represented as a superposition of symmetric and antisymmetric solutions (Figure1). In the symmetric solution, the fluid-filled crack opens and closes symmetrically with a wave-type pressure perturbation, but the antisymmetric solution exhibits a flexural deformation of the fluid layer with a relatively constant thickness. It is of considerable interest that previous studies have focused only on the symmetric solution, and this is probably because the symmetric Krauklis wave exhibits the interesting characteristic of a phase velocity that approaches zero when the wavelength becomes very long (Krauklis1962; Ferrazzini and Aki1987), which has been thought to possibly explain the long-period nature of volcanic tremors. In fact, Balmforth et al. (2005) and Dunham and Ogden (2012) found that the critical flow speed for the onset of instability could approach zero in a long-wave symmetric solution. However, it is well known that antisymmetric disturbances are first destabilized in the Orr-Sommerfeld problem with rigid plane boundaries when the flow speed is increased. Therefore, in our analysis, we also focus on the antisymmetric instability.

Figure 1

Illustration of our model. A plane fluid layer of a constant thickness, 2a, is sandwiched between two semi-infinite elastically deformable bodies. The fluid inside flows in one direction. Deformation of the elastic surrounding and the fluid motion can be generally represented by the superposition of symmetric and antisymmetric solutions.

Numerical calculations were performed for linearly unstable wave solutions. The equations and boundary conditions are described in Section ‘Methods,’ and the results are presented in Section ‘Results’ with some physical interpretations of instability. In Section ‘Discussion,’ we attempt to apply our model to the excitation of volcanic tremors. Subsequently, we discuss the flow of magma through a thin sheet-like dike, and we conclude that instability in a classical Poiseuille flow with a parabolic velocity profile cannot explain volcanic tremors, but when the magma flow has a velocity profile that is slightly deviated from the parabolic one, a moderate flow speed can cause instability with antisymmetric fluid and solid motion, which may be a source of some volcanic tremors.

Methods

Fundamental equations

We consider a compressible barotropic fluid of constant viscosity, ζ^∗, and sound velocity,${\alpha }_0^{\ast }$, sandwiched between two semi-infinite elastic media, |x^∗| ≥ a^∗, in which the density is uniformly${\rho }_s^{\ast }$ and the P- and S-wave velocities are respectively${\alpha }_{\mathrm{s}}^{\ast }=\sqrt{\left({\lambda }_{\mathrm{s}}^{\ast }+2{\mu }_{\mathrm{s}}^{\ast }\right)/{\rho }_{\mathrm{s}}^{\ast }}$ and${\beta }_{\mathrm{s}}^{\ast }=\sqrt{{\mu }_{\mathrm{s}}^{\ast }/{\rho }_{\mathrm{s}}^{\ast }}$ (${\lambda }_{\mathrm{s}}^{\ast }$ and${\mu }_{\mathrm{s}}^{\ast }$ represent the Lamé constants). The quantity with an asterisk indicates that it has a dimension. Because we focus on a two-dimensional fluid motion that is largely along the z^∗-direction, we neglect the y^∗-dependence of any physical quantities. In the fluid layer, we employ the Navier-Stokes equation and the equation of continuity:

$$\begin{array}{ll}{\rho }^{\ast }\left[\frac{\partial {\mathit{V}}^{\ast }}{\partial t^{\ast }}+\left({\mathit{V}}^{\ast }·{\mathit{\nabla }}^{\ast }\right){\mathit{V}}^{\ast }\right]& =-{\mathit{\nabla }}^{\ast }{P}^{\ast }+{\zeta }^{\ast }\left[{\nabla {\ast }^{}}^2{\mathit{V}}^{\ast }\right.\\ +\left.\frac{1}{3}{\mathit{\nabla }}^{\ast }\left({\mathit{\nabla }}^{\ast }·{\mathit{V}}^{\ast }\right)\right]+{\mathit{K}}^{\ast },\end{array}$$

(1)

$$\begin{array}{l}\frac{\partial {\rho }^{\ast }}{\partial t^{\ast }}+{\mathit{V}}^{\ast }·{\mathit{\nabla }}^{\ast }{\rho }^{\ast }=-{\rho }^{\ast }{\mathit{\nabla }}^{\ast }·{\mathit{V}}^{\ast },\end{array}$$

(2)

where${\mathit{V}}^{\ast }=V_x^{\ast }\left(x^{\ast },z^{\ast },t^{\ast }\right){\mathit{e}}_x+V_z^{\ast }\left(x^{\ast },z^{\ast },t^{\ast }\right){\mathit{e}}_z$ denotes the flow velocity; P^∗(x^∗,z^∗,t^∗), the pressure; ρ^∗(x^∗,z^∗,t^∗), the density; and K^∗ represents an external body force. The displacement vector,${\mathit{U}}^{\ast }=U_x^{\ast }\left(x^{\ast },z^{\ast },t^{\ast }\right){\mathit{e}}_x+U_z^{\ast }\left(x^{\ast },z^{\ast },t^{\ast }\right){\mathit{e}}_z$, in the elastic medium satisfies the equation of motion:

$$\frac{{\partial }^2{\mathit{U}}^{\ast }}{\partial {t^{\ast }}^2}={{\alpha }_s^{\ast }}^2{\mathit{\nabla }}^{\ast }\left({\mathit{\nabla }}^{\ast }·{\mathit{U}}^{\ast }\right)-{{\beta }_{\mathrm{s}}^{\ast }}^2{\mathit{\nabla }}^{\ast }\times {\mathit{\nabla }}^{\ast }\times {\mathit{U}}^{\ast }.$$

(3)

The fluid and solid motions are connected by boundary conditions at the fluid-solid interface. We adopt the no-slip condition and continuity of traction across the fluid-solid boundary, the latter of which is derived from the relation between the stress components in the fluid:

$$\left\{\begin{array}{l}{T}_{\mathit{\text{xx}}}^{\ast }=-P^{\ast }+2{\zeta }^{\ast }\left(\frac{\partial V_x^{\ast }}{\partial x^{\ast }}-\frac{1}{3}{\mathit{\nabla }}^{\ast }·{\mathit{V}}^{\ast }\right),\\ T_{\mathit{\text{zx}}}^{\ast }={\zeta }^{\ast }\left(\frac{\partial V_z^{\ast }}{\partial x^{\ast }}+\frac{\partial V_x^{\ast }}{\partial z^{\ast }}\right),\end{array}\right.$$

(4)

and those in the solid:

$$\left\{\begin{array}{l}{\mathrm{\Sigma }}_{\mathit{\text{xx}}}^{\ast }=\left({\lambda }_{\mathrm{s}}^{\ast }+2{\mu }_{\mathrm{s}}^{\ast }\right){\mathit{\nabla }}^{\ast }·{\mathit{U}}^{\ast }-2{\mu }_{\mathrm{s}}^{\ast }\frac{\partial U_z^{\ast }}{\partial z^{\ast }},\\ {\mathrm{\Sigma }}_{\mathit{\text{zx}}}^{\ast }={\mu }_{\mathrm{s}}^{\ast }\left(\frac{\partial U_z^{\ast }}{\partial x^{\ast }}+\frac{\partial U_x^{\ast }}{\partial z^{\ast }}\right).\end{array}\right.$$

(5)

Basic state

We assume a basic state in which the fluid density is uniformly${\rho }_0^{\ast }$ and the flow is driven by a body force${\mathit{K}}^{\ast }={\bar{K}}^{\ast }\left(x^{\ast }\right){\mathit{e}}_z$ exerted along the z^∗-direction. When${\bar{K}}^{\ast }$ is uniform, the balance between the driving force and the viscous force leads to a solution known as the Poiseuille flow:

$${\mathit{V}}^{\ast }={\bar{V}}^{\ast }\left(x^{\ast }\right){\mathit{e}}_z=V_0^{\ast }\left(1-\frac{{x^{\ast }}^2}{{a^{\ast }}^2}\right){\mathit{e}}_z,$$

(6)

where$V_0^{\ast }$ denotes the flow speed at x^∗ = 0. The basic flow involves a surface force acting on the fluid-solid boundary, and this should be balanced by the elastic force due to deformation of the surrounding solid. However, it is intrinsically impossible for the elastic deformation to sustain the static surface force because there is no fixed boundary in the solid medium. The remedy we adopt is to assume that a no-deformation state, U^∗ = 0, is also sustained by a prescribed external surface force acting on the bottom of the solid medium. We assume that${\bar{\mathrm{\Sigma }}}_{\mathit{\text{zx}}}^{\ast }\left(x^{\ast }\right)$ in the basic state is entirely zero at |x^∗| > a^∗, but it coincides with${\bar{T}}_{\mathit{\text{zx}}}^{\ast }={\zeta }^{\ast }\mathrm{d}{\bar{V}}^{\ast }/\mathrm{d}{x}^{\ast }$ at x^∗ = ± a^∗. This assumption makes the fluid layer thickness in the basic state uniformly 2a^∗, and the displacement vector can be treated as a small-amplitude variable in the subsequent linear analysis. The problem with use of the basic state has been recognized in previous studies (Balmforth et al.2005), and an artificial approximation seems to be inevitable when proceeding to perform a simplified wave analysis.

In Section ‘Results’ we show that a slight change in the functional form of the driving force brings about a drastic change in the nature of the stability. Here, we consider a generalized model for the basic flow:

$${\bar{V}}^{\ast }\left(x^{\ast }\right)=V_0^{\ast }\left(1-\frac{{x^{\ast }}^2}{{a^{\ast }}^2}\right)\left[1-\frac{\delta {x^{\ast }}^2}{(6-\delta ){a^{\ast }}^2}\right],$$

(7)

which is driven by the body force:

$${\bar{K}}^{\ast }\left(x^{\ast }\right)=-\frac{\mathrm{d}{\bar{T}}_{\mathit{\text{zx}}}^{\ast }}{\mathrm{d}{x}^{\ast }}=\frac{12{\zeta }^{\ast }{V}_0^{\ast }}{(6-\delta ){a^{\ast }}^2}\left(1-\frac{\delta {x^{\ast }}^2}{{a^{\ast }}^2}\right).$$

(8)

The parameter δ specifies the degree of the deviation of the body force from the constant profile, and the classical Poiseuille flow (6) with a parabolic velocity profile can be represented as a special case with δ = 0. When δ = 1, the external body force vanishes at the fluid-solid boundary. The deviation from the classical Poiseuille flow occurs not only when the driving force is nonuniform, but when the fluid viscosity depends on x^∗. In a magma conduit, the latter situation is highly probable because the viscosity of the cooler magma near the bedrock is likely to be higher than that of the internal hot magma. However, in this paper, we simplify the model using the flow (7) with a constant viscosity, and we investigate the effect of the deviation from the Poiseuille flow (6) by changing the flow parameter δ from zero to around unity.

Linearized equations

The governing equations (1) to (3) are linearized about the basic state. The density perturbation,${\rho }^{\ast }\left(x^{\ast },z^{\ast },t^{\ast }\right)-{\rho }_0^{\ast }$, relates to the pressure perturbation, p^∗(x^∗,z^∗,t^∗), via the equation of state:

$$p^{\ast }={{\alpha }_0^{\ast }}^2\left({\rho }^{\ast }-{\rho }_0^{\ast }\right).$$

(9)

In this and the following sections, we use lower-case letters to represent small-amplitude perturbation variables, and we mainly discuss the fluid and solid motion with nondimensional variables, taking${\alpha }_0^{\ast }$, a^∗,$a^{\ast }/{\alpha }_0^{\ast }$, and${\rho }_0^{\ast }{{\alpha }_0^{\ast }}^2$ as the units of velocity, length, time, and pressure, respectively. Therefore, the nondimensional fluid layer thickness is 2, the nondimensional perturbation velocity is$\mathit{v}(x,z,t)=\mathit{V}(x,z,t)-\bar{V}(x){\mathit{e}}_z$, and the basic flow speed (7) is reduced to

$$\bar{V}(x)=M\left(1-x^2\right)\left(1-\frac{\delta x^2}{6-\delta }\right),$$

(10)

where

$$M=\frac{V_0^{\ast }}{{\alpha }_0^{\ast }}$$

(11)

denotes the Mach number.

As we have neglected the elastic deformation in the basic state, the displacement is simply treated as a small-amplitude variable and it can generally be expressed as u = ∇ f+∇ × (g e _y ). The governing equations (1) to (3) are reduced to nondimensional and linearized forms:

$$\begin{array}{l}{S}^{-1}\left(\frac{\partial v_x}{\partial t}+\bar{V}\frac{\partial v_x}{\partial z}+\frac{\partial p}{\partial x}\right)={\nabla }^2{v}_x+\frac{1}{3}\frac{\partial }{\partial x}(\mathit{\nabla }·\mathit{v}),\end{array}$$

(12)

$$\begin{array}{l}{S}^{-1}\left(\frac{\partial v_z}{\partial t}+\bar{V}\frac{\partial v_z}{\partial z}+\frac{\mathrm{d}\bar{V}}{\mathrm{d}x}{v}_x+\frac{\partial p}{\partial z}\right)={\nabla }^2{v}_z+\frac{1}{3}\frac{\partial }{\partial z}(\mathit{\nabla }·\mathit{v}),\end{array}$$

(13)

$$\begin{array}{l}\frac{\partial p}{\partial t}+\bar{V}\frac{\partial p}{\partial z}=-\mathit{\nabla }·\mathit{v},\end{array}$$

(14)

$$\begin{array}{l}\frac{{\partial }^{2}f}{\partial t^2}={\alpha }^2{\nabla }^{2}f,\end{array}$$

(15)

$$\begin{array}{l}\frac{{\partial }^{2}g}{\partial t^2}={\beta }^2{\nabla }^{2}g,\end{array}$$

(16)

where

$$S=\frac{{\zeta }^{\ast }}{{\rho }_0^{\ast }{\alpha }_0^{\ast }{a}^{\ast }},\alpha =\frac{{\alpha }_{\mathrm{s}}^{\ast }}{{\alpha }_0^{\ast }},\beta =\frac{{\beta }_s^{\ast }}{{\alpha }_0^{\ast }}.$$

(17)

In this paper, the first nondimensional parameter, S, in (17) is called the fluid viscosity parameter. The Reynolds number can be expressed as R = S^-1M. When the working fluid is a magma of${\zeta }^{\ast }\sim 10^4\text{Pa}·\mathrm{s},{\rho }_0^{\ast }\sim 3\times 10^3$ kg/m ³, and${\alpha }_0^{\ast }\sim 10^3$ m/s, the viscosity parameter is smaller than unity, unless the fluid layer thickness is less than 1 cm. If the working fluid was a less viscous fluid such as high-temperature mafic magma and water, S would be considerably smaller. Therefore, in this paper, we are mainly concerned with cases in which S is smaller than unity.

Boundary conditions

Let us suppose that a particle is located at the fluid-solid boundary (±1, z) in the basic state and it moves to (X,Z) at a certain time t. It follows by definition that

$$\left\{\begin{array}{l}X(z,t)=\pm 1+u_x(\pm 1,z,t),\\ Z(z,t)=z+u_z(\pm 1,z,t).\end{array}\right.$$

(18)

The no-slip condition can be expressed as follows:

$$\left\{\begin{array}{l}{V}_x(X,Z,t)=\frac{\partial X}{\partial t}=\frac{\partial u_x}{\partial t},\\ V_z(X,Z,t)=\frac{\partial Z}{\partial t}=\frac{\partial u_z}{\partial t}.\end{array}\right.$$

(19)

The assumption that the perturbation variables are of a small amplitude gives the following relation:

$$\left\{\begin{array}{ll}{V}_x(X,Z,t)& \simeq V_x(\pm 1,z,t)+\frac{\partial V_x}{\partial x}{u}_x(\pm 1,z,t)+\frac{\partial V_x}{\partial z}{u}_z(\pm 1,z,t)\\ \simeq v_x(\pm 1,z,t),\\ V_z(X,Z,t)& \simeq V_z(\pm 1,z,t)+\frac{\partial V_z}{\partial x}{u}_x(\pm 1,z,t)+\frac{\partial V_z}{\partial z}{u}_z(\pm 1,z,t)\\ \simeq v_z(\pm 1,z,t)+\frac{\mathrm{d}\bar{V}}{\mathrm{d}x}{u}_x(\pm 1,z,t),\end{array}\right.$$

(20)

and therefore, the no-slip condition at x = ±1 is reduced to the following:

$$\begin{array}{l}{v}_x=\frac{\partial u_x}{\partial t},\end{array}$$

(21)

$$\begin{array}{l}{v}_z=\frac{\partial u_z}{\partial t}-\frac{\mathrm{d}\bar{V}}{\mathrm{d}x}{u}_x.\end{array}$$

(22)

In evaluating the traction across the fluid-solid boundary, we need to consider the existence of the initial stress field. The fluid stress across the deformed boundary can be evaluated in terms of the quantities defined in the coordinate system fixed in space in the following manner (see, for example, Eq. (3.36) in Dahlen and Tromp1998):

$$\begin{array}{l}{f}_n={\tau }_{\mathit{\text{xx}}}-{\bar{T}}_{\mathit{\text{zx}}}\frac{\partial u_x}{\partial z},\end{array}$$

(23)

$$\begin{array}{l}{f}_t={\bar{T}}_{\mathit{\text{zx}}}+{\tau }_{\mathit{\text{zx}}}+{\bar{T}}_{\mathit{\text{zx}}}\frac{\partial u_z}{\partial z}+\frac{\mathrm{d}{\bar{T}}_{\mathit{\text{zx}}}}{\mathrm{d}x}{u}_x,\end{array}$$

(24)

where f _n and f _t respectively denote the normal and tangential forces on the deformed boundary per unit undeformed area, and (τ _xx ,τ _zx ) denotes the perturbed Cauchy stress in the fluid defined in a manner similar to (4). The second term in the right-hand side of (23) arises because of a small-angle rotation of the boundary. The third and fourth terms in the right-hand side of (24) respectively represent the change in the unit area and the displacement of the boundary along the x-direction. In the solid side, these forces can be simply expressed as f _n  = σ _xx and$f_t={\bar{\mathrm{\Sigma }}}_{\mathit{\text{zx}}}+{\sigma }_{\mathit{\text{zx}}}$, where (σ _xx ,σ _zx ) is, strictly speaking, the first Piola-Kirchhoff stress due to the small-amplitude elastic deformation and can be defined in a manner similar to (5). Because we assumed${\bar{T}}_{\mathit{\text{zx}}}={\bar{\mathrm{\Sigma }}}_{\mathit{\text{zx}}}$ at x = ±1 and${\bar{T}}_{\mathit{\text{zx}}}=S(\mathrm{d}\bar{V}/\mathrm{d}x)$, the linearized stress boundary condition at x = ±1 can be written as follows:

$$\begin{array}{ll}-p+2S\left(\frac{\partial v_x}{\partial x}-\frac{1}{3}\mathit{\nabla }·\mathit{v}\right)& =\varrho \left({\alpha }^2\mathit{\nabla }·\mathit{u}-2{\beta }^2\frac{\partial u_z}{\partial z}\right)\\ +S\frac{\mathrm{d}\bar{V}}{\mathrm{d}x}\frac{\partial u_x}{\partial z},\end{array}$$

(25)

$$\begin{array}{ll}S\left(\frac{\partial v_z}{\partial x}+\frac{\partial v_x}{\partial z}\right)& =\varrho {\beta }^2\left(\frac{\partial u_z}{\partial x}+\frac{\partial u_x}{\partial z}\right)\\ -S\left(\frac{\mathrm{d}\bar{V}}{\mathrm{d}x}\frac{\partial u_z}{\partial z}+\frac{{\mathrm{d}}^2\bar{V}}{\mathrm{d}{x}^2}{u}_x\right),\end{array}$$

(26)

where$\varrho ={\rho }_{\mathrm{s}}^{\ast }/{\rho }_0^{\ast }$ denotes the density contrast. In the above condition, the terms proportional to S appear to be of secondary importance because the viscosity parameter is generally small. In fact, Balmforth et al. (2005) derived a general boundary condition similar to (25) and (26), but ignored the terms relevant to viscous stress in their linear analysis. However, as shown in Section ‘Antisymmetric solution,’ the last term on the right-hand side of (26) has a significant effect on the flow stability because it is the only term proportional to the displacement, not to the gradient of the displacement, in (25) and (26). Even when S is small, this term may be significant in a long wavelength solution.

Wave analysis

We focus on a wave-type solution wherein the wavenumber is denoted by k (>0), the wavelength is given by L = 2π/k, and the complex phase velocity is denoted as c. Adopting the expressions

$$\left\{\begin{array}{l}{v}_x(x,z,t)=i{\widehat{v}}_x(x)e^{\mathit{\text{ik}}(z-\mathit{\text{ct}})}+\mathit{\text{c.c.}},\\ v_z(x,z,t)={\widehat{v}}_z(x)e^{\mathit{\text{ik}}(z-\mathit{\text{ct}})}+\mathit{\text{c.c.}},\\ p(x,z,t)=\widehat{p}(x)e^{\mathit{\text{ik}}(z-\mathit{\text{ct}})}+\mathit{\text{c.c.}},\\ f(x,z,t)=\widehat{f}(x)e^{\mathit{\text{ik}}(z-\mathit{\text{ct}})}+\mathit{\text{c.c.}},\\ g(x,z,t)=iĝ(x)e^{\mathit{\text{ik}}(z-\mathit{\text{ct}})}+\mathit{c.c.},\end{array}\right.$$

(27)

we obtain the ordinary differential equations:

$$\begin{array}{l}i{S}^{-1}\left[k(\bar{V}-c){\widehat{v}}_x-\frac{\mathrm{d}\widehat{p}}{\mathrm{d}x}\right]=\left(\frac{4}{3}\frac{{\mathrm{d}}^2}{\mathrm{d}{x}^2}-k^2\right){\widehat{v}}_x+\frac{k}{3}\frac{\mathrm{d}{\widehat{v}}_z}{\mathrm{d}x},\end{array}$$

(28)

$$\begin{array}{l}i{S}^{-1}\left[k(\bar{V}-c){\widehat{v}}_z+\frac{\mathrm{d}\bar{V}}{\mathrm{d}x}{\widehat{v}}_x+k\widehat{p}\right]=\left(\frac{{\mathrm{d}}^2}{\mathrm{d}{x}^2}-\frac{4k^2}{3}\right){\widehat{v}}_z-\frac{k}{3}\frac{\mathrm{d}{\widehat{v}}_x}{\mathrm{d}x},\end{array}$$

(29)

$$\begin{array}{l}k(\bar{V}-c)\widehat{p}=-\frac{\mathrm{d}{\widehat{v}}_x}{\mathrm{d}x}-k{\widehat{v}}_z.\end{array}$$

(30)

for the fluid motion, where$i=\sqrt{-1}$, and c.c. denotes the complex conjugate of the preceding term. We assume that the phase velocity, Re(c), is less than the P- and S-wave speeds in the solid. Therefore, the solid motion is expressed as an elastic surface wave concentrated around the fluid layer. In the range of x ≥ 1, for example, the wave equations (15) and (16) permit the general solution:

$$\left\{\begin{array}{l}\widehat{f}(x)=Aexp\left[-k(x-1)\sqrt{1-\frac{c^2}{{\alpha }^2}}\right],\\ ĝ(x)=Bexp\left[-k(x-1)\sqrt{1-\frac{c^2}{{\beta }^2}}\right],\end{array}\right.$$

(31)

where A and B denote complex constants. The boundary conditions, (21), (22), (25), and (26), at x = 1 are subsequently reduced to the following:

$${\widehat{v}}_x(1)=k^{2}c\left(\sqrt{1-\frac{c^2}{{\alpha }^2}}A-B\right),$$

(32)

$${\widehat{v}}_z(1)=k^{2}c\left(A-\sqrt{1-\frac{c^2}{{\beta }^2}}B\right)+k{\bar{V}}^{\prime }(1)\left(\sqrt{1-\frac{c^2}{{\alpha }^2}}A-B\right),$$

(33)

$$\begin{array}{ll}-\widehat{p}(1)+2\mathit{\text{iS}}\left({\widehat{v}}_x^{\prime }(1)-\frac{\mathit{\text{kc}}}{3}\widehat{p}(1)\right)& =2k^2\varrho {\beta }^2\left[\left(1-\frac{c^2}{2{\beta }^2}\right)A\right.\\ -\left.\sqrt{1-\frac{c^2}{{\beta }^2}}B\right]-\mathit{\text{iS}}{k}^2{\bar{V}}^{\prime }(1)\\ \times \left[\sqrt{1-\frac{c^2}{{\alpha }^2}}A-B\right],\end{array}$$

(34)

$$\begin{array}{ll}\mathit{\text{iS}}\left({\widehat{v}}_z^{\prime }(1)-k{\widehat{v}}_x(1)\right)& =2k^2\varrho {\beta }^2\left[\sqrt{1-\frac{c^2}{{\alpha }^2}}A-\left(1-\frac{c^2}{2{\beta }^2}\right)B\right]\\ +\mathit{\text{iSk}}\left\{k{\bar{V}}^{\prime }(1)\left[A-\sqrt{1-\frac{c^2}{{\beta }^2}}B\right]\right.\\ +{\bar{V}}^{\prime \prime }(1)\left.\left[\sqrt{1-\frac{c^2}{{\alpha }^2}}A-B\right]\right\},\end{array}$$

(35)

where the prime denotes the x-derivative.

Numerical method

In the case of the classical Orr-Sommerfeld problem with rigid boundaries, the complex phase velocity, c, can be obtained by solving an eigenvalue problem, provided that k and other nondimensional parameters are given. The system is linearly unstable if the imaginary part, Im(c), of the eigenvalue is positive. However, in our problem, the calculation of c is not straightforward because the boundary conditions depend on c. Here, we adopt a shooting method that is widely used to solve boundary value problems. It has also been used in solving the Orr-Sommerfeld equation (Scott and Watts1977). On choosing either the symmetric or antisymmetric solutions, we can specify symmetry conditions at x = 0 and integrate the ordinary differential equations (28) to (30) from x = 0 to x = 1. As the continuity equation (30) is not suitable for the shooting method, we take the divergence of the Navier-Stokes equation, and for the pressure perturbation, we solve the resulting second-order equation:

$$\begin{array}{l}\left[i{S}^{-1}-\frac{4}{3}k(\bar{V}-c)\right]\left(\frac{{\mathrm{d}}^2}{\mathrm{d}{x}^2}-k^2\right)\widehat{p}\\ =\mathit{\text{ik}}{S}^{-1}\left[2\frac{\mathrm{d}\bar{V}}{\mathrm{d}x}{\widehat{v}}_x-k{(\bar{V}-c)}^2\widehat{p}\right]+\frac{4}{3}k\left(\frac{{\mathrm{d}}^2\bar{V}}{\mathrm{d}{x}^2}\widehat{p}+2\frac{\mathrm{d}\bar{V}}{\mathrm{d}x}\frac{\mathrm{d}\widehat{p}}{\mathrm{d}x}\right).\end{array}$$

(36)

The continuity equation (30) is treated as a boundary condition at x = 1. Defining a six-dimensional vector variable,$\mathit{q}(x){=}^{\mathrm{t}}\left({\widehat{v}}_x{\widehat{v}}_x^{\prime }{\widehat{v}}_z{\widehat{v}}_z^{\prime }\widehat{p}{\widehat{p}}^{\prime }\right)$, we can readily integrate the system of differential equations:

$$\frac{\mathrm{d}\mathit{q}}{\mathrm{d}x}=W(x)\mathit{q}(x),$$

(37)

where the 6 × 6 matrix W(x) can be derived from (28), (29), and (36). As the number of symmetry conditions at x = 0 is 3, we seek the remaining three independent solutions, q₁(x), q₂(x), and q₃(x), and we express the general solution as q = C₁q₁ + C₂q₂ + C₃q₃. The boundary conditions, (30) and (32) to (35), at x = 1 are subsequently reduced to a linear homogeneous system of algebraic equations for C₁, C₂, C₃, A, and B. The matching condition is that the determinant of the 5 × 5 coefficient matrix of the above linear equations becomes zero. In the case of the symmetric solution, the symmetry conditions at x = 0 should be${\widehat{v}}_x(0)={\widehat{v}}_z^{\prime }(0)={\widehat{p}}^{\prime }(0)=0$. Therefore, we simply choose three independent initial vectors as

$${\mathit{q}}_1(0)=\left(\begin{array}{c}0\\ 1\\ 0\\ 0\\ 0\\ 0\end{array}\right),{\mathit{q}}_2(0)=\left(\begin{array}{c}0\\ 0\\ 1\\ 0\\ 0\\ 0\end{array}\right),{\mathit{q}}_3(0)=\left(\begin{array}{c}0\\ 0\\ 0\\ 0\\ 1\\ 0\end{array}\right).$$

(38)

In the case of the antisymmetric solution, these conditions can be changed to match the opposite symmetry. We apply a standard Gram-Schmidt orthogonalization procedure to these three vectors to perform numerical integration regularly.

We use the fourth-order Runge-Kutta method for integration of (37). The computational grid we adopt is not an equidistant one, but it is finer near the fluid-solid boundary in case a viscous boundary layer develops. The number of grid points is 400. We confirm that the computer program produced a correct stability condition of a plane Poiseuille flow with rigid walls for which the boundary conditions at x = 1 were (30) and${\widehat{v}}_x(1)={\widehat{v}}_z(1)=0$. As the differential equations themselves are unchanged from the above well-known problem, the validity of our study lies in the implementation of boundary conditions (32) to (35). In the next section, we compare our results to those of Dunham and Ogden (2012) for a validity check, although their study adopted a long-wave approximation.

Results

Krauklis waves

As the solutions relevant to the excitation of volcanic tremors have a close relation with elastic surface waves confined in the fluid layer, here, we introduce some of the characteristics of the Krauklis (crack) wave before continuing on to the main part of our paper. Figure2 summarizes the dispersion relation of the Krauklis wave, which can be obtained as a special case of zero viscosity with no basic flow (for details, see Krauklis1962; Ferrazzini and Aki1987). When the wavelength is shorter than the fluid layer thickness, both the symmetric and antisymmetric fundamental modes approach the Stoneley wave (the surface wave trapped near a fluid-solid interface). There are also several higher-order modes if${\alpha }_0^{\ast }<{\beta }_s^{\ast }$, and the phase velocities of these higher-order solutions approach the fluid sound speed in the short wavelength limit. When the wavelength is longer, the antisymmetric solution approaches the nondispersive Rayleigh wave because the fluid layer thickness does not change and the fluid pressure has a negligible effect on the solid motion. In other words, the solid motion takes place as if the fluid-solid boundary is a free surface. The phase speed can be written as c = χ β (k ≪ 1), where χ denotes a constant slightly less than unity, satisfying the algebraic equation

$${\left(1-\frac{{\chi }^2}{2}\right)}^2-\sqrt{\left(1-{\chi }^2\right)\left(1-\frac{{\chi }^2{{\beta }_{\mathrm{s}}^{\ast }}^2}{{{\alpha }_{\mathrm{s}}^{\ast }}^2}\right)}=0.$$

(39)

Figure 2

Dispersion relation of Krauklis wave. The phase speed, c, of the Krauklis wave scaled by the S-wave velocity, β, of the surrounding solid is plotted as a function of wavelength, L = 2π/k. Shown are four cases wherein ϱ = 1,$\alpha /\beta =\sqrt{3}$ and (a) β = 0.5, (b) β = 1, (c) β = 1.5, and (d) β = 5. At a certain wavenumber, there is more than one solution. In this case, we show the 12 smallest phase velocities at most to simplify the plot. The solid and broken lines respectively represent the symmetric and antisymmetric solutions. The gray horizontal line represents the Rayleigh wave speed. The gray oblique line shows the asymptotic phase speed (40) of the symmetric solution in the long-wave limit.

For example, χ is about 0.9194 when${\alpha }_{\mathrm{s}}^{\ast }/{\beta }_{\text{s}}^{\ast }=\sqrt{3}$. In contrast, the phase velocity of the symmetric solution decreases to zero in the long-wave limit. The phase speed changes in inverse proportion to the square root of the wavelength and the asymptotic dispersion relation can be derived as follows:

$$c=\beta \sqrt{2k\varrho \left(1-\frac{{\beta }^2}{{\alpha }^2}\right)}(k\ll 1),$$

(40)

or in the dimensional form,

$$c^{\ast }=\sqrt{\frac{2\pi {\tilde{\mu }}_{\mathrm{s}}^{\ast }{a}^{\ast }}{{\rho }_0^{\ast }{L}^{\ast }}}(L^{\ast }\gg a^{\ast }),$$

(41)

where${\tilde{\mu }}_{\mathrm{s}}^{\ast }={\mu }_{\mathrm{s}}^{\ast }/(1-{\nu }_{\mathrm{s}})$ and ν_s denotes the solid Poisson ratio. This implies that the crucial parameter of the symmetric Krauklis wave is the speed calculated from the ratio of the solid shear modulus to the fluid density and that the elastic response is quasi-static where the solid density is irrelevant (see Dunham and Ogden2012).

In a narrow sense, the Krauklis wave refers to this symmetric, long-wave solution, (40) or (41), in which the crack slowly opens and closes symmetrically, but we use this term to mean both symmetric and antisymmetric long-wave solutions. Even when the fluid sound speed is slower than the elastic wave speeds of the solid, the higher-order modes are cut off in the long-wave limit and only one mode is permitted for each symmetry. It should be noted that the Krauklis wave dispersion relation is, in some sense, contrary to that of the Lamb wave propagating along an elastic plate (Lamb1917). In the case of the Lamb wave, the phase speed of the fundamental antisymmetric mode approaches zero when the wavelength becomes very long.

Figure3 shows an example of the wave motion of the Krauklis wave. When the wavelength is sufficiently long, the fluid motion associated with the symmetric wave is largely parallel to the z-axis. The small amplitude of v _x is a natural consequence of mass conservation (30) in a long-wave solution, and this is the reason why the solid particle velocity and hence, the wave speed, are small in the symmetric solution. The solid particle at the boundary moves along an elliptical path significantly elongated in the x-direction, rather than a circular path. In the antisymmetric wave, the fluid motion is largely parallel to the x-axis, following the wall motion. The solid particle motion is almost the same as that of the Rayleigh wave, which is more circular at the boundary than in the symmetric case (e.g., see Chapter 5 in Aki and Richards2002).

Figure 3

Fluid and solid motion of Krauklis wave along z -direction projected onto the xz -plane. In these examples, L = 25, β = 1,$\alpha =\sqrt{3}$, ϱ = 1, and the phase speeds are respectively 0.46 and 0.88 in (a) symmetric and (b) antisymmetric cases. A solid particle (open circle) sticking to the fluid-solid boundary (gray solid line) moves along an ellipse (closed dotted line). The black arrow represents the direction and magnitude of the particle motion. The white arrows represent the velocity field in the fluid layer. The arrow lengths are scaled by a common factor so that the speed can be compared between the symmetric and antisymmetric waves.

General properties of solutions

The solutions can be generally classified into two groups. The first group of solutions is intrinsic to the instability of a plane Poiseuille flow within rigid boundaries. In this group, the phase velocity, Re(c), is always positive and less than the Mach number. The wave propagating at a speed slightly slower than the peak speed of the main flow is known as the Tollmien-Schlichting wave, which is generally excited in the instability of a viscous shear layer (e.g., Schlichting and Gersten2000). The growth rate, Im(k c), largely depends on the Reynolds number, R = S^-1M. It is known that linear instability occurs at R ≃ 5800 and k ≃ 1 with antisymmetric disturbances. Figure4 illustrates an example of the distribution of phase velocities in the complex plane for cases in which the wavelength, L, is equal to 30,$\beta =\alpha /\sqrt{3}=1.5$, ϱ = 1, S = 10^-3, δ = 0.5, and M = 0, 0.4, and 0.8. The solutions belonging to the first category are seen in the range of 0 ≤ Re(c) ≤ M (the gray shaded part in Figure4). When M = 0.4, for example, there are three symmetric solutions and two antisymmetric solutions in the range shown, but they are all stable. In this parameter range, the growth rate decreases with increase in the Reynolds number.

Figure 4

Distribution of phase velocities in the complex plane. Red and blue lines respectively indicate the zero contour lines of the real and imaginary parts of the determinant of the 5 × 5 matrix representing the boundary conditions at x = 1. The solution, c, is represented by an intersection of red and blue lines and shown by a black dot. (a) The symmetric case where L = 20, β = 1.5,$\alpha /\beta =\sqrt{3}$, ϱ = 1, S = 10^-3, δ = 0.5, and M = 0, 0.4, and 0.8. The gray shading shows the range of 0 ≤ Re(c) ≤ M. (b) Same but for the antisymmetric case. Two solutions, denoted by A and B, originating from the antisymmetric Krauklis waves are observed around Re(c) = ±1.3. Close-ups of these solutions are shown on the right-hand side.

The second group of solutions emerges only when the surrounding region is elastic. As shown later, the flow instability in this second group can be interpreted as the instability of the Krauklis wave aided by the main laminar flow. Figure2 indicates that the speeds of the symmetric and antisymmetric Krauklis waves are respectively about c = 0.39β = 0.59 and c = 0.89β = 1.3, when β = 1.5 and L = 30. These surface waves appear in Figure4 as four solutions located at Re(c) ≃ ±0.59 and ±1.3. The growth rate is negative when M = 0 because of finite fluid viscosity. The decay in the symmetric wave is more rapid because the z-component of the flow perturbation is so great that a large amount of friction occurs at the fluid-solid boundary (see Figure3). As M increases, the symmetric solution traveling against the main flow (Re(c) ≃ -0.59) becomes more stabilized, and the wave speed slightly increases. In contrast, the growth rate of the symmetric solution traveling with the main flow (Re(c) ≃ 0.59) increases with increase in M, and it becomes positive at M = 0.8 (Figure4a). Instability also occurs in the antisymmetric solution, but the growing solution has a negative-phase speed. Figure4b shows that the growth rate of the antisymmetric solution at M = 0.4, traveling against the main flow at a speed of about -1.3, is still negative, but this solution becomes almost neutrally stable at M = 0.8. Therefore, it is evident that instability occurs with a smaller flow speed for the symmetric solution for this parameter set.

In general, as the flow speed increases, the growth rate, Im(k c), changes from negative to positive at a certain Mach number, which we call the critical Mach number, M_c. In the following subsections, we mainly discuss M_c and the corresponding phase speed of the neutrally stable solution as functions of the wavenumber and other dimensionless parameters. For computational reasons, the solutions are separately obtained for forward waves (Re(c) > 0) and backward waves (Re(c) < 0). We no longer discuss the solutions belonging to the first group because they require large Reynolds numbers to be destabilized and appear to be unimportant in explaining the excitation of volcanic tremors.

Symmetric solution

Figure5a,b,c,d illustrates the critical Mach number and the corresponding phase velocity of the symmetric solution, calculated for S = 10^-2, ϱ = 1, and$\alpha /\beta =\sqrt{3}$. When the wavelength is considerably more than unity, only one solution can be unstable. The unstable mode exhibits a wave traveling with the main flow. Both the Mach number and the phase velocity have a magnitude similar to the Krauklis wave speed (40), decreasing in inverse proportion to the square root of the wavelength and increasing in proportion to the shear modulus of the surrounding solid (see (41)). This suggests a close connection between the unstable mode and the slow Krauklis (crack) wave. As the wavelength increases, the critical Mach number increases with increase in the wavelength, and the unstable solution disappears in the very long wavelength range (Figure5a,b). When the wavelength is of the order of unity, M_c is generally large and more than one solution may become unstable. This is why the curve of M_c (and the phase speed) in Figure5a,b,c,d is discontinuous at some wavenumbers. At a particular parameter range, a solution propagating against the main flow is also destabilized. For example, there is a backward wave solution when β = 5 and 40 < L < 100, but it disappears when L > 100. This solution is probably related to the second mode of the symmetric Krauklis wave traveling at a speed comparable to the Rayleigh wave speed, as seen in the range of 10 < L < 100 in Figure2d. However, this backward-wave solution is of less importance because M_c is generally large. It is evident from Figure5a,b,c,d that a change in the functional form of the basic flow (the parameter δ) has no substantial impact on the symmetric solutions.

Figure 5

Plot of the critical Mach number and the corresponding wave phase speed. The Mach number, M _c , and the phase speed, |c|, are respectively represented by red and blue lines, and plotted as functions of the wavelength. Shown are four cases where S = 10^-2, ϱ = 1,$\alpha =\sqrt{3}\beta$, β = 1 and 5, and δ = 0 and 0.5. The wave solution propagating with (against) the main flow is represented by a solid (dotted) line. The orange curve indicates the phase speed of the fundamental mode of the crack wave (see Figure2). (a to d) Symmetric solutions. The gray line indicates the curve corresponding to the analytic expression (50) for M_c in the long-wave limit estimated from a slot-averaging approximation. The cross symbols indicate the results of Dunham and Ogden (2012), obtained from Figure 4 in their study. (e to h) Antisymmetric solutions. The gray lines indicate the curve corresponding to the analytic expressions (58), (63), and (64) for M_c in the long-wave limit. The Mach numbers (58) and (63) in a longer wavelength range are estimated from a slot-averaging approximation, but the steeper line is estimated empirically.

Figure6a,b shows the effect of fluid viscosity on the neutrally stable solution with a relatively long wavelength. The viscosity parameter, S, appears to have almost no relationship with the symmetric solutions when S is smaller than$O(\sqrt{k})$, particularly in the long-wave limit. This implies that the fluid viscosity is unimportant in the neutrally stable solution when S is small. As the critical Mach number and the phase velocity are both proportional to in the long-wave limit, it is natural to define a small parameter,$\epsilon =\sqrt{k}$, and express c and M_c as follows:

$$c=\epsilon c_0,M_{\mathrm{c}}={\xi }_{0}c=\epsilon {\xi }_0{c}_0,$$

(42)

Figure 6

Effect of fluid viscosity parameter on neutrally stable solution. The critical Mach number (red lines) and the phase speed (blue lines) are plotted as functions of the fluid viscosity parameter, S. Three cases are shown in each panel where ϱ = 1, β = 1,$\alpha =\sqrt{3}$, and L = 10² (solid lines); L = 10³ (broken line); and L = 10⁴ (dotted lines). (a, b) Symmetric solutions with δ = 0 and 0.5. (c, d) Antisymmetric solutions with δ = 0 and 0.5. The phase velocity of the antisymmetric solution does not vary significantly in this parameter range (|c| ≃ 0.9194).

where c₀ and ξ₀ denote O(1) constants to be determined. On eliminating A and B from (32) to (35) in which the viscous stress is neglected (i.e., S < ε) and taking the leading-order term with respect to ε, we obtain the following relations:

$$\begin{array}{l}{\epsilon }^2{\widehat{v}}_z(1)=-\frac{4(3-\delta )}{6-\delta }{\xi }_0{\widehat{v}}_x(1),\end{array}$$

(43)

$$\begin{array}{l}\epsilon \widehat{p}(1)=-\frac{2\varrho {\beta }^2}{c_0}\left(1-\frac{{\beta }^2}{{\alpha }^2}\right){\widehat{v}}_x(1).\end{array}$$

(44)

The above boundary conditions imply that${\widehat{v}}_z(x)\sim O(1)$,$\widehat{p}(x)\sim O(\epsilon )$, and${\widehat{v}}_x(x)\sim O({\epsilon }^2)$. Since the viscous force in the momentum equation is also negligible, the leading-order differential equations reduce to the following forms:

$$\frac{\mathrm{d}\widehat{p}}{\mathrm{d}x}=0,$$

(45)

$$\begin{array}{l}{\epsilon }^2\left[{\xi }_0\left(1-x^2\right)\left(1-\frac{\delta x^2}{6-\delta }\right)-1\right]{\widehat{v}}_z-\frac{4{\xi }_{0}x\left(3-\delta x^2\right)}{6-\delta }{\widehat{v}}_x+\frac{\epsilon \widehat{p}}{c_0}=0,\end{array}$$

(46)

$$\frac{\mathrm{d}{\widehat{v}}_x}{\mathrm{d}x}+{\epsilon }^2{\widehat{v}}_z=0.$$

(47)

The pressure is evidently a constant function. Although${\widehat{v}}_x(x)$ and${\widehat{v}}_z(x)$ can be obtained exactly by solving (46) and (47) simultaneously, we simply assume that${\widehat{v}}_z(x)$ is also a constant function; that is, we adopt a slot-averaging approximation and assume that$\widehat{p}(x)=\widehat{p}(1)$ and${\widehat{v}}_z(x)={\widehat{v}}_z(1)$. On integrating (47), we obtain the following:

$${\widehat{v}}_x(x)=-{\epsilon }^2{\widehat{v}}_z(1)x={\widehat{v}}_x(1)x.$$

(48)

The solutions derived from the above never satisfy the remaining condition (46), but we assume that the x-average of Equation 46, shown below, should be taken into account:

$$\begin{array}{l}\left[\frac{4(5-\delta )}{5(6-\delta )}{\xi }_0-1\right]{\epsilon }^2{\widehat{v}}_z(1)-\frac{4(5-\delta )}{5(6-\delta )}{\xi }_0{\widehat{v}}_x(1)+\frac{\epsilon \widehat{p}(1)}{c_0}=0.\end{array}$$

(49)

Finally, we find that four conditions, (43), (44), (48), and (49), for three unknowns, are satisfied only when

$$\left\{\begin{array}{l}{M}_{\mathrm{c}}=\frac{6-\delta }{4(3-\delta )}c,\\ c=\beta \sqrt{2k\varrho \left(1-\frac{{\beta }^2}{{\alpha }^2}\right)·\frac{5(3-\delta )}{5-3\delta }},\end{array}\right.$$

(50)

and k ≪ 1 and$S<O\left(\sqrt{k}\right)$. In the classical Poiseuille flow (δ = 0), the phase velocity is twice the peak speed of the main flow, and it is$\sqrt{3}\simeq 1.73$ times greater than the speed of the symmetric Krauklis wave, which can be inferred from (40) or (41). In other words, the wave associated with the long-wave and small-S symmetric instability can be interpreted as a destabilized Krauklis wave in a reference frame moving with a mean flow speed of$2(1-1/\sqrt{3})V_0^{\ast }\simeq 0.84V_0^{\ast }$. Figure5a,b,c,d indicates that the phase velocity in (50) is slightly overestimated, probably because${\widehat{v}}_x(x)$ and${\widehat{v}}_z(x)$ derived above are merely rough approximations. In contrast, the relation between M_c and c in (50) appears to be more valid. It is evident from (47) that compressibility is not important in the leading-order solution. Therefore, the fluid sound speed is not an important parameter in this regime. The phase speed in the long-wave symmetric instability can be scaled by the quantity,$\sqrt{{\tilde{\mu }}_{\mathrm{s}}^{\ast }/{\rho }_0^{\ast }}$, as seen in (41).

The above solution is essentially the same as that derived in the previous studies. Balmforth et al. (2005), assuming incompressibility for the working fluid, showed that the dimensional form of the critical flow speed is

$$V_0^{\ast }\sim {\beta }_{\mathrm{s}}^{\ast }\sqrt{\frac{\epsilon {\rho }_{\mathrm{s}}^{\ast }}{{\rho }_0^{\ast }}}=\sqrt{\frac{\epsilon {\mu }_{\mathrm{s}}^{\ast }}{{\rho }_0^{\ast }}},$$

(51)

where ε denotes the aspect ratio of the fluid-filled crack, which is comparable to our nondimensional wavenumber k. They explained that the relation (51) can be simply derived from a balance between the z-components of the inertial force of the order of${\rho }_0^{\ast }{V_0^{\ast }}^2/L^{\ast }$ and a pressure gradient of the order of p^∗/L^∗, wherein the pressure variation, p^∗, is equal to the elastic force of the order of${\mu }_{\mathrm{s}}^{\ast }{a}^{\ast }/L^{\ast }$. Dunham and Ogden (2012) pointed out that this long wavelength limit involves a quasi-static elastic response. They also showed that the solution that propagates with the main flow becomes unstable and that the solution that propagates against the main flow becomes more stabilized as the Mach number increases. Figure5a also shows the results of Dunham and Ogden (2012) taken from Figure4 in their study. We note that the phase speed is in good agreement with our result. However, the critical Mach number is slightly lower than our value, probably because an approximation was made in Dunham and Ogden (2012) when averaging the fluid momentum equation.

When S is greater than$O\left(\sqrt{k}\right)$, the long-wave solution no longer exists probably because in the stress boundary condition (26), the magnitudes of the fluid shear stress, τ _zx , and the last term,$-S{\bar{V}}^{\prime \prime }{u}_x$, are so large that no solution can satisfy this stress balance. The absence of the long-wave solution means that there is an optimal wavenumber, k ∼ O(S²), at which the most unstable mode grows. The existence of an overall minimum in the critical Mach number for the symmetric instability has not been reported in previous studies.

Antisymmetric solution

Our most remarkable finding is that the antisymmetric solution can be unstable, as illustrated in Figure5e,f,g,h for cases in which S = 10^-2, ϱ = 1, and$\alpha /\beta =\sqrt{3}$. When the wavelength is considerably more than the fluid layer thickness, there exists an unstable mode that travels against the main flow at a speed comparable to the antisymmetric Krauklis wave speed, which is almost the same as the nondispersive Rayleigh wave speed. When the basic laminar flow is driven by a uniform body force (δ = 0; Figure5e,g), this backward-wave solution yields a constant critical Mach number, M_c, which is almost equal to c/2, as has been seen in the symmetric instability. It is evident from comparing Figure5a,e or Figure5c,g that the overall minimum of M_c is found in the symmetric solution when δ = 0.

However, this situation becomes different when the driving force is deviated from the uniform profile even by a small degree. Figure5f,h indicates that in the presence of a driving force with δ = 0.5, the critical Mach number for the backward-wave instability decreases in inverse proportion to the wavelength, in contrast to the case of the symmetric solution wherein M_c decreases in inverse proportion to the square root of the wavelength. This difference is crucial in long-wave solutions. For example, the critical Mach number for the antisymmetric solution is nearly two orders of magnitude smaller than that for the symmetric solution when L = 10⁴. When the wavelength is shorter but still remains much greater than unity, M_c appears to be inversely proportional to L². The wave speed in this intermediate wavelength range is almost equal to the Rayleigh wave speed, and the propagation direction is also opposite to the main flow. When the wavelength is reasonably short (L ≲ O(10)), the backward-wave solution disappears in agreement with the reduction in the phase speed of the antisymmetric Krauklis wave (see Figure2). Instead, the solution propagating with the main flow becomes most unstable, though the critical Mach number is generally larger than unity. The emergence of the forward-wave solution in the smaller wavenumber range is also seen when δ = 0.

Figure6c,d shows the dependence of the long-wave antisymmetric solution on the viscosity parameter. The unstable mode with the constant critical Mach number, M_c = c/2, appears in the large- L and large-S regime in the case of δ = 0 (Figure6c). Our calculations show that M_c is almost constant when$S\gtrsim k\beta$, and that M _c increases with decrease in S when S ≲ k β. A similar tendency is observed in the solutions with δ ≠ 0 (Figure6d). The critical Mach number is independent of S when$S\gtrsim k\beta$ but increases with decreasing S when S ≲ k β. Figure6d shows that the critical Mach number is proportional to S^-1 in the small-S regime when δ ≠ 0. This confirms the fact that the antisymmetric instability with nonparabolic basic flow yields M_c values proportional to k in the long wavelength regime, M_c proportional to k² in the intermediate wavelength regime and that the regime boundary exists around S = k β. The phase speed is always coincident with the Rayleigh wave speed, c = χ β, irrespective of the magnitude of S and δ.

Figure7 illustrates the sensitivity of the long-wave solutions to the nonuniformity in the driving force. The critical Mach number generally decreases proportional to δ^-1. When L = 10³ and S = 10^-2, for example, the critical Mach number for the antisymmetric instability becomes significantly less than that for the symmetric instability at δ = 0.05, which means that the driving force at the fluid-solid boundary is 5% weaker than the force around the central part of the fluid layer. The effect of the reduction in the driving force near the boundary is more striking for longer wavelength solutions.

Figure 7

Effect of deviation from parabolic basic flow profile on neutrally stable solution. The critical Mach number is plotted as functions of the parameter, δ, calculated for ϱ = 1, β = 1,$\alpha =\sqrt{3}$; and (a) L = 10²; (b) L = 10³; and (c) L = 10⁴. Black (red) lines indicate the value at S = 10^-2 (S = 10^-4). Symmetric (antisymmetric) solutions are shown as solid (broken) lines.

For the condition that$S\gtrsim k\beta$, k ≪ 1, and δ ≠ 0, the critical Mach number is proportional to k. In this large-S (or very long wave) regime, the left-hand sides of (28) and (29) are essentially negligible. Therefore, both inertia and advection are not important for the neutrally stable solution. In the long-wave limit, the differential equations (29) and (30) yield the following leading-order solutions:

$${\widehat{v}}_x(x)={\widehat{v}}_x(1),{\widehat{v}}_z(x)={\widehat{v}}_z(1)x.$$

(52)

The pressure perturbation can be written as$\widehat{p}(x)=\widehat{p}(1)x$ in the slot-averaging approximation, but the differential equations (28) and (30) are satisfied only when${\widehat{v}}_x(1)\sim {\widehat{v}}_z(1)\sim O(1)$ and$\widehat{p}(1)$ is of the order of k. Under the stress boundary conditions, (34) and (35), the O(1) terms lead to the dispersion relation (39) of the Rayleigh wave because A ∼ B ∼ O(k^-2). The second largest O(S) terms in (34) and (35), which are greater than the O(k) terms, yield the relation

$${\widehat{v}}_z(1)=-{\bar{V}}^{\prime \prime }(1){û}_x(1),$$

(53)

where we define$u_x(x,z,t)={û}_x(x)exp\left[\mathit{\text{ik}}(z-\mathit{\text{ct}})\right]+\mathit{\text{c.c.}}$ and${û}_x(1)\sim O(k^{-1})$. The leading-order terms in the no-slip boundary conditions, (32) and (33), yield the following:

$$\begin{array}{l}{\widehat{v}}_x(1)=-\mathit{\text{kc}}{û}_x(1),\end{array}$$

(54)

$$\begin{array}{l}{\widehat{v}}_z(1)=\mathit{\text{kc}}{û}_z(1)-{\bar{V}}^{\prime }(1){û}_x(1),\end{array}$$

(55)

where we define$u_z(x,z,t)=i{û}_z(x)exp\left[\mathit{\text{ik}}(z-\mathit{\text{ct}})\right]+\mathit{\text{c.c.}}$ and c = χ β denotes the Rayleigh wave speed. As the Rayleigh wave solution implies the relation

$${û}_x(1):{û}_z(1)=1-\frac{1}{2}{\chi }^2:\sqrt{1-{\chi }^2},$$

(56)

we obtain from (53) to (55) that

$${\bar{V}}^{\prime \prime }(1)={\bar{V}}^{\prime }(1)-\mathit{\text{kc}}\frac{2\sqrt{1-{\chi }^2}}{2-{\chi }^2},$$

(57)

and

$$M_{\mathrm{c}}=-\frac{(6-\delta )\chi \sqrt{1-{\chi }^2}}{4\delta \left(2-{\chi }^2\right)}k\beta \left(k\ll 1,S\gtrsim k\beta ,\delta \ne 0\right).$$

(58)

The negative sign in the above expression indicates that the wave should propagate against the basic flow (i.e., c = χ β should be negative). In contrast to the case of the symmetric solution, the density contrast has no relation to the antisymmetric instability in the long wave and large-S limit, essentially because the pressure perturbation has almost no relation to the neutrally stable solution. In addition, the fluid sound velocity does not explicitly enter the expression for the critical flow speed,$V_0^{\ast }$. Therefore, it can be said that the stability condition for the large-S antisymmetric solution does not depend on the physical properties of the working fluid at all. As seen in Figure5f,h, our numerical calculations validate the asymptotic relation (58).

When the driving force is uniform throughout the fluid layer (δ = 0), the basic flow (6) exactly satisfies${\bar{V}}^{\prime \prime }(1)={\bar{V}}^{\prime }(1)$, and this peculiarity makes the condition (57) meaningless. In this case, the correct asymptotic behaviors are found to be${\widehat{v}}_x\sim O(1)$,${\widehat{v}}_z\sim {û}_x\sim {û}_z\sim O\left(k^{-1}\right)$, and$\widehat{p}\sim O(S)$. The left-hand sides of (28) and (29) are no longer negligible, and the leading-order solutions can be written as follows:

$$\left\{\begin{array}{l}{\widehat{v}}_x(x)={\widehat{v}}_x(1)\left[1-\frac{M}{c}\left(1-x^2\right)\right],\\ {\widehat{v}}_z(x)={\widehat{v}}_z(1)x=\left(-\frac{2M}{\mathit{\text{kc}}}{\widehat{v}}_x(1)+F\right)x,\\ \widehat{p}(x)=\widehat{p}(1)x=\frac{2\mathit{\text{iSM}}}{c}{\widehat{v}}_x(1)x.\end{array}\right.$$

(59)

The O(1) constant, F, included as a higher-order term in the expression for${\widehat{v}}_z$ is needed to correctly evaluate the boundary conditions. The leading-order terms in the stress boundary conditions (34) and (35) are of the order of S/k and require that

$${\widehat{v}}_z^{\prime }(1)\simeq -\frac{2M}{\mathit{\text{kc}}}{\widehat{v}}_x(1)=2M{û}_x(1),$$

(60)

but this condition is exactly the same as the velocity boundary condition for${\widehat{v}}_x$ (see (54)). The second largest O(1) terms lead to the dispersion relation of the Rayleigh wave. The third largest O(S) terms in (34) yield the same information as (60), but the boundary condition (35) leads to the relation

$$F=-2\mathit{\text{Mk}}{û}_z(1).$$

(61)

Finally, we need the no-slip condition (33) for${\widehat{v}}_z$. The leading-order O(k^-1) terms yield the same condition as (60), but the O(1) terms give the following:

$$F=\mathit{\text{kc}}{û}_z(1).$$

(62)

These two conditions (61) and (62) for the higher-order coefficient in${\widehat{v}}_z$ are compatible only when

$$M_{\mathrm{c}}=-\frac{c}{2}=-\frac{\chi \beta }{2}\left(k\ll 1,S\gtrsim k\beta ,\delta =0\right).$$

(63)

The negative sign above confirms that this is also a backward-wave solution.

When the viscosity parameter is considerably less than k β and the basic flow profile is deviated from the parabolic one (δ ≠ 0), it is suggested from Figure6d that M_c is proportional to k² and is inversely proportional to S. Our calculations for various dimensionless parameters indicate that

$$M_{\mathrm{c}}\simeq 0.3\times \frac{(6-\delta )\chi \sqrt{1-{\chi }^2}}{4\delta \left(2-{\chi }^2\right)}\frac{k^2{\beta }^2}{S}$$

(64)

when${\alpha }_{\mathrm{s}}^{\ast }/{\beta }_{\mathrm{s}}^{\ast }=\sqrt{3}$, k ≪ 1, and (k β)² ≲ S ≲ k β. However, despite our best efforts, we could not find an analytic form of small-S (or intermediate wavelength) antisymmetric solutions satisfying the governing equations (28) to (30) and the boundary conditions (32) to (35). In this case, the advection terms in (28) and (29) become more important than the viscous terms. We encountered the problem that the inviscid solution never satisfied the no-slip boundary condition. Therefore, we expect that a thin boundary layer needs to be formed in order for the solution to match the no-slip condition. Obtaining the analytic form of the small-S solution still remains our future work.

Physical mechanisms of instability

Because of finite viscous friction primarily arising at the fluid-solid boundary, the Krauklis wave is damped in time without the basic laminar flow. In the symmetric wave, the propagation speed is so slow that the wall motion is considerably slower than the perturbed fluid velocity (Figure3a). The viscous friction is primarily caused by the large amplitude of the z-component, v _z , of the perturbed fluid velocity, which is almost constant along the x-direction, and it is positive (negative) at the point at which the fluid layer is wider (narrower), provided that the wave travels along the z-direction. The friction in the symmetric Krauklis wave is therefore reduced by adding a basic flow in the same direction as the wave propagation (Figure8a). At the point at which the fluid layer is narrower, the basic flow velocity near the boundary is positive and can compensate the negative v _z . At the point at which the fluid layer is wider, the basic flow near the boundary is directed along the negative z-direction because the driving force is assumed to be unchanged in time and it needs to be balanced by the viscous force while maintaining the parabolic or quartic profile, (6) or (7). Therefore, the laminar basic flow and the perturbed fluid velocity are opposite near the boundary, and the friction is also reduced. As the basic flow speed increases further, the Krauklis wave is further amplified. The difference between parabolic and nonparabolic velocity profiles does not affect the mechanism of the symmetric instability because the important factor causing the instability is the magnitude of${\bar{V}}^{\prime }(\pm 1)$ and the magnitude of higher-order derivatives is irrelevant.

Figure 8

Schematic explanation for instability of a Krauklis wave traveling to the right. A solid particle (open circle) at the boundary moves along an ellipse (dotted line). White arrows represent the perturbed fluid and solid velocities associated with the Krauklis wave. The basic flow is represented by slender black arrows with its parabolic (or quartic) envelop. (a) The symmetric Krauklis wave inherently involves considerable friction at the boundary because of large amplitude of v _z , but it can be reduced if the basic flow is in the same direction as the wave propagation. (b) In the antisymmetric Krauklis wave, v _z is negligibly small but the solid particle motion at the boundary causes friction, which can be reduced if the basic flow is in the opposite direction to the wave propagation. This mechanism works only when the basic flow profile is nonparabolic.

It is readily confirmed that the above explanation leads to the relation between M_c and c in (50). As the solid particle at the boundary moves along an elliptical path significantly elongated in the x-direction (see Figure3a), the dimensional oscillation period T^∗ can be approximately expressed as

$$T^{\ast }=\frac{L^{\ast }}{c^{\ast }}\simeq \frac{4h^{\ast }}{\left|{\widehat{v}}_x^{\ast }\left(a^{\ast }\right)\right|},$$

(65)

where h^∗ ≪ a^∗ denotes the amplitude of the wall deformation. Mass conservation yields the following approximate equation:

$$\left|{\widehat{v}}_z^{\ast }\left(a^{\ast }\right)\right|·2a^{\ast }\simeq \left|{\widehat{v}}_x^{\ast }\left(a^{\ast }\right)\right|·\frac{L^{\ast }}{2},$$

(66)

because${\widehat{v}}_z^{\ast }\left(x^{\ast }\right)$ is almost constant along the x^∗-direction. The condition when the z-component of the flow velocity is zero at boundaries where the fluid layer is widest or narrowest is

$$\left|{\widehat{v}}_z^{\ast }\left(a^{\ast }\right)\right|=\left|{\bar{V}}^{\ast }\left(a^{\ast }\pm h^{\ast }\right)\right|\simeq \frac{2V_0^{\ast }{h}^{\ast }}{a^{\ast }},$$

(67)

where we assume that${\bar{V}}^{\ast }$ is parabolic in x^∗ for the sake of simplicity. From (65) to (67), it follows that$V_0^{\ast }\simeq c^{\ast }/2$. This is basically the same relation as that determined for the symmetric instability, irrespective of the viscosity parameter.

The physical mechanism of the antisymmetric instability with a nonparabolic basic flow profile (δ ≠ 0) is of particular interest because the critical Mach number can be significantly small and it can be more simply understood than the above symmetric case. The antisymmetric Krauklis wave comprises two parallel Rayleigh waves (Figure3b). The particle motion at the boundary is more circular than the symmetric wave, and the z-component of the perturbed fluid velocity is negligibly small. Let us suppose that the antisymmetric Krauklis wave travels along the z-direction. At the point at which the boundary wall comes toward the crack center, the particle motion is along the negative z-direction (opposite to the propagation direction of the Rayleigh wave). Therefore, the viscous friction can be reduced provided that the laminar basic flow is along the negative z-direction (same direction as the particle motion). At the point at which the boundary wall moves away from the crack center, the particle motion at the boundary is along the z-direction. If there is a basic flow in the direction opposite to the wave propagation, the flow direction near the boundary turns along the z-direction (same direction as the particle motion) so that the viscous friction can be reduced (Figure8b). In other words, the main flow can be sustained as if the moving wall is frictionless. If the flow speed is further increased, the flow can accelerate the wall motion and amplify the Rayleigh wave because of the no-slip boundary condition. In terms of energy, the work exerted on the fluid to maintain the plane laminar flow can be directly transferred to the energy of the Rayleigh wave via viscous drag. If the flow is in the same direction as the wave propagation, the wall motion decelerates because the friction increases. A remarkable point in this antisymmetric instability is that the flow and pressure perturbations merely play a passive role. Instead, the main flow itself drives the amplification of the Rayleigh wave.

Following the above argument, we can easily confirm the stability condition (58) observed in the very long wavelength range. If the dimensional amplitude of the wall deformation is represented by h^∗ and the angular velocity of the nearly circular particle motion is denoted by ω^∗, the speed of the particle motion is approximately h^∗ω^∗. The main flow speed at x^∗ = a^∗ ± h^∗ is of the order of$V_0^{\ast }·\left(h^{\ast }/a^{\ast }\right)$. Equating these two speeds yields

$$V_0^{\ast }\simeq a^{\ast }{\omega }^{\ast }=\frac{2\pi \chi a^{\ast }}{L^{\ast }}{\beta }_{\mathrm{s}}^{\ast },$$

(68)

because the phase speed of the Rayleigh wave is$c^{\ast }=L^{\ast }{\omega }^{\ast }/2\pi =\chi {\beta }_{\mathrm{s}}^{\ast }$. This result is almost identical to the stability condition (58) determined in the parameter space of S > k β and δ ≠ 0. Because ω^∗ = 2π/T^∗, the condition (68) indicates that in the large-S and long-wave antisymmetric instability with a nonparabolic velocity profile, the basic state is destabilized if the flow speed is so fast that the fluid particles move against the traveling Rayleigh wave through a distance longer than the fluid layer thickness in one oscillation period.

The above mechanism of antisymmetric instability with δ ≠ 0 is intrinsic to the characteristic of the Rayleigh wave particle motion at the boundary. In principle, the same mechanism could work in the symmetric solution too. The symmetric backward-wave solution, seen at β = 5 and propagating at a speed close to the Rayleigh wave speed (Figure5c,d), is probably caused by this mechanism, though the parameter range that allows it is quite limited.

Both mechanisms of antisymmetric and symmetric instabilities rely on the fact that friction of surface-wave propagation arising at the fluid-solid boundary can be reduced by adding a plane laminar flow and that a sufficiently fast fluid flow even amplifies the surface wave because of the no-slip boundary condition. This implies that although the fluid viscosity is insensitive to the stability conditions for the small-S symmetric solution and the large-S antisymmetric solution, the viscosity is an important factor regarding the growth of amplitude, as was also noticed by Dunham and Ogden (2012). In fact, in numerical calculations, we found that the magnitude of the growth rate, Im(k c), was sensitive to the viscosity parameter, S, and that it was generally large when S was large. This might also explain why the critical flow speed is generally a decreasing function of S, as can be observed in Figure6.

Similar mechanisms of instability are found in fluid dynamical problems such as a film flow down a slope (Kelly et al.1989) and Couette flows of two semi-infinite viscous fluids (Hooper and Boyd1983; Hinch1984). In such problems, a laminar shear flow is in contact with another medium and an unstable wave arises at the interface. In a film-flow problem, a thin fluid layer is in contact with air, and two fluids that have different viscosities are in contact in the latter example. It is also a common characteristic among these problems that the critical flow speed can be zero, at least in an ideal situation. When the interface slightly moves in the direction perpendicular to the laminar shear flow, discontinuity in the velocity field occurs at the interface, and a secondary flow is created to compensate this. In problems such as these, the secondary flow pattern propagates along the direction of the main flow due to advection. Instability in our symmetric solution can be analogously explained by this general idea, as has been pointed out by Balmforth and Mandre (2004). However, our antisymmetric solution exhibits a unique instability mechanism, in that the advection of the flow perturbation due to the main flow is not an important process, but the main flow itself drives acceleration of solid-particle motions. It would be more appropriate to state that our antisymmetric solution is a result of instability of the Rayleigh wave with the aid of the plane laminar flow, not a fluid dynamical instability aided by elastic wall deformation. The antisymmetric instability determined in this study may have been overlooked in the previous studies probably because the wave propagation was unexpectedly against the main flow, contrary to other known examples.

Rust et al. (2008) experimentally showed that a gas flow through a sheet-like elastically deformable conduit of a finite length, L^∗, is destabilized at a speed greater than$V_0^{\ast }\sim L^{\ast }{\omega }^{\ast }$. They hypothesized that the gas flow destabilizes the Rayleigh waves that exist at the fluid-solid boundaries as standing waves in the finite conduit, although the precise wave function is not fully described. We expect that the instability observed by Rust et al. (2008) corresponds to our solutions (63) obtained in the case of a purely parabolic basic flow (δ = 0), because we predict that$V_0^{\ast }=L^{\ast }{\omega }^{\ast }/4\pi$ in this solution. This correspondence may justify our mathematical and computational treatment, although there are still uncertainties as to whether the wave observed by Rust et al. (2008) travels against the main gas flow and how the wavelength is chosen in the finite conduit. As our study treats linear stability of an infinitely extended fluid layer, the relationship between the mechanisms of our antisymmetric instability and their studies’ instability of elastic normal modes remains to be discovered.

Discussion

Physical parameters relevant to volcanic tremors

In this section, we discuss a possible relation between our results and typical volcanic tremors. However, our intention is not to explain any specific tremor events because our model is too simplistic to be applied to a complicated volcanic setting. Our objective is to elucidate the fundamental physics of the plane-layer model and to present a general idea to explain some of the characteristics of volcanic tremors. We consider that a magma flow through a sheet-like dike induces an oscillatory instability. It is generally believed that a relatively less viscous basaltic magma, with a viscosity of 10¹ to 10⁴ Pa ·s, tends to intrude into bedrock from a reservoir, thus creating a planar crack. Geological observations indicate that the dike thickness, 2a^∗, varies from 0.1 m to more than 100 m, and the length can be up to 10 km or more (see Rubin1995and references therein), although the length is generally more difficult to measure than the thickness. The finiteness of the dike length means that our model may be applicable only when the wavelength, L^∗ = L a^∗, of the unstable solution is smaller than the dike length. Therefore, the upper limit of the nondimensional wavelength, L, is probably about O(10⁴) in the present study.

The magma flow should continue for a sufficiently long time period to explain the volcanic tremors. When hot magma intrudes into a cold crust, the flow ceases because of solidification of the magma, unless sufficient heat is supplied from the incoming flow. Theoretical consideration shows that the dike thickness needs to be greater than about 1 m in order for a continuous magma flow to be maintained (Bruce and Huppert1989). The threshold thickness depends on the temperature of the bedrock, the dike length, and the other physical properties of magma. In particular, the more viscous the intruding magma is, the wider the dike thickness should be (Petford et al.1993). An estimation of viscosities of various dike samples indicates that a dike of 1-m thickness is typically created by a magma with a viscosity of 10 to 100 Pa ·s and that thicker dikes tend to be formed from more viscous magmas (Wada1994). Therefore, in this section, we consider a typical dike thickness of about 1 m; however, we argue the effect of size on the wave generation of our model.

The physical properties of magma vary greatly with its chemical composition, temperature, and other factors such as the existence of bubbles and crystals. The magma density and the sound speed are respectively assumed to be${\rho }_0^{\ast }=2700$ kg/m³ and${\alpha }_0^{\ast }=10^3$ m/s, though the latter may vary by a factor of a few units or more if the magma contains bubbles. However, as shown later, we consider that the fluid sound speed is not of great importance to our discussion. The elastic wave speed of the surrounding bedrock also varies, basically increasing with the depth in the uppermost crust. We adopt${\beta }_{\mathrm{s}}^{\ast }\sim 2.5$ km/s and${\alpha }_{\mathrm{s}}^{\ast }=\sqrt{3}{\beta }_{\mathrm{s}}^{\ast }\sim 4.3$ km/s in this section.

The flow speed$V_0^{\ast }$ at the center of the dike is not an observationally well-constrained quantity, but it is probably not faster than several meters per second, and in most cases it is in the range of 0.1 to 1 m/s (Rubin1995; Best2003). In our model, the flow speed at the center of the dike is given by

$$V_0^{\ast }=\frac{K_0^{\ast }(6-\delta ){a^{\ast }}^2}{12{\zeta }^{\ast }},$$

(69)

where$K_0^{\ast }$ denotes the magnitude of K^∗ at the center of the dike (see (8)). If the driving force is simply magma buoyancy along the vertical direction and the fluid density is uniform,$K_0^{\ast }$ is equal to$\left({\rho }_{\mathrm{s}}^{\ast }-{\rho }_0^{\ast }\right)g^{\ast }$, where g^∗ denotes the acceleration due to gravity. Adopting parameters such as${\rho }_{\mathrm{s}}^{\ast }-{\rho }_0^{\ast }=300$ kg/m ³, ζ^∗ = 10³ Pa ·s, and 2a^∗ = 1 m, for example, we obtain$V_0^{\ast }=0.37$ m/s.

The observation quantity that we consider to be most relevant to volcanic tremors is the oscillation period T^∗. Observations indicate that T^∗ mostly ranges from 0.2 to 2 s (Chouet1996). One of the most unexplained features of volcanic tremors is that the oscillation period lies in such a narrow range at nearly every volcano. Another important index quantifiable in our model is the timescale of the onset. In general, volcanic tremors often show an emergent onset in the sense that the signal gradually emerges from the background noise, and determining the instant of the wave arrival is almost impossible (Konstantinou and Schlindwein2002). Based on previously published reports on volcanic tremors that occurred during eruptions, McNutt and Nishimura (2008) showed that about one half of the events exhibit an exponential increase in the oscillation amplitude at the beginning of eruption and that the e-folding time, τ^∗, ranges from several minutes to several hours. In our linear stability analysis, the e-folding time can be evaluated by τ^∗ = L^∗/(2π Im(c^∗)).

Applicability of symmetric solution

Under the constraints described in the last subsection, the symmetric instability cannot suitably explain the excitation of volcanic tremors. In this subsection, we discuss the case where δ = 0 because the critical Mach number is almost independent of δ in the symmetric solutions (see Figure7). When the wavelength, L^∗, is sufficiently greater than the fluid layer thickness, the critical flow speed can be expressed as

$$V_0^{\ast }=\frac{L^{\ast }}{2T^{\ast }}$$

(70)

because M_c is almost equal to c/2 when δ = 0, irrespective of the magnitude of the viscosity parameter. Therefore, the wavelength needs to be around 1 m to explain long-period tremors of T^∗ ∼ 1 s with a moderate flow speed of$V_0^{\ast }\sim 1$ m/s. As the nondimensional wavelength L must be sufficiently long to reduce the critical flow speed, we need a very thin fluid-filled crack. In our problem, the critical flow speed is generally a decreasing function of the wavelength so that instability will occur with the longest wavelength as permitted by the finite dike length. red Therefore, L^∗ cannot reasonably be around 1 m. If the symmetric solution explained volcanic tremors, the oscillation period would be considerably longer than 1 s, or the flow speed would be considerably faster than 1 m/s. This is the primary reason for the difficulty in using symmetric solutions, which has also been pointed out in previous studies (Balmforth et al.2005; Rust et al.2008; Dunham and Ogden2012).

When the viscosity parameter S is smaller than$O(\sqrt{k})$, the asymptotic result (50) gives

$$\begin{array}{l}{V}_0^{\ast }\simeq 210\times {\varrho }^{\frac{1}{3}}{\left(\frac{{\beta }_{\mathrm{s}}^{\ast }}{2.5\text{km/s}}\right)}^{\frac{2}{3}}{\left(\frac{2a^{\ast }}{1\text{m}}\right)}^{\frac{1}{3}}{\left(\frac{1\text{s}}{T^{\ast }}\right)}^{\frac{1}{3}}\left[\text{m/s}\right],\end{array}$$

(71)

where we assume that${\alpha }_{\mathrm{s}}^{\ast }=\sqrt{3}{\beta }_{\mathrm{s}}^{\ast }$. This implies that long-period volcanic tremors of T^∗ ∼ 1 s can be excited by the symmetric instability if a magma flows through a dike of moderate thickness at a speed of more than 100 m/s, or if the flow speed is slowed down to 1 m/s but the dike thickness is micron-scale. As$V_0^{\ast }$ only weakly depends on${\beta }_{\mathrm{s}}^{\ast }$, some reduction in the S wave speed of the surrounding bedrock does not alter the conclusion.

Applicability of antisymmetric solution

When the basic flow profile is parabolic in x, driven by a uniform forcing (δ = 0), the critical Mach number for the antisymmetric instability is generally larger than in the symmetric case. Therefore, it is more unlikely that the antisymmetric instability is an origin of volcanic tremors. However, we show in Section ‘Antisymmetric solution’ that a nonparabolic basic flow can be destabilized at an anomalously slow flow speed with perturbations of antisymmetric fluid and solid motion. In this section, we discuss the antisymmetric instability of such nonparabolic basic flows and we mainly use δ = 0.5, which means that the forcing${\bar{K}}^{\ast }\left(\pm a^{\ast }\right)$ at the fluid-solid boundary is${\bar{K}}_0^{\ast }/2$.

As the phase velocity associated with the long-wave antisymmetric instability is almost equal to the nondispersive Rayleigh wave speed, the oscillation period is simply given by

$$T^{\ast }=\frac{L^{\ast }}{\chi {\beta }_{\mathrm{s}}^{\ast }}.$$

(72)

Therefore, if χ = 0.9194 and${\beta }_{\mathrm{s}}^{\ast }=2.5$ km/s, a long-period wave of T^∗ = 1 s is created by the antisymmetric instability of a magma-flowing dike with a wavelength of 2.3 km, which would be acceptable because there are many geological examples in which the dike length reaches several kilometers. The long-wave antisymmetric instability changes its characteristic at S ∼ k β, or equivalently at${\zeta }^{\ast }\sim {\zeta }_0^{\ast }$. Equating (58) and (64) gives the threshold viscosity

$$\begin{array}{ll}{\zeta }_0^{\ast }& \simeq 1400\times \left(\frac{{\rho }_0^{\ast }}{2700{\text{kg/m}}^3}\right)\left(\frac{{\beta }_s^{\ast }}{2.5\text{km/s}}\right)\left(\frac{2.3\text{km}}{L^{\ast }}\right)\\ \times {\left(\frac{2a^{\ast }}{1\text{m}}\right)}^2\left[\text{Pa}·\text{s}\right],\end{array}$$

(73)

where we adopt χ = 0.9194 and δ = 0.5. From the analytic expressions (58) and (64), we can generally estimate the critical flow speed:

$$V_0^{\ast }=2.9\times \left\{\begin{array}{ll}\left(\frac{{\beta }_s^{\ast }}{2.5\text{km/s}}\right)\left(\frac{2.3\text{km}}{L^{\ast }}\right)\left(\frac{2a^{\ast }}{1\text{m}}\right)\left[\text{m/s}\right]& \left({\zeta }^{\ast }>{\zeta }_0^{\ast }\right),\\ \left(\frac{1400\text{Pa}·\text{s}}{{\zeta }^{\ast }}\right){\left(\frac{{\beta }_s^{\ast }}{2.5\text{km/s}}\right)}^2{\left(\frac{2.3\text{km}}{L^{\ast }}\right)}^2{\left(\frac{2a^{\ast }}{1\text{m}}\right)}^3\left[\text{m/s}\right]& \left({\zeta }^{\ast }<{\zeta }_0^{\ast }\right),\end{array}\right.$$

(74)

or, using (72), we can obtain the following:

$$V_0^{\ast }=2.9\times \left\{\begin{array}{ll}\left(\frac{2a^{\ast }}{1\text{m}}\right)\left(\frac{1\text{s}}{T^{\ast }}\right)\left[\text{m/s}\right]& \left({\zeta }^{\ast }>{\zeta }_0^{\ast }\right),\\ \left(\frac{1400\text{Pa}·\text{s}}{{\zeta }^{\ast }}\right){\left(\frac{2a^{\ast }}{1\text{m}}\right)}^3{\left(\frac{1\text{s}}{T^{\ast }}\right)}^2\left[\text{m/s}\right]& \left({\zeta }^{\ast }<{\zeta }_0^{\ast }\right).\end{array}\right.$$

(75)

The above estimation (illustrated in Figure9) supports the speculation that a long-period volcanic tremor is caused by the antisymmetric instability of a magma with viscosity of approximately 10³ Pa ·s, flowing at a speed of approximately 3 m/s, in a dike with a thickness of approximately 1 m and a length of several kilometers. If the magma has a smaller viscosity in the range of 10 to 100 Pa ·s, a large-scale dike with a length of 10 km or longer is needed to cause instability, provided that the dike thickness remains at 1 m and the flow speed is not greater than 3 m/s. In such a case, the resulting wave would have a very long period of approximately 5 s.

Figure 9

Dimensional estimate of critical flow speed and oscillation period for antisymmetric instability. Our analytical expressions of the critical flow speed ($V_0^{\ast }$ [m/s]; solid lines) and the resulting oscillation period (T^∗ [s]; broken line) for the antisymmetric instability are plotted as functions of the wavelength, L^∗ [m], wherein we assume that${\beta }_{\text{s}}^{\ast }={\alpha }_{\text{s}}^{\ast }/\sqrt{3}=2.5$ km/s,${\rho }_0^{\ast }=2700$ kg/m ³, and δ = 0.5. Four cases are shown wherein the thickness 2a^∗ of the fluid-filled crack is (a) 0.2 m, (b) 0.5 m, (c) 1 m, and (d) 5 m. In each case, the fluid viscosity is varied from 10 to 10⁵ P a ·s. There is a tendency for a more viscous fluid to be destabilized at a slower flow speed. The oscillation period does not depend on the viscosity and the crack thickness. The gray shaded part denotes the wavelength range that produces a period of 0.2 s < T^∗ < 2 s.

If the dike is thinner, instability can occur at a slower flow speed and a shorter wavelength, and hence, the oscillation period can be lesser (Figure9a,b). For example, a wavelength of approximately 500 m is possible in a dike with a thickness of 20 cm, if the magma viscosity is greater than 10³ Pa ·s and the flow speed at the center of the dike is 3 m/s. Even if the magma has a smaller viscosity of 10 to 100 Pa ·s, flow instability occurs at the same flow speed but the required wavelength is about 3 km. If the wavelength is 5 km and the viscosity is greater than 100 Pa ·s, a flow speed of 0.3 m/s is sufficient to cause instability (Figure9a). In contrast, if the dike is thicker, the antisymmetric instability generally requires a faster flow speed and a longer wavelength. Even if the wavelength can possibly range up to 5 km, the required flow speed is 10 m/s when the dike’s thickness is 5 m and the magma viscosity is greater than 10⁴ Pa ·s. The required flow speed unreasonably increases when the viscosity is smaller (Figure9d). It should be noticed that the above argument is based on the assumption of${\beta }_{\mathrm{s}}^{\ast }=2.5$ km/s and δ = 0.5. If the S-wave speed of the elastic surrounding was slower or the driving force near the dike wall was more significantly reduced, the resulting oscillation period would be longer and the critical flow speed would be slower.

If our model with a nonparabolic basic flow is applicable, the observational constraint of an upper bound of approximately 2 s in the oscillation period of volcanic tremors indicates that there should be an upper bound in the wavelength of the Rayleigh wave traveling along a dike (see (72)). Considering that the dike length should at least be longer than the wavelength, we would expect that the upper bound of the dike length would be approximately 10 km in most cases. The lower bound, approximately 0.2 s, of the tremor period would simply mean that there is an upper bound of a few meters per second in the magma flow speed and this might imply that there is a lower bound, approximately 0.2 m, in the dike thickness. These upper and lower bounds for dike length and thickness and magma flow speed are reasonably in agreement with observationally and theoretically expected values. We mentioned in Section ‘Physical parameters relevant to volcanic tremors’ that one of the most mysterious features of volcanic tremors is that the oscillation period falls in a narrow range. In this respect, an important factor in the narrowing of the period range would be that the critical flow speed in the antisymmetric solutions strongly depends on the wavelength (proportional to L^∗-2) when the wavelength is intermediate but weakly depends on the wavelength (proportional to L^∗-1) when the wavelength is very long.

A constraint on the tremor period can also be approached using a different point of view. Provided that a constant driving force$K_0^{\ast }$ at the center of the dike is given and the flow speed given by (69) meets the critical speed (64), we obtain the relationship:

$$T^{\ast }=\frac{L^{\ast }}{\chi {\beta }_{\mathrm{s}}^{\ast }}\simeq 3.6\times \sqrt{\frac{{\rho }_0^{\ast }{a}^{\ast }}{\delta K_0^{\ast }}}.$$

(76)

The fact that the oscillation period of the small-S (or intermediate wavelength) antisymmetric instability, which occurs in a wide parameter range as shown in Figure9, does not depend on the magma viscosity and the wavelength and only weakly depends on the dike thickness (and the density too) may explain why tremor events in different places and on different occasions tend to produce similar periods. It is evident from (73) to (75) that the fluid sound speed,${\alpha }_0^{\ast }$, and the density,${\rho }_{\mathrm{s}}^{\ast }$, of the surrounding bedrock are not related to the antisymmetric instability when the wavelength is sufficiently long. In particular, the magma sound speed can be significantly reduced with the existence of any small amount of gas phase (Kieffer1977). The absence of${\alpha }_0^{\ast }$ in the expressions for the critical flow speed and the oscillation period may also explain why long-period volcanic tremors of T^∗ ∼ 1 s occur at many different places.

Because the antisymmetric Krauklis wave is a surface wave concentrated around the fluid layer, its amplitude decays exponentially at distances far from the fluid layer. However, the decay is not very rapid because the wave relevant to volcanic tremors has a long wavelength of several kilometers. The solution (31) indicates that the e-folding distance is about 0.4 × L^∗ when χ = 0.9194. Therefore, the surface wave would be detectable unless the source is embedded deeper than several kilometers. Even if the source depth is greater, the wave might be detectable because of a body wave emission from the edge of the finite dike.

We have found that the dimensionless growth rate, Im(kc), is roughly proportional to S in our numerical calculations for the antisymmetric instability. As seen in Figure10, the e-folding time, τ^∗ = L^∗/(2π Im(c^∗)), calculated for various parameters obeys an approximate relation:

$${\tau }^{\ast }\sim \frac{{\rho }_0^{\ast }{a^{\ast }}^2{\beta }_{\mathrm{s}}^{\ast }}{{\zeta }^{\ast }{V}_0^{\ast }},$$

(77)

Figure 10

Time constant of exponential growth in wave amplitude. The e-folding time for the amplitude of linearly unstable antisymmetric solution to grow is plotted as a function of the peak flow speed at the center of the fluid layer. Calculation results are shown for three wavelengths of L^∗ = 3.2 km (solid line), 1.6 km (broken line), and 0.8 km (dotted line) and for four fluid viscosities of ζ^∗ = 10⁴ Pa ·s (black line), 10³ Pa ·s (red line), 10² Pa ·s (blue line), and 10 Pa ·s (orange line). In these examples, we set${\alpha }_0^{\ast }=1$ km/s,${\beta }_{\mathrm{s}}^{\ast }=2.5$ km/s,${\alpha }_{\mathrm{s}}^{\ast }=4.3$ km/s,${\rho }_{\mathrm{s}}^{\ast }={\rho }_0^{\ast }=2700$ kg/m ³, δ = 0.5, and the fluid layer thickness, 2a^∗, is (a) 0.2 m, (b) 0.5 m, and (c) 1 m.

though the wavelength is also an important parameter when the flow speed is close to the critical speed. The calculated e-folding time basically lies within the observationally acceptable range of several minutes to several hours (McNutt and Nishimura2008), supporting the applicability of our model. The basic flow speed is given as (69) in our model. Therefore, the e-folding time is essentially proportional to the inverse of$K_0^{\ast }$, thereby implying that the timescale of the onset of instability is largely controlled by the driving force of the magma flow. This is not surprising because the energy transfer rate from the basic laminar flow to the Rayleigh wave in the antisymmetric instability is identical to the rate of the work done by the viscous drag at the boundary wall. Some tremors show an abrupt onset. This could be explained by a wave emission due to dike propagation or wall fracture, or alternatively, this may suggest that our model cannot necessarily explain every volcanic tremor.

By studying both Figures9 and10, we conclude that the antisymmetric motion caused by flow-induced instability in a magma-filled dike reasonably explains some characteristics of volcanic tremors, such as the oscillation period and the timescale of the onset. The dike thickness should be 1 m at most and the length should be several kilometers. Magma viscosities of 10 to 10⁴ Pa ·s lie in the appropriate range, but a less viscous magma will flow in a thinner dike. Higher flow speed is desirable, but a speed of a few meters per second is sufficient to cause instability.

The most important assumption we made is that the basic flow profile is quartic in x, not parabolic as seen in the classical Poiseuille flow. This means that the driving force around the center of the dike is significantly greater than near the boundary. The nonparabolic flow profile may be caused by strongly temperature-dependent magma viscosity. Because the magma flow through a dike is generally cooled from the outside bedrock, the viscosity is higher near the boundary than in the central part, and the magma even solidifies in part near the bedrock. We assumed that the viscosity is uniformly constant but the driving force is nonuniform in order to simplify the governing equations and the boundary conditions. However, a greater forcing around the center of the dike filled with uniform-viscosity magma and a reduction in magma viscosity around the central part of the dike with a uniform forcing have a mathematically similar effect on the functional form of the basic flow profile. If the viscosities at the central part and near the boundary differ by an order of magnitude, our assumption of δ = 0.5 in this section would not be so unreasonable.

Volcanic tremors still have some important characteristics that cannot be explained by our model. Most importantly, we cannot explain the excitation of two or more simultaneous spectral peaks because we predict that the observed oscillations are primarily composed of the Rayleigh wave excited by the most unstable (i.e., the longest wavelength) mode. One possibility is that two or more magma flows occur at once, such as in an event that creates what is known geologically as a ‘dike swarm.’ Another possibility is that resonance occurs at a different location from the dike. In the antisymmetric instability, the Rayleigh wave propagates against the magma flow and hence toward a magma reservoir in the deeper part, which vibrates because of a continuous energy input through the Rayleigh wave, creating harmonic oscillations prescribed by the shape of the reservoir. The fluid flow in our model extends infinitely along the y- and z-directions. When the flow is confined in a finite dike, characteristic higher-order oscillations may be excited. In some tremors, the spectral peaks gradually shift to higher or lower frequencies in concert (Konstantinou and Schlindwein2002; McNutt2005). However, the above possibilities could still be incapable of explaining this gliding of the spectral peaks. It is possible that we now need to consider the effects of the finite dike dimensions and the finite amplitude disturbances, as well as the nonlinear time-dependent behaviors of the model, to explain every aspect of volcanic tremors.

Conclusions

We have discussed the possibility of flow-induced tremors in active volcanoes using a model mimicking a magma flow through a sheet-like dike. A plane Poiseuille flow is linearly destabilized with an extraordinarily slow flow speed when the surrounding infinite medium is elastically deformable, and the wavelength of the elastic deformation is considerably longer than the fluid layer thickness. We have found that the critical flow speed can be significantly reduced when the main flow has a nonparabolic velocity profile, driven by nonuniform forcing in the fluid. In this case, the most unstable mode exhibits an antisymmetric (flexural) deformation that propagates against the flow direction at the speed of an elastic Rayleigh wave. The critical flow speed decreases inversely proportional to the wavelength, and the asymptotic solution is derived when the wavelength is very long (or the viscosity is high). The physical mechanism of this instability can be understood as a simple friction drag on the Rayleigh wave particle motion caused by the main flow. We have shown that magma flowing at a speed of a few meters per second is sufficient to produce an amplification of the Rayleigh wave in an acceptable timescale, provided that the magma viscosity is 10 to 10⁴ Pa ·s and that the dike has a thickness of 0.2 to 1 m and a length of several kilometers. Our model predicts that the resulting oscillation period falls into a relatively narrow range because of natural constraints on the dike length and magma flow speed, and the predicted timescales substantially agree with the observations.

Problems to be examined in future works are summarized as follows. Our model has a serious problem of incapability in explaining oscillations containing several distinct spectral peaks. We recognize that this cannot be resolved within the framework of linear analysis of infinitely extending flows. It is required to discuss the effects of the finite dimensions of the fluid layer and also nonlinear behaviors of the model. Laboratory experiments may be useful in elucidating this complicated wave generation system. We have not obtained analytic solutions associated with the antisymmetric instability in an intermediate wavelength (or a low-viscosity) regime and with the symmetric instability in a high-viscosity regime. The physical mechanism of the former instability is particularly important because we predict that some volcanic tremors could be explained in that parameter regime. A discussion of the wave function that could not be obtained in our numerical calculations is at least needed. A comparison with observational data is not deemed to be sufficient. We also need to study the consistency of our model, making a comparison with the primary oscillation period, the e-folding time at the exponentially growing initial phase, magma viscosity, magma flow rate, and the seismic radiation properties in actual volcanoes. This paper has not dealt with the possible relation between tremors and long-period volcanic earthquakes that decay in a short timescale. To discover whether these two phenomena are created by an identical process is an important subject for future research. We hope that this paper triggers a fruitful discussion toward an understanding of the physical processes occurring under volcanoes.

Finally, we would like to mention the applicability of our model to other fields. It is well known that tremor-like seismic waves are observed, for example, in hydraulic fracturing of rock formations and in glacial icequakes, which may have a link to the propagation of fluid-filled cracks and fluid flows through a conduit embedded in an elastic medium (e.g., see Ferrick et al.1982 and references therein). As water has a smaller viscosity of approximately 10^-3 Pa ·s, the viscosity parameter S tends to be small, and hence the Reynolds number, R = S^-1M, would be sufficiently large to cause turbulence. For example, even a slow water flow of 0.1 m/s through a thin crack with a thickness of 1 cm has a large Reynolds number of approximately 10³. Even though the main fluid flow is already turbulent, the long wavelength (symmetric or antisymmetric) instability discussed in this paper could be applicable because viscous drag is its primary physical mechanism. An attractive possibility would be that a flow-induced instability explains long-period nonvolcanic tremors.