2.1 Linearized dynamics formulation

We consider a Boussinesq, 2-D flow slice model in the zonal (streamwise) vertical (cross-stream) plane ( $x{-}z$ ), with a zonally uniform basic state (denoted by overbars) which varies with height and is in hydrostatic balance. The momentum and continuity equations, linearized around this base state are:

(2.1) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{D}u}{\text{D}t}=-w\overline{U}_{z}-\frac{1}{\unicode[STIX]{x1D70C}_{0}}\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}x}, & \displaystyle\end{eqnarray}$$

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{D}w}{\text{D}t}=b-\frac{1}{\unicode[STIX]{x1D70C}_{0}}\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}z}, & \displaystyle\end{eqnarray}$$

(2.3) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{D}b}{\text{D}t}=-wN^{2}, & \displaystyle\end{eqnarray}$$

(2.4) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}x}+\frac{\unicode[STIX]{x2202}w}{\unicode[STIX]{x2202}z}=0. & \displaystyle\end{eqnarray}$$

Here $\text{D}/\text{D}t\equiv \unicode[STIX]{x2202}/\unicode[STIX]{x2202}t+\overline{U}\unicode[STIX]{x2202}/\unicode[STIX]{x2202}x$ is the linearized material derivative; $\boldsymbol{u}=(u,w)$ is the perturbation velocity vector; $\overline{U}(z)$ is the zonal mean flow; $p$ is the perturbation pressure; $\unicode[STIX]{x1D70C}$ is the perturbation density; $\unicode[STIX]{x1D70C}_{0}$ is a constant reference density and $b=-\unicode[STIX]{x1D70C}g/\unicode[STIX]{x1D70C}_{0}$ is the perturbation buoyancy. The squared buoyancy (or Brunt–Väisälä) frequency is defined as $N^{2}\equiv -(g/\unicode[STIX]{x1D70C}_{0})\,\text{d}\overline{\unicode[STIX]{x1D70C}}/\text{d}z=\overline{b}_{z}$ , where $\overline{\unicode[STIX]{x1D70C}}$ is the mean density, $g$ denotes gravity, $\overline{b}\equiv -\overline{\unicode[STIX]{x1D70C}}g/\unicode[STIX]{x1D70C}_{0}$ is the mean buoyancy and the subscript $z$ denotes the vertical derivative.

Equation set (2.1)–(2.4) can be transformed straightforwardly into a single equation in terms of the perturbation vorticity $q\equiv \unicode[STIX]{x2202}w/\unicode[STIX]{x2202}x-\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}z$ and vertical displacement $\unicode[STIX]{x1D701}$ fields:

(2.5) $$\begin{eqnarray}\displaystyle \frac{\text{D}}{\text{D}t}(q+\overline{q}_{z}\unicode[STIX]{x1D701})=-\overline{b}_{z}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D701}}{\unicode[STIX]{x2202}x}, & & \displaystyle\end{eqnarray}$$

where $\text{D}\unicode[STIX]{x1D701}/\text{D}t=w$ (from the kinematic condition) and $\overline{q}_{z}=-\overline{U}_{zz}$ is the mean vorticity gradient. For homogeneous fluids ( $\overline{b}_{z}=0$ ) the left-hand side indicates that vorticity perturbation is generated by vertical advection of the mean vorticity (the Rossby mechanism, sketched in figure 1 a,b), whereas in density stratified plane Couette flows (flows with constant shear, i.e.  $\overline{q}_{z}=0$ ), the right-hand side indicates that vorticity is generated due to the buoyancy restoring mechanism illustrated in figure 1(c,d).

2.2 Pseudo-momentum and pseudo-energy conservation

We assume zonal boundary conditions to be periodic and the vertical velocity to vanish at the upper and lower horizontal boundaries. Hence the domain integration of the cross-stream vorticity flux vanishes, i.e. $\langle wq\rangle =0$ . Hence, multiplying (2.5) by $\unicode[STIX]{x1D701}$ and integrating by parts yields the conservation of pseudo-momentum ${\mathcal{P}}$ :

(2.6) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}{\mathcal{P}}=0, & & \displaystyle\end{eqnarray}$$

where

(2.7) $$\begin{eqnarray}\displaystyle {\mathcal{P}}\equiv \left\langle \unicode[STIX]{x1D701}\left(q+\frac{\overline{q}_{z}}{2}\unicode[STIX]{x1D701}\right)\right\rangle . & & \displaystyle\end{eqnarray}$$

For the homogeneous case of zero stratification $q=-\overline{q}_{z}\unicode[STIX]{x1D701}$ , the familiar expression

(2.8) $$\begin{eqnarray}\displaystyle {\mathcal{P}}_{h}=\frac{1}{2}\langle \unicode[STIX]{x1D701}q\rangle =-\left\langle \frac{q^{2}}{2\overline{q}_{z}}\right\rangle =-\left\langle \frac{\overline{q}_{z}}{2}\unicode[STIX]{x1D701}^{2}\right\rangle , & & \displaystyle\end{eqnarray}$$

is recovered and leads to the Rayleigh inflection point condition. For the stratified Couette-like flow of constant shear we simply obtain

(2.9) $$\begin{eqnarray}\displaystyle {\mathcal{P}}_{c}=\langle \unicode[STIX]{x1D701}q\rangle . & & \displaystyle\end{eqnarray}$$

Thus, modal instability, for which ${\mathcal{P}}=0$ , can be satisfied by pairs of waves with opposite $(\unicode[STIX]{x1D701},q)$ sign relations that mutually amplify each other in accordance with figure 2(c,d).

Defining $P\equiv \unicode[STIX]{x1D701}(q+\overline{q}_{z}\unicode[STIX]{x1D701}/2)$ as the integrand of ${\mathcal{P}}$ , pseudo-energy conservation is obtained from (2.1)–(2.4) when noting that

(2.10) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}{\mathcal{E}}=-\langle \overline{U}wq\rangle =-\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}\langle \overline{U}P\rangle , & & \displaystyle\end{eqnarray}$$

where ${\mathcal{E}}\equiv \langle E\rangle$ is the sum of the eddy kinetic and available potential energies:

(2.11) $$\begin{eqnarray}\displaystyle \langle E\rangle \equiv \langle EKE\rangle +\langle APE\rangle =\left\langle \frac{1}{2}(u^{2}+w^{2})\right\rangle +\left\langle \frac{\overline{b}_{z}}{2}\unicode[STIX]{x1D701}^{2}\right\rangle . & & \displaystyle\end{eqnarray}$$

The pseudo-energy conservation then reads

(2.12) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}{\mathcal{H}}=0, & & \displaystyle\end{eqnarray}$$

where

(2.13a,b ) $$\begin{eqnarray}\displaystyle {\mathcal{H}}\equiv \langle H\rangle ;\quad H=E+\overline{U}P. & & \displaystyle\end{eqnarray}$$

For modal instability the pseudo-energy integral vanishes, therefore $\langle \overline{U}P\rangle =-\langle E\rangle <0$ . Hence, in both the homogeneous and the Couette-like cases, this implies that $\langle \overline{U}\unicode[STIX]{x1D701}q\rangle <0$ , which is the condition for counter-propagation (figure 2 c,d). For the homogeneous case we obtain the familiar Fjørtoft condition $\langle \overline{U}q^{2}/(2\overline{q}_{z})\rangle >0$ indicating positive correlation between the mean flow $\overline{U}$ and the mean vorticity gradient $\overline{q}_{z}$ . For the more general set-up of stratified shear flow we show in appendix A that the Howard–Miles criterion (Drazin & Reid Reference Drazin and Reid2004) for modal instability is related to the vanishing of PE. The simultaneous satisfaction of the two conditions of ${\mathcal{P}}=0$ and ${\mathcal{H}}=0$ are therefore in agreement with the conditions of mutual amplification and phase locking between pairs of vorticity waves (figure 2 d).

3.1 Plane waves in constant stratification and zero shear

For the sake of simplicity let us first consider the case of constant stratification and mean flow ( $N^{2}=\overline{b}_{z}=\text{const}_{1}$ , $\overline{U}=\text{const}_{2}$ ). Equation set (2.1)–(2.4) admits the familiar plane wave solution of the form of $\text{e}^{-\text{i}\unicode[STIX]{x1D703}}\text{e}^{\text{i}(kx+mz)}$ , where the phase: $\unicode[STIX]{x1D703}=\unicode[STIX]{x1D714}t+\unicode[STIX]{x1D703}_{0}$ ( $k$ and $m$ are the zonal and vertical wavenumber components and $k$ is assumed positive, $\unicode[STIX]{x1D714}$ is the wave frequency and $\unicode[STIX]{x1D703}_{0}$ is the initial phase). The dispersion relation is given by

(3.1a,b ) $$\begin{eqnarray}\displaystyle c=\frac{\unicode[STIX]{x1D714}}{k}=\overline{U}+{\hat{c}};\quad {\hat{c}}=\frac{\hat{\unicode[STIX]{x1D714}}}{k}=\pm \frac{N}{\sqrt{k^{2}+m^{2}}}, & & \displaystyle\end{eqnarray}$$

where $c$ is the phase speed in the zonal direction and $\hat{\unicode[STIX]{x1D714}}$ and ${\hat{c}}$ denote the intrinsic frequency and phase speed in the reference frame of the mean flow. It is straightforward to show that the zonal-averaged wave energy is equi-partitioned between its kinetic and potential counterparts so that $\overline{E}=\overline{b_{z}\unicode[STIX]{x1D701}^{2}}$ (Bühler Reference Bühler2009). Defining the zonal-averaged wave action (here after WA) as $\overline{A}\equiv \overline{E}/\hat{\unicode[STIX]{x1D714}}$ (note that action in classical mechanics is traditionally referred by the symbols $J$ or $I$ , however in fluid mechanics, wave action is usually referred by the symbol $A$ , see Bühler (Reference Bühler2009)) we obtain its relations with the zonal-averaged PM and PE:

(3.2a,b ) $$\begin{eqnarray}\displaystyle \overline{P}=k\overline{A}=\frac{\overline{E}}{{\hat{c}}};\quad \overline{H}=\overline{E+UP}=\unicode[STIX]{x1D714}\overline{A}=\frac{\unicode[STIX]{x1D714}}{\hat{\unicode[STIX]{x1D714}}}\overline{E}. & & \displaystyle\end{eqnarray}$$

We note then that the sign of the contribution of the wave to PM is equal to the sign of the intrinsic frequency, whereas the contribution to PE is positive unless ${\hat{c}}$ has the opposite sign of $\overline{U}$ and $|{\hat{c}}|<|\overline{U}|$ , that is when the wave propagates counter the mean flow, however with an intrinsic phase speed that is not large enough to overcome the mean flow advection.

From the second equation in (3.2) we can deduce that

(3.3) $$\begin{eqnarray}\displaystyle \overline{H}=\dot{\unicode[STIX]{x1D703}}\overline{A}\Rightarrow \frac{\unicode[STIX]{x2202}\overline{H}}{\unicode[STIX]{x2202}\overline{A}}=\dot{\unicode[STIX]{x1D703}}, & & \displaystyle\end{eqnarray}$$

and since $\overline{H}$ is time independent

(3.4) $$\begin{eqnarray}\displaystyle \dot{\overline{H}}=\frac{\unicode[STIX]{x2202}\overline{H}}{\unicode[STIX]{x2202}\overline{A}}\dot{\overline{A}}+\frac{\unicode[STIX]{x2202}\overline{H}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}\dot{\unicode[STIX]{x1D703}}=0\Rightarrow \frac{\unicode[STIX]{x2202}\overline{H}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}=-\dot{\overline{A}}. & & \displaystyle\end{eqnarray}$$

Equations (3.3) and (3.4) provide the canonical action ( $\overline{A}$ ) – angle ( $\unicode[STIX]{x1D703}$ ) representation of the linearized plane wave dynamics. (As mentioned previously, in classical mechanics the action is obtained from the circulation integral $J=\oint \boldsymbol{u}\,\text{d}\boldsymbol{x}$ . For stratified shear flow the circulation (of the mean flow plus perturbation) should be evaluated on an undulating plane of constant density where the baroclinic torque is zero and circulation is conserved. Under linearized wave dynamics Bühler (Reference Bühler2009) showed that the Eulerian representation of the zonal component of $J$ has a constant contribution from the initial steady mean flow and a contribution that is proportional to the negative of PM. The latter is the $O(\unicode[STIX]{x1D716}^{2})$ mean flow response. Hence $J$ is proportional to the negative of the wave action, so in principle we should define the wave action as $-\overline{A}$ and the angle as $-\unicode[STIX]{x1D703}$ . However, since in the context of fluid dynamics wave action is defined as $\overline{A}$ , we follow this convention in the text.) Obviously in this simple case $\overline{H}$ is not a function of the wave phase (angle) and indeed $\overline{A}$ is time independent. This is the standard A-A formulation in which action is a constant of motion. This simple example has been brought in order to pave the road for the generalized A-A description of the interaction between interfacial vorticity waves in stratified shear flows, where the total action is conserved but the action of each wave is generally time dependent.

3.2 Single interfacial vorticity wave

The presence of shear distorts the structure of plane waves. Thus, let us now consider interfacial waves that are resilient to shear and preserve their untilted structures. For a stratified shear flow, where both the mean vorticity and the stratification are discontinuous at some level $z=z_{0}$ :

(3.5a,b ) $$\begin{eqnarray}\displaystyle \overline{q}_{z}=\unicode[STIX]{x0394}\overline{q}_{0}\unicode[STIX]{x1D6FF}(z-z_{0}),\quad \overline{b}_{z}=\unicode[STIX]{x0394}\overline{b}_{0}\unicode[STIX]{x1D6FF}(z-z_{0}), & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x0394}\overline{q}_{0}\equiv \overline{q}(z>z_{0})-\overline{q}(z<z_{0})$ , and $\unicode[STIX]{x0394}\overline{b}_{0}\equiv \overline{b}(z>z_{0})-\overline{b}(z<z_{0})$ , (2.5) then implies that the vorticity perturbation should as well have the form of a $\unicode[STIX]{x1D6FF}$ -function at $z=z_{0}$ . Thus, for a normal mode solution we may write

(3.6) $$\begin{eqnarray}\displaystyle q=\tilde{q}_{0}\text{e}^{\text{i}k(x-ct)}\unicode[STIX]{x1D6FF}(z-z_{0}). & & \displaystyle\end{eqnarray}$$

We can find the wave dispersion relation and structure on combining (2.5) with the kinematic condition $\text{D}\unicode[STIX]{x1D701}/\text{D}t=w$ at $z=z_{0}$ and express $w$ in terms of $q$ via the Green’s function formulation (e.g. Harnik et al. Reference Harnik, Heifetz, Umurhan and Lott2008):

(3.7) $$\begin{eqnarray}\displaystyle w(z)=\int _{z^{\prime }}q(z^{\prime })G(z^{\prime },z,k)\,\text{d}z^{\prime }, & & \displaystyle\end{eqnarray}$$

where $-k^{2}G+\unicode[STIX]{x2202}^{2}G/\unicode[STIX]{x2202}z^{2}=\text{i}k\unicode[STIX]{x1D6FF}(z-z^{\prime })$ . The Green’s function $G(z^{\prime },z,k)$ depends on the boundary conditions and on the zonal wavenumber $k$ . For open flows, which we assume here for simplicity

(3.8) $$\begin{eqnarray}\displaystyle G(z^{\prime },z,k)=-\frac{\text{i}}{2}\text{e}^{-k|z-z^{\prime }|}, & & \displaystyle\end{eqnarray}$$

so that $w(z=z_{0})=-\text{i}(\tilde{q}/2)\text{e}^{\text{i}k(x-ct)}$ . Different Green’s functions for different boundary conditions can be found in Heifetz & Methven (Reference Heifetz and Methven2005). The dispersion relation obtained satisfies

(3.9a,b ) $$\begin{eqnarray}\displaystyle c^{\pm }=\overline{U}_{0}+{\hat{c}}^{\pm };\quad {\hat{c}}^{\pm }=-\frac{\unicode[STIX]{x0394}\overline{q}_{0}}{4k}\pm \sqrt{\left(\frac{\unicode[STIX]{x0394}\overline{q}_{0}}{4k}\right)^{2}+\frac{\unicode[STIX]{x0394}\overline{b}_{0}}{2k}}. & & \displaystyle\end{eqnarray}$$

Thus for stable stratification ( $\unicode[STIX]{x0394}\overline{b}_{0}>0$ ), ${\hat{c}}^{+}$ is always positive and ${\hat{c}}^{-}$ is always negative. The terms associated with $\unicode[STIX]{x0394}\overline{q}_{0}$ capture the Rossby wave propagation mechanism (figure 1 a,b) and the term with $\unicode[STIX]{x0394}\overline{b}_{0}$ results from the buoyancy restoring force (figure 1 c,d). When $\unicode[STIX]{x0394}\overline{q}_{0}=0$ we recover the familiar deep water internal wave dispersion relation, whereas when $\unicode[STIX]{x0394}\overline{b}_{0}=0$ we recover the interfacial Rossby wave together with a degenerated solution of zero vorticity perturbation. The asymmetry between ${\hat{c}}^{+}$ and ${\hat{c}}^{-}$ results from the Rossby wave mechanism propagating the wave to the left of the mean vorticity gradient, whereas the buoyancy restoring force is even for both rightward and leftward propagation. The two interfacial waves at the interface satisfy

(3.10a,b ) $$\begin{eqnarray}q=[\tilde{q}_{0}^{+}\text{e}^{-\text{i}kc^{+}t}+\tilde{q}_{0}^{-}\text{e}^{-\text{i}kc^{-}t}]\unicode[STIX]{x1D6FF}(z-z_{0})\text{e}^{\text{i}kx},\quad \unicode[STIX]{x1D701}(z=z_{0})=[\unicode[STIX]{x1D701}_{0}^{+}\text{e}^{-\text{i}kc^{+}t}+\unicode[STIX]{x1D701}_{0}^{-}\text{e}^{-\text{i}kc^{-}t}]\text{e}^{\text{i}kx},\end{eqnarray}$$

where

(3.11) $$\begin{eqnarray}\displaystyle \tilde{q}_{0}^{\pm }=2k{\hat{c}}^{\pm }\unicode[STIX]{x1D701}_{0}^{\pm }. & & \displaystyle\end{eqnarray}$$

These are in agreement with figure 1, indicating that the wave whose vorticity and displacement structures are in (anti) phase propagates to the (left) right with respect to the local mean flow.

Despite the presence of a mean shear (generally $\overline{q}=-\overline{U}_{z}(z)\neq 0$ ) these interfacial waves preserve their untilted structure throughout the domain. Defining the perturbation streamfunction $\unicode[STIX]{x1D713}$ , so that $u=-\unicode[STIX]{x1D713}_{z}$ , $w=\unicode[STIX]{x1D713}_{x}$ , $q=\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D713}$ , and recall that $\langle EKE\rangle =-\langle \unicode[STIX]{x1D713}q\rangle /2$ , the domain integrated wave energy for each wave can be evaluated solely from the interface using (3.5), (3.6), (3.8) and (3.11):

(3.12) $$\begin{eqnarray}\displaystyle {\mathcal{E}}^{\pm }=[\langle EKE\rangle +\langle APE\rangle ]^{\pm }=\left[-\left\langle \frac{\unicode[STIX]{x1D713}q}{2}\right\rangle +\left\langle \frac{\overline{b}_{z}}{2}\unicode[STIX]{x1D701}^{2}\right\rangle \right]^{\pm }=\frac{k}{2}\left[({\hat{c}}^{\pm })^{2}+\frac{\unicode[STIX]{x0394}\overline{b}_{0}}{2k}\right]|{\unicode[STIX]{x1D701}_{0}}^{\pm }|^{2}, & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

where the shear prevents the equi-partition between the kinetic and the potential energies of the waves. Substitution of (3.11) in (3.12) indicates that the two interfacial waves are orthogonal to each other under the energy norm, i.e. $\langle E\rangle =\langle E\rangle ^{+}+\langle E\rangle ^{-}$ . This is due to the fact that the integrated energy results solely from the waves’ signatures at the interface. Generally it is PE (which is the sum of the wave energy and the second-order mean flow response), rather than the wave energy itself, that is conserved under linearized dynamics. The contribution of each wave to PM is obtained by substituting the displacement and vorticity of the waves at the interface (i.e.  $\tilde{q}_{0}^{\pm }$ and $\unicode[STIX]{x1D701}_{0}^{\pm }$ , and using (3.11)) in (2.7):

(3.13a,b ) $$\begin{eqnarray}\displaystyle {\mathcal{P}}^{\pm }=\pm \frac{1}{2}|{\unicode[STIX]{x1D701}_{0}}^{\pm }|^{2}\sqrt{\left(\frac{\unicode[STIX]{x0394}\overline{q}_{0}}{2}\right)^{2}+2k\unicode[STIX]{x0394}\overline{b}_{0}};\quad {\mathcal{P}}={\mathcal{P}}^{+}+{\mathcal{P}}^{-}. & & \displaystyle\end{eqnarray}$$

Thus, the contribution of the rightward (leftward) propagating interfacial wave to PE is positive (negative). Furthermore, as expected, neutral normal modes with different phase speeds are orthogonal with respect to a norm that is a conserved quantity (Held Reference Held1985). Defining the domain integrated WA as ${\mathcal{A}}\equiv {\mathcal{E}}/\hat{\unicode[STIX]{x1D714}}$ , it is straightforward to verify for the interfacial waves that

(3.14a,b ) $$\begin{eqnarray}\displaystyle {\mathcal{A}}^{\pm }=\frac{{\mathcal{P}}^{\pm }}{k};\quad {\mathcal{A}}={\mathcal{A}}^{+}+{\mathcal{A}}^{-}. & & \displaystyle\end{eqnarray}$$

Hence the interfacial waves are orthogonal as well with respect to WA. To complete the analogy with plane waves in the absence of shear, we note as well that

(3.15) $$\begin{eqnarray}\displaystyle {\mathcal{H}}={\mathcal{H}}^{+}+{\mathcal{H}}^{-}=(\unicode[STIX]{x1D714}{\mathcal{A}})^{+}+(\unicode[STIX]{x1D714}{\mathcal{A}})^{-}. & & \displaystyle\end{eqnarray}$$

Defining the phase angle for the interfacial waves $\unicode[STIX]{x1D703}^{\pm }$ to satisfy $\dot{\unicode[STIX]{x1D703}}^{\pm }=\unicode[STIX]{x1D714}^{\pm }$ , we obtain the canonical A-A formulation

(3.16) $$\begin{eqnarray}\displaystyle {\mathcal{H}}=(\dot{\unicode[STIX]{x1D703}}{\mathcal{A}})^{+}+(\dot{\unicode[STIX]{x1D703}}{\mathcal{A}})^{-}, & & \displaystyle\end{eqnarray}$$

where

(3.17a,b ) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}{\mathcal{H}}}{\unicode[STIX]{x2202}{\mathcal{A}}^{\pm }}=\dot{\unicode[STIX]{x1D703}}^{\pm },\quad \frac{\unicode[STIX]{x2202}{\mathcal{H}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}^{\pm }}=-\dot{{\mathcal{A}}}^{\pm }=0. & & \displaystyle\end{eqnarray}$$

3.3 Multiple interfaces

For the interfacial wave dynamics we discretize the continuous mean flow into piecewise linear profiles of vorticity and density (buoyancy):

(3.18a,b ) $$\begin{eqnarray}\displaystyle \overline{q}_{z}=\mathop{\sum }_{n=1}^{N}\unicode[STIX]{x0394}\overline{q}_{n}\unicode[STIX]{x1D6FF}(z-z_{n}),\quad \overline{b}_{z}=\mathop{\sum }_{n=1}^{N}\unicode[STIX]{x0394}\overline{b}_{n}\unicode[STIX]{x1D6FF}(z-z_{n}), & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x0394}\overline{q}_{n}=\overline{q}(z_{n+1}>z>z_{n})-\overline{q}(z_{n}>z>z_{n-1})$ , and $\unicode[STIX]{x0394}\overline{b}_{n}=\overline{b}(z_{n+1}>z>z_{n})-\overline{b}(z_{n}>z>z_{n-1})$ . This formulation may include interfaces with only density jumps, only vorticity jumps or both. Thus, it may be applied to different basic set-ups such as Rayleigh, Holmboe and Taylor–Caulfield profiles. In the next section we show that the linearized interfacial wave dynamics can be presented in the compact A-A form

(3.19a-d ) $$\begin{eqnarray}{\mathcal{H}}=\mathop{\sum }_{n=1}^{N}[(\dot{\unicode[STIX]{x1D703}}{\mathcal{A}})^{+}+(\dot{\unicode[STIX]{x1D703}}{\mathcal{A}})^{-}]_{n};\quad {\mathcal{A}}=\mathop{\sum }_{n=1}^{N}[{\mathcal{A}}^{+}+{\mathcal{A}}^{-}]_{n};\quad \frac{\unicode[STIX]{x2202}{\mathcal{H}}}{\unicode[STIX]{x2202}{\mathcal{A}}^{\pm }}_{n}=\dot{\unicode[STIX]{x1D703}}_{n}^{\pm };\quad \frac{\unicode[STIX]{x2202}{\mathcal{H}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}^{\pm }}_{n}=-\dot{{\mathcal{A}}}_{n}^{\pm },\end{eqnarray}$$

where ${\mathcal{H}}$ and ${\mathcal{A}}=k{\mathcal{P}}$ are the two constant of motion of the linearized monochromatic wave interaction dynamics. As opposed to (3.17) and to textbook examples in classical mechanics, the WA components, ${\mathcal{A}}_{n}^{\pm }$ , are generally non-zero due to interaction at a distance between remote interfacial waves. While it may be argued that the last two equations in (3.19) straightforwardly indicate canonical Hamilton equations, we have however chosen to refer to (3.19) as the ‘generalized A-A equations’. This is due to the additional constraint in our system – the sum of all the individual wave actions, ${\mathcal{A}}$ is constant. For modal instability of the form of $\text{e}^{\text{i}[kx-(\unicode[STIX]{x1D714}_{r}+\text{i}\unicode[STIX]{x1D714}_{i})t )]}$ , ${\mathcal{H}}=\unicode[STIX]{x1D714}_{r}{\mathcal{A}}=0$ , however

(3.20a,b ) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}{\mathcal{H}}}{\unicode[STIX]{x2202}{\mathcal{A}}^{\pm }}_{n}=\dot{\unicode[STIX]{x1D703}}_{n}^{\pm }=\unicode[STIX]{x1D714}_{r};\quad -\frac{1}{{\mathcal{A}}_{n}^{\pm }}\frac{\unicode[STIX]{x2202}{\mathcal{H}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}_{n}^{\pm }}=\left(\frac{\dot{{\mathcal{A}}}_{n}^{\pm }}{{\mathcal{A}}_{n}^{\pm }}\right)=2\unicode[STIX]{x1D714}_{i}. & & \displaystyle\end{eqnarray}$$

4.1 PM, energy, WA and PE for two interfaces

Consider now two interfaces located at $z_{1}$ and $z_{2}=(z_{1}+\unicode[STIX]{x0394}z)$ , so that $N=2$ in (3.18). Let us decompose the displacement and the vorticity perturbations at those interfaces into the interfacial waves discussed in § 3.2:

(4.1) $$\begin{eqnarray}\displaystyle & \displaystyle q_{1,2}=[\tilde{q}^{+}(t)+\tilde{q}^{-}(t)]_{1,2}\unicode[STIX]{x1D6FF}(z-z_{1,2})\text{e}^{\text{i}kx}=[Q^{+}(t)\text{e}^{-\text{i}\unicode[STIX]{x1D703}^{+}(t)}+Q^{-}(t)\text{e}^{-\text{i}\unicode[STIX]{x1D703}^{-}(t)}]_{1,2}\unicode[STIX]{x1D6FF}(z-z_{1,2})\text{e}^{\text{i}kx}, & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

(4.2) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D701}(z=z_{1,2})=[\tilde{\unicode[STIX]{x1D701}}^{+}(t)+\tilde{\unicode[STIX]{x1D701}}^{-}(t)]_{1,2}\text{e}^{\text{i}kx}=[Z^{+}(t)\text{e}^{-\text{i}\unicode[STIX]{x1D703}^{+}(t)}+Z^{-}(t)\text{e}^{-\text{i}\unicode[STIX]{x1D703}^{-}(t)}]_{1,2}\text{e}^{\text{i}kx}, & \displaystyle\end{eqnarray}$$

where

(4.3) $$\begin{eqnarray}\displaystyle \tilde{q}_{1,2}^{\pm }=[(2k{\hat{c}}^{\pm })\tilde{\unicode[STIX]{x1D701}}^{\pm }]_{1,2}, & & \displaystyle\end{eqnarray}$$

and

(4.4) $$\begin{eqnarray}\displaystyle {\hat{c}}_{1,2}^{\pm }=\left[-\frac{\unicode[STIX]{x0394}\overline{q}}{4k}\pm \sqrt{\left(\frac{\unicode[STIX]{x0394}\overline{q}}{4k}\right)^{2}+\frac{\unicode[STIX]{x0394}\overline{b}}{2k}}\right]_{1,2}. & & \displaystyle\end{eqnarray}$$

Here ${\hat{c}}_{1,2}^{\pm }$ is the intrinsic phase speed of the four waves (two at each interface) in the absence of interaction. To avoid confusion we emphasize that $[\dot{\unicode[STIX]{x1D703}}^{\pm }/k-\overline{U}]_{1,2}\neq {\hat{c}}_{1,2}^{\pm }$ , since $\unicode[STIX]{x1D703}_{1,2}^{\pm }$ is affected by the interaction at a distance with the waves at the opposed interface. All fields of each wave propagate in concert with the same instantaneous frequency $\dot{\unicode[STIX]{x1D703}}_{1,2}^{\pm }$ . Similarly the vorticity and the displacement wave amplitudes change in time due to this interaction, however since in this partition the waves’ structures are preserved, the ratio between vorticity and displacement, $(2k{\hat{c}}^{\pm })_{1,2}$ , always remains constant (either with or without interaction with the opposed interfacial waves). We therefore refer to these waves as the ‘building blocks’ of the linearized interfacial dynamics. While the structure of each wave is untilted (for instance, their far field cross-stream velocity $w$ remains untilted, as can be verified from (3.7), (3.8)) their superposition may yield a complex tilted structure.

We immediately note that PM preserves the same simple structure of (3.13):

(4.5a,b ) $$\begin{eqnarray}\displaystyle {\mathcal{P}}_{1,2}^{\pm }=\left[\pm \frac{1}{2}(Z^{\pm })^{2}\sqrt{\left(\frac{\unicode[STIX]{x0394}\overline{q}}{2}\right)^{2}+2k\unicode[STIX]{x0394}\overline{b}}\right]_{1,2};\quad {\mathcal{P}}=\mathop{\sum }_{n=1}^{2}({\mathcal{P}}^{+}+{\mathcal{P}}^{-})_{n}. & & \displaystyle\end{eqnarray}$$

The energy (2.11) however includes mixed terms between the interfaces since the streamfunction at each interface is contributed to both from the waves in situ and from the waves of the opposed interface. Using the Green function (3.8) we find that

(4.6) $$\begin{eqnarray}\displaystyle {\mathcal{E}} & = & \displaystyle \frac{k}{2}\mathop{\sum }_{n=1}^{2}\left[\left(({\hat{c}}^{+})^{2}+\frac{\unicode[STIX]{x0394}\overline{b}}{2k}\right)(Z^{+})^{2}+\left(({\hat{c}}^{-})^{2}+\frac{\unicode[STIX]{x0394}\overline{b}}{2k}\right)(Z^{-})^{2}\right]_{n}\nonumber\\ \displaystyle & & \displaystyle +\,k\text{e}^{-k|\unicode[STIX]{x0394}z|}\left[{\hat{c}}_{1}^{+}{\hat{c}}_{2}^{+}Z_{1}^{+}Z_{2}^{+}\cos (\unicode[STIX]{x1D703}_{2}^{+}-\unicode[STIX]{x1D703}_{1}^{+})+{\hat{c}}_{1}^{-}{\hat{c}}_{2}^{-}Z_{1}^{-}Z_{2}^{-}\cos (\unicode[STIX]{x1D703}_{2}^{-}-\unicode[STIX]{x1D703}_{1}^{-})\right.\nonumber\\ \displaystyle & & \displaystyle +\,\left.{\hat{c}}_{1}^{+}{\hat{c}}_{2}^{-}Z_{1}^{+}Z_{2}^{-}\cos (\unicode[STIX]{x1D703}_{2}^{-}-\unicode[STIX]{x1D703}_{1}^{+})+{\hat{c}}_{1}^{-}{\hat{c}}_{2}^{+}Z_{1}^{-}Z_{2}^{+}\cos (\unicode[STIX]{x1D703}_{2}^{+}-\unicode[STIX]{x1D703}_{1}^{-})\right].\end{eqnarray}$$

We write the terms in (4.6) symbolically as

(4.7) $$\begin{eqnarray}\displaystyle {\mathcal{E}}=\mathop{\sum }_{n=1}^{2}({\mathcal{E}}^{+}+{\mathcal{E}}^{-})_{n}+\mathop{\sum }_{i=1}^{2}\mathop{\sum }_{j=1}^{2}({\mathcal{E}}^{++}+{\mathcal{E}}^{--}+{\mathcal{E}}^{+-}+{\mathcal{E}}^{-+})_{i,j(i\neq j)}, & & \displaystyle\end{eqnarray}$$

where the first sum includes the waves’ self-contributions to the energy and the second double sum includes the mixed contributions. We define the interfacial WA (only for the self-contribution energy terms) as

(4.8) $$\begin{eqnarray}\displaystyle {\mathcal{A}}_{1,2}^{\pm }\equiv \left(\frac{{\mathcal{E}}}{\hat{\unicode[STIX]{x1D714}}}\right)_{1,2}^{\pm }=\left[\left({\hat{c}}^{\pm }+\frac{\unicode[STIX]{x0394}\overline{b}}{2k{\hat{c}}^{\pm }}\right)\frac{(Z^{\pm })^{2}}{2}\right]_{1,2}=\frac{{\mathcal{P}}_{1,2}^{\pm }}{k}, & & \displaystyle\end{eqnarray}$$

so that (3.14) holds. Then we can write PE as

(4.9) $$\begin{eqnarray}\displaystyle {\mathcal{H}}=k\mathop{\sum }_{n=1}^{2}[(\overline{U}+{\hat{c}}^{+}){\mathcal{A}}^{+}+(\overline{U}+{\hat{c}}^{-}){\mathcal{A}}^{-}]_{n}+\mathop{\sum }_{i=1}^{2}\mathop{\sum }_{j=1}^{2}({\mathcal{E}}^{++}+{\mathcal{E}}^{--}+{\mathcal{E}}^{+-}+{\mathcal{E}}^{-+})_{i,j(i\neq j)}. & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

4.2 Wave interaction equations

The following derivation is based on Harnik et al. (Reference Harnik, Heifetz, Umurhan and Lott2008), hereafter referred to as H08. For clarity we re-derive the essence of it with the notation of the current paper.

We wish to describe the wave interaction dynamics solely in terms of the waves’ displacements at the interfaces. Toward this end we first take the vorticity equation (2.5) and use (4.2), (4.3) to write

(4.10) $$\begin{eqnarray}\displaystyle \left[{\hat{c}}^{+}\frac{\text{D}\unicode[STIX]{x1D701}^{+}}{\text{D}t}+{\hat{c}}^{-}\frac{\text{D}\unicode[STIX]{x1D701}^{-}}{\text{D}t}=-\frac{1}{2k}(\unicode[STIX]{x0394}\overline{q}w+\text{i}k\unicode[STIX]{x0394}\overline{b}\unicode[STIX]{x1D701})\right]_{1,2}. & & \displaystyle\end{eqnarray}$$

After this, we implement the kinematic condition at the interfaces:

(4.11) $$\begin{eqnarray}\displaystyle \left[\frac{\text{D}}{\text{D}t}(\unicode[STIX]{x1D701}^{+}+\unicode[STIX]{x1D701}^{-})=w\right]_{1,2}, & & \displaystyle\end{eqnarray}$$

and express the vertical velocity using (3.7), (3.8), (4.3):

(4.12) $$\begin{eqnarray}w_{1,2}=-\frac{\text{i}}{2}(\tilde{q}_{1,2}+\tilde{q}_{2,1}\text{e}^{-k|\unicode[STIX]{x0394}z|})=-\text{i}k\left[({\hat{c}}^{+}\unicode[STIX]{x1D701}^{+}+{\hat{c}}^{-}\unicode[STIX]{x1D701}^{-})_{1,2}+({\hat{c}}^{+}\unicode[STIX]{x1D701}^{+}+{\hat{c}}^{-}\unicode[STIX]{x1D701}^{-})_{2,1}\text{e}^{-k|\unicode[STIX]{x0394}z|}\right].\end{eqnarray}$$

Then we substitute (4.12) in (4.10), (4.11), and after some algebra obtain the equations for the time variation of displacement of each interfacial wave:

(4.13) $$\begin{eqnarray}\displaystyle \dot{\unicode[STIX]{x1D701}}_{1,2}^{\pm }=-\text{i}k\left[(\overline{U}+{\hat{c}})_{1,2}\unicode[STIX]{x1D701}_{1,2}^{\pm }\,\pm \text{e}^{-k|\unicode[STIX]{x0394}z|}\left(\frac{{\hat{c}}^{\pm }}{{\hat{c}}^{+}-{\hat{c}}^{-}}\right)_{1,2}({\hat{c}}^{+}\unicode[STIX]{x1D701}^{+}+{\hat{c}}^{-}\unicode[STIX]{x1D701}^{-})_{2,1}\right], & & \displaystyle\end{eqnarray}$$

where the first term on the right-hand side is due to the self-interfacial wave dynamics and the second results from the interaction at a distance with the two waves at the opposed interface. Taking the imaginary part of (4.13) we get

(4.14) $$\begin{eqnarray}\displaystyle \dot{\unicode[STIX]{x1D703}}_{1,2}^{\pm } & = & \displaystyle k\left(\overline{U}+{\hat{c}}\right)_{1,2}^{\pm }\pm \frac{k\text{e}^{-k|\unicode[STIX]{x0394}z|}}{Z_{1,2}^{\pm }}\left(\frac{{\hat{c}}^{\pm }}{{\hat{c}}^{+}-{\hat{c}}^{-}}\right)_{1,2}\nonumber\\ \displaystyle & & \displaystyle \times \,[({\hat{c}}^{+}Z^{+})_{2,1}\cos (\unicode[STIX]{x1D703}_{2,1}^{+}-\unicode[STIX]{x1D703}_{1,2}^{\pm })+({\hat{c}}^{-}Z^{-})_{2,1}\cos (\unicode[STIX]{x1D703}_{2,1}^{-}-\unicode[STIX]{x1D703}_{1,2}^{\pm })].\end{eqnarray}$$

The first term on the right-hand side contains the advection of the mean flow at the interface and the self-propagation mechanism in the absence of interaction. The last two terms represent the interaction with the two waves at the opposed interface. The dependence of the interaction on the cosine of their phase difference indicates the mechanism of interaction. When the waves are in (out of) phase (the cosine is $1$ ( $-1$ )), the self and the induced cross-stream velocity are in (out of) phase, hence the waves help (hinder) each other to propagate (for more details the reader is kindly referred to the review paper (Carpenter et al. Reference Carpenter, Tedford, Heifetz and Lawrence2013)).

Next we multiply each of the four equations, represented by (4.14), with their counterparts at (4.8) and compare their product with (4.9). After some algebra we obtain the remarkable result:

(4.15) $$\begin{eqnarray}\displaystyle \mathop{\sum }_{n=1}^{2}[(\dot{\unicode[STIX]{x1D703}}{\mathcal{A}})^{+}+(\dot{\unicode[STIX]{x1D703}}{\mathcal{A}})^{-}]_{n}={\mathcal{H}}, & & \displaystyle\end{eqnarray}$$

from which

(4.16) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}{\mathcal{H}}}{\unicode[STIX]{x2202}{\mathcal{A}}^{\pm }}_{1,2}=\dot{\unicode[STIX]{x1D703}}_{1,2}^{\pm }. & & \displaystyle\end{eqnarray}$$

Finally, taking the real part of (4.13), multiplying it by $[({\hat{c}}^{\pm }+(\unicode[STIX]{x0394}\overline{b}/2k{\hat{c}}^{\pm }))Z^{\pm }]_{1,2}$ , and noting that $[({\hat{c}}^{\pm }+(\unicode[STIX]{x0394}\overline{b}/2k{\hat{c}}^{\pm }))=\pm ({\hat{c}}^{+}-{\hat{c}}^{-})]_{1,2}$ , we obtain

(4.17) $$\begin{eqnarray}\dot{{\mathcal{A}}}_{1,2}^{\pm }=-k\text{e}^{-k|\unicode[STIX]{x0394}z|}{\hat{c}}_{1,2}^{\pm }[({\hat{c}}^{+}Z^{+})_{2,1}\sin (\unicode[STIX]{x1D703}_{2,1}^{+}-\unicode[STIX]{x1D703}_{1,2}^{\pm })+({\hat{c}}^{-}Z^{-})_{2,1}\sin (\unicode[STIX]{x1D703}_{2,1}^{-}-\unicode[STIX]{x1D703}_{1,2}^{\pm })]Z_{1,2}^{\pm }.\end{eqnarray}$$

The above equation indicates that the growth of WA of each interfacial wave is solely due to the interaction with the opposed interfacial waves. The dependence of the interaction on the sine of their phase difference results from the mechanism of growth – the induced cross-stream velocity amplify the wave displacement of the opposed wave and this amplification is maximized when the waves are in quadrature (as in figure 2 d). For more details the reader is referred again to Carpenter et al. (Reference Carpenter, Tedford, Heifetz and Lawrence2013).

If we now differentiate (4.9) with respect to $\unicode[STIX]{x1D703}_{1,2}^{\pm }$ and compare with (4.17), after some algebra we indeed find that

(4.18) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}{\mathcal{H}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}^{\pm }}_{1,2}=-\dot{{\mathcal{A}}}_{1,2}^{\pm }. & & \displaystyle\end{eqnarray}$$

This completes the explicit derivation of the generalized A-A formulation for two interfaces. The procedure can be carried on systematically (not shown here) for multiple interfaces to obtain (3.19) for any $N$ integer number of interfaces.

4.3 The subset of counter-propagating waves

As indicated from the heuristic arguments in the introduction, we expect that the instability mechanism will be obtained mainly through action at a distance between the counter-propagating interfacial vorticity waves. Indeed, Rabinovich et al. (Reference Rabinovich, Umurhan, Harnik, Lott and Heifetz2011) analysed the Taylor–Caulfield instability (Caulfield Reference Caulfield1994) and showed that for a large range of Richardson numbers the growth rates, corresponding to the most unstable modes, are practically unaffected when the pro-propagating waves are neglected. Similar results have been obtained by Heifetz & Mak (Reference Heifetz and Mak2015) for the more complex non-Boussinesq dynamics of stratified shear flows (Heifetz & Mak Reference Heifetz and Mak2015), for swirling flow instability in rotating cylinders (Yellin-Bergovoy, Heifetz & Umurhan Reference Yellin-Bergovoy, Heifetz and Umurhan2017), for gravity–capillary waves (Biancofiore et al. Reference Biancofiore, Heifetz, Hoepffner and Gallaire2017) and even for Alfvén waves in magneto-hydrodynamic shear flows (Heifetz et al. Reference Heifetz, Mak, Nycander and Umurhan2015). Therefore, next we consider the subset dynamics of the counter propagating interfacial vorticity waves.

First we need to identify the counter-propagating waves. We note that the conservation of pseudo-momentum and pseudo-energy are unaffected by Galilean transformation. Hence in the frame of reference of the mean zonal mean velocity $\overline{U}_{m}\equiv (\overline{U}_{1}+\overline{U}_{2})/2$ the new Hamiltonian ${\mathcal{H}}_{m}\equiv \langle E+(\overline{U}-\overline{U}_{m})P\rangle$ is also a constant of motion (since both $\langle P\rangle$ and $\overline{U}_{m}$ are constant). Define the zonal mean flow at the reference frame of $\overline{U}_{m}$ as $\overline{U}_{1,2}^{\ast }\equiv (\overline{U}_{1,2}-\overline{U}_{m})=(\overline{U}_{2,1}-\overline{U}_{1,2})/2$ , we define the counter-propagating waves as the ones whose intrinsic phase speed has the opposite sign of $\overline{U}^{\ast }$ at their interface. For instance, let us assume that the shear at the two interfaces is as in figure 2 (level 1 is below level 2) so that the counter-propagating waves will be the ones of figure 2(d), i.e. $(\unicode[STIX]{x1D701}_{1}^{+},\unicode[STIX]{x1D701}_{2}^{-})$ .

Equation (4.17) indicates immediately that even if we begin with a pair of counter-propagating waves, the pro-propagating ones $(\unicode[STIX]{x1D701}_{1}^{-},\unicode[STIX]{x1D701}_{2}^{+})$ will be generated immediately due to the interaction. Hence the counter-propagating wave subset is just an approximation to the dynamics which is valid only when the waves’ amplitudes satisfy $Z_{pro}\ll Z_{counter}$ . Under this approximation the energy and PE become

(4.19) $$\begin{eqnarray}\displaystyle {\mathcal{E}} & = & \displaystyle \frac{k}{2}\left[\left(({\hat{c}}^{+})^{2}+\frac{\unicode[STIX]{x0394}\overline{b}}{2k}\right)(Z^{+})^{2}\right]_{1}+\frac{k}{2}\left[\left(({\hat{c}}^{-})^{2}+\frac{\unicode[STIX]{x0394}\overline{b}}{2k}\right)(Z^{-})^{2}\right]_{2}\nonumber\\ \displaystyle & & \displaystyle +\,k\text{e}^{-k|\unicode[STIX]{x0394}z|}{\hat{c}}_{1}^{+}{\hat{c}}_{2}^{-}Z_{1}^{+}Z_{2}^{-}\cos (\unicode[STIX]{x1D703}_{2}^{-}-\unicode[STIX]{x1D703}_{1}^{+}),\end{eqnarray}$$

(4.20) $$\begin{eqnarray}\displaystyle {\mathcal{H}} & = & \displaystyle k[(\overline{U}+{\hat{c}}^{+}){\mathcal{A}}^{+}]_{1}+k[(\overline{U}+{\hat{c}}^{-}){\mathcal{A}}^{-}]_{2}\nonumber\\ \displaystyle & & \displaystyle +\,k\text{e}^{-k|\unicode[STIX]{x0394}z|}{\hat{c}}_{1}^{+}{\hat{c}}_{2}^{-}Z_{1}^{+}Z_{2}^{-}\cos (\unicode[STIX]{x1D703}_{2}^{-}-\unicode[STIX]{x1D703}_{1}^{+}).\end{eqnarray}$$

The counter-propagating wave interaction equations become

(4.21) $$\begin{eqnarray}\displaystyle & \displaystyle \dot{\unicode[STIX]{x1D703}}_{1}^{+}=k(\overline{U}+{\hat{c}})_{1}^{+}+\frac{k}{2}\text{e}^{-k|\unicode[STIX]{x0394}z|}{\hat{c}}_{1}^{+}{\hat{c}}_{2}^{-}\frac{Z_{1}^{+}Z_{2}^{-}}{{\mathcal{A}}_{1}^{+}}\cos (\unicode[STIX]{x1D703}_{2}^{-}-\unicode[STIX]{x1D703}_{1}^{+}), & \displaystyle\end{eqnarray}$$

(4.22) $$\begin{eqnarray}\displaystyle & \displaystyle \dot{\unicode[STIX]{x1D703}}_{2}^{-}=k\left(\overline{U}+{\hat{c}}\right)_{2}^{-}+\frac{k}{2}\text{e}^{-k|\unicode[STIX]{x0394}z|}{\hat{c}}_{1}^{+}{\hat{c}}_{2}^{-}\frac{Z_{1}^{+}Z_{2}^{-}}{{\mathcal{A}}_{2}^{-}}\cos (\unicode[STIX]{x1D703}_{2}^{-}-\unicode[STIX]{x1D703}_{1}^{+})\,, & \displaystyle\end{eqnarray}$$

(4.23) $$\begin{eqnarray}\displaystyle & \displaystyle \dot{{\mathcal{A}}}_{1}^{+}=-k\text{e}^{-k|\unicode[STIX]{x0394}z|}{\hat{c}}_{1}^{+}{\hat{c}}_{2}^{-}Z_{1}^{+}Z_{2}^{-}\sin (\unicode[STIX]{x1D703}_{2}^{-}-\unicode[STIX]{x1D703}_{1}^{+})=-\dot{{\mathcal{A}}}_{2}^{-}, & \displaystyle\end{eqnarray}$$

from which it is clear that

(4.24) $$\begin{eqnarray}\displaystyle (\dot{\unicode[STIX]{x1D703}}{\mathcal{A}})_{1}^{+}+(\dot{\unicode[STIX]{x1D703}}{\mathcal{A}})_{2}^{-}={\mathcal{H}}, & & \displaystyle\end{eqnarray}$$

(4.25a,b ) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}{\mathcal{H}}}{\unicode[STIX]{x2202}{\mathcal{A}}^{+}}_{1}=\dot{\unicode[STIX]{x1D703}}_{1}^{+},\quad \frac{\unicode[STIX]{x2202}{\mathcal{H}}}{\unicode[STIX]{x2202}{\mathcal{A}}^{-}}_{2}=\dot{\unicode[STIX]{x1D703}}_{2}^{-} & & \displaystyle\end{eqnarray}$$

(4.26a,b ) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}{\mathcal{H}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}^{+}}_{1}=-\dot{{\mathcal{A}}}_{1}^{+},\quad \frac{\unicode[STIX]{x2202}{\mathcal{H}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}^{-}}_{2}=-\dot{{\mathcal{A}}}_{2}^{-}. & & \displaystyle\end{eqnarray}$$

The same procedure can be carried on naturally for any number ( ${\geqslant}2$ ) of interfaces, where the counter-propagating wave at interface $n$ is the one whose intrinsic phase speed ${\hat{c}}_{n}$ is of opposite sign to $\overline{U}_{n}^{\ast }=(\overline{U}-\overline{U}_{m})_{n}$ , $\overline{U}_{m}$ being the mean zonal mean velocity at all of the interfaces. This counter-propagation interaction approximation reduces the complexity of the dynamics by a factor of two. As shown by Rabinovich et al. (Reference Rabinovich, Umurhan, Harnik, Lott and Heifetz2011) and Carpenter, Guha & Heifetz (Reference Carpenter, Guha and Heifetz2017), a dense enough grid of interfaces can accurately resolve complex shear flow dynamics, including the rapid changes across critical layers.