1 Introduction
Internal gravity waves move vertically through a continuously stratified fluid, transporting momentum and irreversibly accelerating the background flow where they break. Even before breaking, however, localized internal wavepackets induce transient flows that migrate with the group velocity of the wavepacket. In part, this is a consequence of the divergence of the momentum flux, which is zero far from the wavepacket and largest near its centre. The result is a ‘divergent-flux-induced flow’
$\boldsymbol{u}_{\mathit{DF}}$
, which scales as the amplitude squared of the waves. Its horizontal component,
$u_{\mathit{DF}}$
, is analogous to the pseudo-momentum per mass originally derived through the principle of wave action for internal waves (Bretherton Reference Bretherton1966, Reference Bretherton1969; Bretherton & Garrett Reference Bretherton and Garrett1969; Acheson Reference Acheson1976), later determined through the generalized Lagrangian mean formulation (Andrews & McIntyre Reference Andrews and McIntyre1978a
,Reference Andrews and McIntyre
b
; Bühler & McIntyre Reference Bühler and McIntyre1998) and Hamiltonian fluid dynamics (Scinocca & Shepherd Reference Scinocca and Shepherd1992), and, perhaps most physically intuitively, through momentum conservation for quasi-monochromatic wavepackets (Dosser & Sutherland Reference Dosser and Sutherland2011a
,Reference Dosser and Sutherland
b
), as reviewed in part in textbooks by Bühler (Reference Bühler2009, Reference Bühler2014) and Sutherland (Reference Sutherland2010).
For one-dimensional wavepackets – which are horizontally periodic (in the
$x$
-direction), spanwise infinite and vertically localized – the total wave-induced flow is just
$u_{\mathit{DF}}\hat{\boldsymbol{i}}$
, a non-divergent unidirectional flow in the
$x$
-direction. However, the structure of the wave-induced flow changes qualitatively for two- and three-dimensional wavepackets, for which the divergent velocity field
$\boldsymbol{u}_{\mathit{DF}}$
necessitates a response flow
$\boldsymbol{u}_{\mathit{RF}}$
, so that the total wave-induced flow
$\boldsymbol{u}_{\mathit{DF}}+\boldsymbol{u}_{\mathit{RF}}$
is non-divergent. Two-dimensional (spanwise-infinite) wavepackets induce long internal waves whose vertical phase speed matches the vertical group velocity of the wavepacket (Bretherton Reference Bretherton1969; Tabaei & Akylas Reference Tabaei and Akylas2007; van den Bremer & Sutherland Reference van den Bremer and Sutherland2014). The induced horizontal velocity field changes sign across the wavepacket itself, oriented with the horizontal group velocity along the leading flank and having opposite sign on the trailing flank. In contrast, three-dimensional wavepackets with comparable horizontal and vertical extents induce a horizontal circulation known as the ‘Bretherton flow’ (Bretherton Reference Bretherton1969; Bühler & McIntyre Reference Bühler and McIntyre2003; Tabaei & Akylas Reference Tabaei and Akylas2007; Bühler Reference Bühler2009). As for one-dimensional wavepackets, the wave-induced flow is strictly horizontal and oriented in the same direction as the horizontal component of the group velocity over the entire vertical extent of the wavepacket. However, the response flow sets up an oppositely oriented recirculation on the spanwise flanks of the wavepacket. For both two- and three-dimensional wavepackets, momentum transport is not necessarily localized to the wavepacket, as in the circumstance of ‘remote recoil’ by the Bretherton flow (Bühler & McIntyre Reference Bühler and McIntyre2003).
Unlike studies of wave beams (Kataoka & Akylas Reference Kataoka and Akylas2013, Reference Kataoka and Akylas2015), for which
$\unicode[STIX]{x1D714}$
is constant and the waves have infinite spatial extent in the along-beam direction (perpendicular to
$\boldsymbol{k}$
), we stress here that this study focuses upon wavepackets for which both the wavenumber and frequency have finite, but small, bandwidth.
The mean flow has important consequences for stability, which can dominate triad resonant interactions, including parametric subharmonic instability (Mied Reference Mied1976; Drazin Reference Drazin1977; Klostermeyer Reference Klostermeyer1991), if the wavetrain is sufficiently locally confined (e.g. Sutherland Reference Sutherland2006a
). For one-dimensional wavepackets, Akylas & Tabaei (Reference Akylas, Tabaei and Litvak2005) and Sutherland (Reference Sutherland2006b
) developed a nonlinear Schrödinger equation to examine the effect of mean flow on stability. In terms of
$\unicode[STIX]{x1D703}$
, the angle that lines of constant phase make with the vertical, Akylas & Tabaei (Reference Akylas, Tabaei and Litvak2005) and Sutherland (Reference Sutherland2006b
) predicted these one-dimensional wavepackets to be modulationally unstable if
$|\unicode[STIX]{x1D703}|<\unicode[STIX]{x1D703}_{c}$
, where
$\unicode[STIX]{x1D703}_{c}=\cos ^{-1}(2/3)\approx 35.3^{\circ }$
, corresponding to waves propagating at the maximum vertical group velocity. The prediction was confirmed using numerical simulations.
More generally, Tabaei & Akylas (Reference Tabaei and Akylas2007, henceforth TA07) examined the nonlinear interaction with the mean for three-dimensional wavepackets, distinguishing ‘round’ (localized in all three dimensions), ‘flat’ (wide in both horizontal directions compared to the vertical) and two-dimensional (spanwise-infinite) wavepackets in the presence of rotation. Without making the small-amplitude approximation a priori, TA07 performed perturbation expansions in all three cases to find coupled sets of differential equations describing the mean flow and the evolution of the packet, as influenced by the mean flow, and provided numerical solutions to these reduced-form equations. For horizontally periodic flat wavetrains and in the absence of rotation, TA07 recovered the stability criterion of Akylas & Tabaei (Reference Akylas, Tabaei and Litvak2005) and Sutherland (Reference Sutherland2006b ) discussed above. For round three-dimensional wavepackets, no resonance with the mean flow is found, and TA07 recover the equations derived previously by Shrira (Reference Shrira1981), with rotation not playing a leading-order role. However, performing a linear stability analysis, TA07 showed that flat three-dimensional wavepackets are always unstable. Crucially, they showed that the most unstable perturbation corresponds to a resonance with the hydrostatic gravity-inertial waves that constitutes the flow induced by the wavepacket.
Bridging the gap between one-, two- and three-dimensional wavepackets, here we examine flows induced by three-dimensional wavepackets of different aspect ratios. We take a step back from TA07: we do not examine the effect of the mean flow on nonlinear evolution of the packet and consider the response to a wavepacket that travels steadily without changing shape due to dispersion. Here, our goal is to derive explicit solutions for the flow induced by a wavepacket of any aspect ratio and thence to quantify how wide a wavepacket must be for induced long waves to dominate over the induced, horizontally recirculating, Bretherton flow and to quantify how long a wavepacket must be for the induced flow to be comparable to that for a one-dimensional wavepacket. These theoretical predictions are corroborated by numerical simulations of Gaussian wavepackets.
The paper is laid out as follows. After introducing the governing equations in § 2, perturbation theory for quasi-monochromatic wavepackets is applied in § 3. Therein, we generalize the results of Bretherton (Reference Bretherton1969) to wavepackets of arbitrary aspect ratio by identifying the divergent-flux-induced flow as the forcing of the total wave-induced flow and decomposing the latter into the sum of a local and horizontal Bretherton circulation and a non-local three-dimensional long-wave field. From asymptotic limits of these solutions, we construct a regime diagram that delineates for which wavepacket aspect ratios the induced flow is dominated by either the Bretherton flow or induced long waves. These predictions are compared with our fully nonlinear numerical simulations in § 4 and with the previous results of TA07 in § 5. Finally, conclusions are drawn in § 6.
2 Governing equations
Ignoring Coriolis, diffusion and viscous effects, and invoking the Boussinesq approximation, the conservation equations for momentum, internal energy and a statement of incompressibility are as follows:
$$\begin{eqnarray}\frac{\text{D}\boldsymbol{u}}{\text{D}t}=-\frac{1}{\unicode[STIX]{x1D70C}_{0}}\unicode[STIX]{x1D735}p-g\frac{\unicode[STIX]{x1D70C}}{\unicode[STIX]{x1D70C}_{0}}\hat{\boldsymbol{k}},\quad \frac{\text{D}\unicode[STIX]{x1D70C}}{\text{D}t}=-w\frac{\text{d}\overline{\unicode[STIX]{x1D70C}}}{\text{d}z},\quad \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0,\end{eqnarray}$$
where
$\boldsymbol{u}=(u,v,w)$
denotes the velocity vector with components in the directions of the respective unit vectors
$\hat{\boldsymbol{i}}$
(
$x$
),
$\hat{\boldsymbol{j}}$
(
$y$
) and
$\hat{\boldsymbol{k}}$
(
$z$
),
$p$
denotes the fluctuation pressure, and gravity
$g$
acts in the negative
$z$
direction. In the Boussinesq approximation,
$\unicode[STIX]{x1D70C}_{0}$
is the (constant) characteristic density. The background density is
$\overline{\unicode[STIX]{x1D70C}}(z)$
and, assuming a uniformly stratified ambient, the squared buoyancy frequency
$N^{2}=-(g/\unicode[STIX]{x1D70C}_{0})(\text{d}\overline{\unicode[STIX]{x1D70C}}/\text{d}z)$
is real and constant. It is convenient to express the density perturbation
$\unicode[STIX]{x1D70C}$
in terms of the vertical displacement field
$\unicode[STIX]{x1D709}$
through
$\unicode[STIX]{x1D709}\equiv -\unicode[STIX]{x1D70C}/\overline{\unicode[STIX]{x1D70C}}^{\prime }(z)$
, where the prime denotes a
$z$
derivative. Then (2.1b
) can be rewritten as
$\text{D}\unicode[STIX]{x1D709}/\text{D}t=w$
. By taking the curl of the momentum equation (2.1a
), the pressure terms can be eliminated and we obtain an equation for the baroclinic generation of vorticity
$\unicode[STIX]{x1D73B}\equiv \unicode[STIX]{x1D735}\times \boldsymbol{u}$
:
$$\begin{eqnarray}\frac{\text{D}\unicode[STIX]{x1D73B}}{\text{D}t}=\unicode[STIX]{x1D73B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}-N^{2}\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}{\unicode[STIX]{x2202}y}\hat{\boldsymbol{i}}-\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}{\unicode[STIX]{x2202}x}\hat{\boldsymbol{j}}\right).\end{eqnarray}$$
Equations (2.1b
) and (2.2) can be combined, by taking
$x$
and
$y$
derivatives of (2.1b
). Substituting the result into (2.2) after taking its temporal derivative gives
$$\begin{eqnarray}\displaystyle & & \displaystyle \underbrace{\left[\begin{array}{@{}ccc@{}}0 & -\unicode[STIX]{x2202}_{t}^{2}\unicode[STIX]{x2202}_{z} & (\unicode[STIX]{x2202}_{t}^{2}+N^{2})\unicode[STIX]{x2202}_{y}\\ \unicode[STIX]{x2202}_{t}^{2}\unicode[STIX]{x2202}_{z} & 0 & -(\unicode[STIX]{x2202}_{t}^{2}+N^{2})\unicode[STIX]{x2202}_{x}\\ -\unicode[STIX]{x2202}_{t}^{2}\unicode[STIX]{x2202}_{y} & \unicode[STIX]{x2202}_{t}^{2}\unicode[STIX]{x2202}_{x} & 0\end{array}\right]}_{\equiv \,\unicode[STIX]{x1D647}}\boldsymbol{u}\nonumber\\ \displaystyle & & \displaystyle \quad =\underbrace{\overline{\unicode[STIX]{x1D735}\boldsymbol{\cdot }[-\unicode[STIX]{x2202}_{t}(\boldsymbol{u}\otimes \unicode[STIX]{x1D73B})-N^{2}[\unicode[STIX]{x2202}_{x}(\boldsymbol{u}\unicode[STIX]{x1D709})\hat{\boldsymbol{j}}-\unicode[STIX]{x2202}_{y}(\boldsymbol{u}\unicode[STIX]{x1D709})\hat{\boldsymbol{i}}]]+\unicode[STIX]{x2202}_{t}(\unicode[STIX]{x1D73B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u})}}_{\equiv \,\boldsymbol{F}},\end{eqnarray}$$
where the linear matrix operator
$\unicode[STIX]{x1D647}$
acts on the velocity vector
$\boldsymbol{u}$
and the terms in the right-hand side vector
$\boldsymbol{F}$
are nonlinear in amplitude. The symbol
$\otimes$
denotes the tensor product, and the overline denotes averaging over the wave scales, as our interest here is in mean-flow terms only. The components of (2.3) are written explicitly in appendix A. We will refer to (2.3) as the forcing equation for the induced mean flow, as linear (in amplitude) solutions substituted into the right-hand side will, after averaging, give rise to a forcing of an order amplitude-squared mean flow on the left-hand side.
3 Perturbation theory
Here we derive the perturbation equations and their solutions, using Fourier transforms, for quasi-monochromatic internal wavepackets that are localized in all three dimensions, but whose amplitude envelope has arbitrary aspect ratio. Without loss of generality, it is assumed that the wavepacket travels in the
$x$
–
$z$
plane, for which the waves with wavenumber vector
$\boldsymbol{k}=(k_{x},0,k_{z})$
are contained within an amplitude envelope having extent
$\unicode[STIX]{x1D70E}_{x}$
,
$\unicode[STIX]{x1D70E}_{y}$
and
$\unicode[STIX]{x1D70E}_{z}$
in the horizontal ‘along-wave’, spanwise and vertical directions, respectively. The waves have a peak frequency given by the dispersion relation for non-rotating Boussinesq internal waves,
$\unicode[STIX]{x1D714}=Nk_{x}/|\boldsymbol{k}|$
.
This section begins in § 3.1 with the definition of the scalings used to define the quasi-monochromatic wavepacket and the presentation of its polarization relations, taking into account its finite spatial extent. Realizing that the solution of (2.3) for general aspect ratio results in both a purely horizontal Bretherton flow and long waves, we formalize this decomposition of the solution in § 3.2 and solve for these two components of the mean flow separately in §§ 3.3 and 3.4, respectively. Finally, specifically for Gaussian wavepackets, relationships are derived that delineate the different regimes of behaviour of the mean flow in § 3.5.
3.1 Quasi-monochromatic wavepackets
We introduce the slow variables
$X$
,
$Y$
and
$Z$
describing the translation of the wavepacket at the group velocity,
$\boldsymbol{c}_{g}\equiv (c_{gx},0,c_{gz})=N(k_{z}^{2},0,-k_{x}k_{z})/|\boldsymbol{k}|^{3}$
:
$$\begin{eqnarray}X=\unicode[STIX]{x1D700}_{x}(x-c_{gx}t),\quad Y=\unicode[STIX]{x1D700}_{y}y,\quad Z=\unicode[STIX]{x1D700}_{z}(z-c_{gz}t),\end{eqnarray}$$
in which the small bandwidth parameters based on the spatial scales of variation of the amplitude envelope are defined as
$\unicode[STIX]{x1D700}_{x}\equiv 1/(k_{x}\unicode[STIX]{x1D70E}_{x})$
,
$\unicode[STIX]{x1D700}_{y}\equiv 1/(k_{x}\unicode[STIX]{x1D70E}_{y})$
and
$\unicode[STIX]{x1D700}_{z}\equiv 1/(k_{x}\unicode[STIX]{x1D70E}_{z})$
. Dispersion of the group occurs on the slow time scale
$T$
, which is understood to be much slower than that which drives the wave-induced flow. Effectively, dispersion is ignored. Formally, we define
$\unicode[STIX]{x1D700}\equiv \max \{\unicode[STIX]{x1D700}_{x},\unicode[STIX]{x1D700}_{y},\unicode[STIX]{x1D700}_{z}\}$
and
$T\equiv \unicode[STIX]{x1D700}^{2}t$
. As a consequence, time derivatives on the right-hand side of (2.3), which act upon the amplitude envelope and not the waves themselves, can be represented at leading order by spatial derivatives on the envelope scale denoting the translation of the wavepacket:
$$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\sim -\unicode[STIX]{x1D700}_{x}c_{gx}\unicode[STIX]{x2202}_{X}-\unicode[STIX]{x1D700}_{z}c_{gz}\unicode[STIX]{x2202}_{Z}.\end{eqnarray}$$
Table 1. Expressions for different linear (in
$\unicode[STIX]{x1D6FC}$
) fields (first column) at the scale of the waves (second column) and with small corrections (third column). This last column accounts for the finite but relatively large extent of the wavepacket in the
$x$
,
$y$
and
$z$
directions, as prescribed through the values of
$\unicode[STIX]{x1D700}_{x}$
,
$\unicode[STIX]{x1D700}_{y}$
and
$\unicode[STIX]{x1D700}_{z}$
, respectively (captured by
$\unicode[STIX]{x1D700}_{i}$
in the column name). Values are given in terms of the amplitude envelope
$A$
of the vertical displacement field, as found through the linearized equations for internal waves propagating in the
$x$
–
$z$
plane. It is understood that the actual fields are the real parts of the tabulated expressions, and the phase is given by
$\unicode[STIX]{x1D711}=k_{x}x+k_{z}z-\unicode[STIX]{x1D714}t$
.
We introduce a second small parameter
$\unicode[STIX]{x1D6FC}=|k_{x}|A_{0}$
, in which
$A_{0}$
is the amplitude of vertical displacements at the centre of the group. At first order in
$\unicode[STIX]{x1D6FC}$
, the vertical displacement field is given by
$$\begin{eqnarray}\unicode[STIX]{x1D709}^{(1)}=\unicode[STIX]{x1D709}_{0}^{(1)}=\text{Re}[A(X,Y,Z,T)\text{e}^{\text{i}(k_{x}x+k_{z}z-\unicode[STIX]{x1D714}t)}],\end{eqnarray}$$
in which the superscript denotes the order in
$\unicode[STIX]{x1D6FC}$
and the subscript denotes the order in the bandwidth parameters
$\unicode[STIX]{x1D700}_{x}$
,
$\unicode[STIX]{x1D700}_{y}$
and
$\unicode[STIX]{x1D700}_{z}$
. As
$\unicode[STIX]{x1D700}_{x}$
,
$\unicode[STIX]{x1D700}_{y}$
and
$\unicode[STIX]{x1D700}_{z}$
can be of different orders depending on the aspect ratio of the wavepacket that we later choose, we emphasize that the subscript does not play the role of a formal ordering parameter, but simply denotes the leading-order correction for the packet structure of the waves. The operator
$\text{Re}[~]$
denotes taking the real part of the enclosed expression. For the vertical displacement field, we set higher-order terms in the bandwidth parameters to be zero without loss of generality. The second column of table 1 reports the usual polarization relationships for plane periodic internal waves, and the third column gives the next-order corrections accounting for the finite size of the wavepacket, which can be obtained from solving the linearized (in
$\unicode[STIX]{x1D6FC}$
) governing equations (2.1).
We also define the aspect ratios
$R_{y}\equiv \unicode[STIX]{x1D70E}_{y}/\unicode[STIX]{x1D70E}_{x}$
and
$R_{z}\equiv \unicode[STIX]{x1D70E}_{z}/\unicode[STIX]{x1D70E}_{x}$
. Table 2 formalizes the terminology used throughout the paper, distinguishing definitions based on the geometry of the wavepacket itself and definitions based on the nature and geometry of the induced flow. We will discuss the latter categorization in § 3.5 when we propose our regime diagrams. In the most general case of a three-dimensional (3D) wavepacket, the linear waves are modulated in all three directions, whereas a two-dimensional (2D) wavepacket is spanwise independent, as obtained by setting
$\unicode[STIX]{x1D70E}_{y}\rightarrow \infty$
and thus
$\unicode[STIX]{x1D700}_{y}\rightarrow 0$
. A one-dimensional (1D) wavepacket is not only spanwise independent or infinitely wide, but also periodic in the
$x$
-direction, resulting in
$\unicode[STIX]{x1D70E}_{x}\rightarrow \infty$
and thus
$\unicode[STIX]{x1D700}_{x}\rightarrow 0$
. Such a 1D wavepacket is only modulated in the
$z$
-direction. Approximately then, a 3D wavepacket is said to be ‘round’ when all three of its spatial scales are comparable:
$0<\unicode[STIX]{x1D700}_{x}\sim \unicode[STIX]{x1D700}_{y}\sim \unicode[STIX]{x1D700}_{z}\ll 1$
. A 3D wavepacket is said to be ‘wide’ when its spanwise extent is larger than its ‘along-wave’ and vertical extents:
$0<\unicode[STIX]{x1D700}_{y}\ll \unicode[STIX]{x1D700}_{x}\sim \unicode[STIX]{x1D700}_{z}\ll 1$
. Finally, a 3D wavepacket is said to be ‘flat’ when both its horizontal scales are comparable and larger than its vertical scale:
$0<\unicode[STIX]{x1D700}_{x}\sim \unicode[STIX]{x1D700}_{y}\ll \unicode[STIX]{x1D700}_{z}\ll 1$
, corresponding to the definition in TA07.
Table 2. Wavepacket- and induced-flow-based definitions. Note that the induced-flow-based definitions are subject to generalizations and restrictions discussed in § 3.5.
3.2 Mean-flow decomposition: Bretherton flow and long waves
It is evident from the outset that the solution to the mean-flow forcing equation (2.3) in combination with conservation of volume (2.1c ), a total of four equations, can give rise to three types of induced flows, leading to the classification into wavepackets that are effectively 1D (primarily inducing unidirectional flows), 2D (primarily inducing long-waves) and fully 3D (primarily inducing horizontally recirculating flows) (cf. table 2). The relative magnitudes of the three types of mean flow depend on the aspect ratios of the wavepackets. We first conceptually distinguish two contributions to the total mean flow of a wavepacket of general aspect ratio, namely the Bretherton flow (BF) and long waves (LW):
$$\begin{eqnarray}\boldsymbol{u}^{(2)}=\boldsymbol{u}_{\mathit{BF}}+\boldsymbol{u}_{\mathit{LW}}.\end{eqnarray}$$
Writing the forcing in (2.3) as
$\boldsymbol{F}=(F_{x},F_{y},F_{z})$
, we define the two components in (3.4) explicitly, so that
$$\begin{eqnarray}\unicode[STIX]{x1D647}\boldsymbol{u}_{\mathit{BF}}=\left[\begin{array}{@{}c@{}}0\\ 0\\ F_{z}\end{array}\right]\quad \text{and}\quad \unicode[STIX]{x1D647}\boldsymbol{u}_{\mathit{LW}}=\left[\begin{array}{@{}c@{}}F_{x}\\ F_{y}\\ 0\end{array}\right],\end{eqnarray}$$
where the linear operator
$\unicode[STIX]{x1D647}$
is defined in (2.3). Both components separately satisfy conservation of volume:
$\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{\mathit{BF}}=0$
and
$\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{\mathit{LW}}=0$
.
We begin in § 3.3 by deriving solutions for the well-known Bretherton flow, and consider the induced long waves in § 3.4. The method of inverse Fourier transforming, taking appropriate branch cuts to remove the singularity in the integrand, is also described in § 3.4. The total wave-induced mean flow is then given by the linear superposition of the two contributions, as defined in (3.4). We emphasize that we only consider up to second order in wave steepness
$\unicode[STIX]{x1D6FC}$
and leading order in the separation-of-scales parameter
$\unicode[STIX]{x1D700}$
, so that the effect of the mean flow on the nonlinear evolution of the wavepacket and the change in shape of the wavepacket due to dispersion, as considered by TA07, are beyond our scope. We return to this in § 5.
3.3 Bretherton flow
Taking the component of the total induced flow that is forced by the vertical vorticity of the linear waves (cf. (3.5a
)), the Bretherton flow
$\boldsymbol{u}_{\mathit{BF}}$
, we further decompose this flow as the sum of the divergent-flux-induced flow
$\boldsymbol{u}_{\mathit{DF}}$
and the response flow
$\boldsymbol{u}_{\mathit{RF}}$
:
$$\begin{eqnarray}\boldsymbol{u}_{\mathit{BF}}=\boldsymbol{u}_{\mathit{DF}}+\boldsymbol{u}_{\mathit{RF}}.\end{eqnarray}$$
The divergent-flux-induced flow is determined explicitly from the advection terms in the momentum equations (§ 3.3.1). The response flow is then found from conditions for incompressibility and irrotational flow (§ 3.3.2).
3.3.1 Divergent-flux-induced flow
Generally, the divergent-flux-induced flow results from the acceleration of the flow resulting from the slow spatial variations in the advection terms of the momentum equations (2.1a ) written in flux form (van den Bremer & Sutherland Reference van den Bremer and Sutherland2014):
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\boldsymbol{u}_{\mathit{DF}}}{\unicode[STIX]{x2202}t}\equiv \overline{-\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\boldsymbol{u}^{(1)}\otimes \boldsymbol{u}^{(1)})},\end{eqnarray}$$
in which the overline denotes averaging over the short spatial scales of the individual waves. At leading order in
$\unicode[STIX]{x1D700}$
, the time evolution of the wavepacket may be recast in terms of spatial derivatives according to (3.2). Hence, we have
$$\begin{eqnarray}\boldsymbol{c}_{g}u_{\mathit{DF}}=\overline{u^{(1)}\boldsymbol{u}^{(1)}},\quad v_{\mathit{DF}}=0,\quad \boldsymbol{c}_{g}w_{\mathit{DF}}=\overline{w^{(1)}\boldsymbol{u}^{(1)}}.\end{eqnarray}$$
Using the polarization relationships in table 1, we obtain from (3.8) the divergent-flux-induced flow, which is oriented in the direction of the group velocity vector
$\boldsymbol{c}_{g}$
(van den Bremer & Sutherland Reference van den Bremer and Sutherland2014):
$$\begin{eqnarray}\boldsymbol{u}_{\mathit{DF}}=u_{\mathit{DF}}\left(1,0,\frac{c_{gz}}{c_{gx}}\right),\end{eqnarray}$$
in which
$$\begin{eqnarray}u_{\mathit{DF}}=\text{sgn}(k_{x}){\textstyle \frac{1}{2}}N|\boldsymbol{k}||A|^{2}.\end{eqnarray}$$
The horizontal component of the divergent-flux-induced flow
$u_{\mathit{DF}}$
is the well-established induced mean flow for 1D wavepackets with amplitude envelope
$A(Z,T)$
except that in (3.10) the amplitude envelope is generally a function of 3D space and time:
$A\equiv A(X,Y,Z,T)$
. The quantity
$u_{\mathit{DF}}$
is equal to the pseudo-momentum per mass, which historically has been represented in terms of the wave action (Bretherton Reference Bretherton1969; Acheson Reference Acheson1976) by
$\overline{E}k_{x}/\unicode[STIX]{x1D714}$
, in which
$\overline{E}$
is the mean energy per unit mass, or by an expression derived from Hamiltonian fluid mechanics (Scinocca & Shepherd Reference Scinocca and Shepherd1992),
$-\overline{\unicode[STIX]{x1D709}\unicode[STIX]{x1D701}_{y}}$
, in which
$\unicode[STIX]{x1D701}_{y}$
is the spanwise vorticity (see also the textbooks by Bühler (Reference Bühler2009, Reference Bühler2014) and Sutherland (Reference Sutherland2010)).
3.3.2 Response and total Bretherton flow
At leading order, we obtain from the first two rows of (3.5a
) that
$\unicode[STIX]{x2202}_{Y}w_{\mathit{BF}}(X,Y,Z)=0$
and
$\unicode[STIX]{x2202}_{X}w_{\mathit{BF}}(X,Y,Z)=0$
, respectively. Hence
$w_{\mathit{BF}}(X,Y,Z)=0$
is the only meaningful solution: the Bretherton flow is purely horizontal. Although arising by assumption here, this corresponds to the classical result first derived by Bretherton (Reference Bretherton1969), and subsequently reproduced and used in studies of long-range influence and transfers of energy to large scales by localized wavepackets (Bühler & McIntyre Reference Bühler and McIntyre1998, Reference Bühler and McIntyre2003; Tabaei & Akylas Reference Tabaei and Akylas2007; Wagner & Young Reference Wagner and Young2015; Xie & Vanneste Reference Xie and Vanneste2015). With
$\boldsymbol{u}_{\mathit{DF}}$
given by (3.9) and the requirement that the flow is purely horizontal (
$w_{\mathit{BF}}=0$
) following from the first two rows of (3.5a
), we must have
$w_{\mathit{RF}}=-w_{\mathit{DF}}$
. The non-zero forcing arises in the third equation of (3.5a
). Substituting the polarization relationships in table 1 into the right-hand side (given explicitly in (A 3)), it is apparent that
$(\unicode[STIX]{x1D735}\times \boldsymbol{u}_{\mathit{BF}})\boldsymbol{\cdot }\hat{\boldsymbol{k}}=(\unicode[STIX]{x1D735}\times \boldsymbol{u}_{\mathit{DF}})\boldsymbol{\cdot }\hat{\boldsymbol{k}}$
. In other words, the response flow
$\boldsymbol{u}_{\mathit{RF}}$
is irrotational in the horizontal, and all the vertical vorticity in the mean flow
$\boldsymbol{u}_{\mathit{BF}}$
is derived from the divergent-flux-induced flow (cf. (3.5a
)), so that
$(\unicode[STIX]{x1D735}\times \boldsymbol{u}_{\mathit{RF},H})\boldsymbol{\cdot }\hat{\boldsymbol{k}}=0$
. Noting that
$v_{\mathit{DF}}=0$
in (3.9), the latter condition is recast in terms of the vertical vorticity of the total horizontal flow being driven by the divergent-flux-induced flow:
$$\begin{eqnarray}\unicode[STIX]{x2202}_{x}v_{\mathit{BF}}-\unicode[STIX]{x2202}_{y}u_{\mathit{BF}}=-\unicode[STIX]{x2202}_{y}u_{\mathit{DF}}.\end{eqnarray}$$
Finally, a fourth equation (in addition to the three in (3.5a )) requires the horizontal Bretherton flow to satisfy incompressibility, so that
$$\begin{eqnarray}\unicode[STIX]{x2202}_{x}u_{\mathit{BF}}+\unicode[STIX]{x2202}_{y}v_{\mathit{BF}}=0.\end{eqnarray}$$
We solve these equations using 3D Fourier transforms with respect to the unscaled translating coordinates
$(\tilde{x},{\tilde{y}},\tilde{z})=(x-c_{gx}t,y,z-c_{gz}t)$
to give transformed velocities in terms of the wavenumber vector
$(\unicode[STIX]{x1D705},\unicode[STIX]{x1D706},\unicode[STIX]{x1D707})$
. In particular, the Fourier transform of the divergent-flux-induced flow given by (3.10) is
$$\begin{eqnarray}\widehat{u_{\mathit{DF}}}=\text{sgn}(k_{x})\frac{1}{8\unicode[STIX]{x03C0}^{3}}\left(\frac{1}{2}N|\boldsymbol{k}|\right)\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }|A(\tilde{\boldsymbol{x}})|^{2}\text{e}^{-\text{i}(\unicode[STIX]{x1D705}\tilde{x}+\unicode[STIX]{x1D706}{\tilde{y}}+\unicode[STIX]{x1D707}\tilde{z})}\,\text{d}\tilde{x}\,\text{d}{\tilde{y}}\,\text{d}\tilde{z}.\end{eqnarray}$$
The algebraic equations resulting from Fourier-transforming (3.12) and (3.11) are solved to give the total induced flow in Fourier space:
$$\begin{eqnarray}\widehat{\boldsymbol{u}_{\mathit{BF}}}=\widehat{u_{\mathit{DF}}}\left[\begin{array}{@{}c@{}}\unicode[STIX]{x1D706}^{2}/\unicode[STIX]{x1D705}_{H}^{2}\\ -\unicode[STIX]{x1D705}\unicode[STIX]{x1D706}/\unicode[STIX]{x1D705}_{H}^{2}\\ 0\end{array}\right],\end{eqnarray}$$
in which
$\unicode[STIX]{x1D705}_{H}=(\unicode[STIX]{x1D705}^{2}+\unicode[STIX]{x1D706}^{2})^{1/2}$
. The inverse transform gives the induced flow in real space:
$$\begin{eqnarray}\boldsymbol{u}_{\mathit{BF}}(\tilde{\boldsymbol{x}})=\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }\widehat{u_{\mathit{DF}}}\left[\begin{array}{@{}c@{}}\unicode[STIX]{x1D706}^{2}/\unicode[STIX]{x1D705}_{H}^{2}\\ -\unicode[STIX]{x1D705}\unicode[STIX]{x1D706}/\unicode[STIX]{x1D705}_{H}^{2}\\ 0\end{array}\right]\text{e}^{\text{i}(\unicode[STIX]{x1D705}\tilde{x}+\unicode[STIX]{x1D706}{\tilde{y}}+\unicode[STIX]{x1D707}\tilde{z})}\,\text{d}\unicode[STIX]{x1D705}\,\text{d}\unicode[STIX]{x1D706}\,\text{d}\unicode[STIX]{x1D707}.\end{eqnarray}$$
Of course, if one can separate the
$\unicode[STIX]{x1D707}$
dependence of
$\widehat{u_{\mathit{DF}}}$
from its
$\unicode[STIX]{x1D705}$
and
$\unicode[STIX]{x1D706}$
dependences, then the integral in
$\unicode[STIX]{x1D707}$
explicitly extracts the
$\tilde{z}$
-dependent part of
$u_{\mathit{DF}}$
,
$N|\boldsymbol{k}||{\mathcal{A}}(\tilde{z})|^{2}/2$
, from the integrals in
$\unicode[STIX]{x1D705}$
and
$\unicode[STIX]{x1D706}$
. The remaining double Fourier integral is the solution of the partial differential equation given by (22) in TA07.
3.4 Long waves
Here, we consider the component of the total induced flow that is forced by the horizontal vorticity of the linear waves (cf. (3.5b
)). We proceed with the following two assumptions. First, we assume that the induced flow has a much larger
$x$
extent than that of the wavepacket itself motivated by considerations similar to the 2D (spanwise-infinite) case (cf. Bretherton Reference Bretherton1969; Tabaei & Akylas Reference Tabaei and Akylas2007; van den Bremer & Sutherland Reference van den Bremer and Sutherland2014). Explicitly, we let the induced long waves on the left-hand side vary horizontally on the new coordinate
$\tilde{X}=\unicode[STIX]{x1D700}^{2}(x-c_{gx}t)$
. For a 2D wavepacket, Bretherton (Reference Bretherton1969) and van den Bremer & Sutherland (Reference van den Bremer and Sutherland2014) replace the forcing on the right-hand side by a Dirac delta function at the centre of the packet of equivalent strength by integrating with respect to
$\tilde{x}$
between
$\tilde{x}\rightarrow -\infty$
and
$\tilde{x}\rightarrow \infty$
. As a consequence, any terms in the forcing that involve slow
$x$
derivatives do not contribute. Second, we assume the wavepacket itself is relatively wide,
$\unicode[STIX]{x1D700}_{y}\sim \unicode[STIX]{x1D700}^{2}$
, because a non-negligibly small contribution from the long waves is only expected then. Explicitly, we set
$\unicode[STIX]{x1D700}_{z}=\unicode[STIX]{x1D700}$
and suppose
$\unicode[STIX]{x1D700}_{y}\sim \unicode[STIX]{x1D700}^{2}$
in the linear operator on the left-hand side and in the forcing on the right-hand side of (2.3). Combining both assumptions, the wavepacket we examine is thus wide, but not necessarily flat (cf. TA07), as we set
$\unicode[STIX]{x1D700}_{x}=\unicode[STIX]{x1D700}$
(or smaller) in the forcing on the right-hand side. From (2.3) and (3.5b
) we obtain
$$\begin{eqnarray}\left[\begin{array}{@{}ccc@{}}0 & -c_{gz}^{2}\unicode[STIX]{x1D700}^{3}\unicode[STIX]{x2202}_{ZZZ} & \unicode[STIX]{x1D700}^{4}c_{gz}^{2}\unicode[STIX]{x2202}_{ZZY}+N^{2}\unicode[STIX]{x1D700}^{2}\unicode[STIX]{x2202}_{Y}\\ c_{gz}^{2}\unicode[STIX]{x1D700}^{3}\unicode[STIX]{x2202}_{ZZZ} & 0 & -\unicode[STIX]{x1D700}^{4}c_{gz}^{2}\unicode[STIX]{x2202}_{ZZ\tilde{X}}-N^{2}\unicode[STIX]{x1D700}^{2}\unicode[STIX]{x2202}_{\tilde{X}}\\ -c_{gz}^{2}\unicode[STIX]{x1D700}^{4}\unicode[STIX]{x2202}_{ZZY} & c_{gz}^{2}\unicode[STIX]{x1D700}^{4}\unicode[STIX]{x2202}_{ZZ\tilde{X}} & 0\end{array}\right]\boldsymbol{u}_{\mathit{LW}}(\tilde{X},Y,Z)=\left[\begin{array}{@{}c@{}}F_{x}\\ F_{y}\\ 0\end{array}\right],\end{eqnarray}$$
where we note that the averaging in (2.3) corresponds to averaging over the fast scales of waves within the wavepacket. It is evident that the entries in the final row, corresponding to operations acting on the vertical vorticity field, are of order
$\unicode[STIX]{x1D700}^{4}$
. In order to have as many equations as unknowns after truncation at
$O(\unicode[STIX]{x1D700}^{3})$
, we replace the final row by the incompressibility condition
$\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{\mathit{LW}}=0$
, which still needs to be satisfied.
The polarization relationships in table 1 allow for explicit evaluation of
$(F_{x},F_{y},0)$
. From our first assumption above, we neglect terms involving slow
$x$
derivatives because the response to the forcing is long compared to the wavepacket. We avoid the cumbersome operation of introducing Dirac delta functions, since neglecting terms involving slow
$x$
derivatives leads to equivalent results. From our second assumption above, we neglect terms involving slow
$y$
derivatives because the wavepacket is relatively wide. The terms involving
$Z$
derivatives alone come from the product of the vertical velocity with the spanwise vorticity, each occurring at order
$\unicode[STIX]{x1D6FC}$
. Thus, we find that the dominant forcing at
$O(\unicode[STIX]{x1D6FC}^{2}\unicode[STIX]{x1D700}^{3})$
is
$$\begin{eqnarray}F_{x}=0\quad \text{and}\quad F_{y}=\frac{N^{3}k_{x}^{2}k_{z}^{2}}{2|\boldsymbol{k}|^{5}}\unicode[STIX]{x2202}_{\tilde{z}\tilde{z}\tilde{z}}|A(\tilde{x},{\tilde{y}},\tilde{z})|^{2}.\end{eqnarray}$$
Combining these results, retaining terms up to
$O(\unicode[STIX]{x1D700}^{3})$
, recasting the result in terms of the unscaled translating coordinate
$(\tilde{x},{\tilde{y}},\tilde{z})$
and using (3.10), we have the following:
$$\begin{eqnarray}\left[\begin{array}{@{}ccc@{}}0 & -c_{gz}^{2}\unicode[STIX]{x2202}_{\tilde{z}\tilde{z}\tilde{z}} & N^{2}\unicode[STIX]{x2202}_{{\tilde{y}}}\\ c_{gz}^{2}\unicode[STIX]{x2202}_{\tilde{z}\tilde{z}\tilde{z}} & 0 & -N^{2}\unicode[STIX]{x2202}_{\tilde{x}}\\ \unicode[STIX]{x2202}_{\tilde{x}} & \unicode[STIX]{x2202}_{{\tilde{y}}} & \unicode[STIX]{x2202}_{\tilde{z}}\end{array}\right]\boldsymbol{u}_{\mathit{LW}}=\left[\begin{array}{@{}c@{}}0\\ c_{gz}^{2}\unicode[STIX]{x2202}_{\tilde{z}\tilde{z}\tilde{z}}u_{\mathit{DF}}\\ 0\end{array}\right].\end{eqnarray}$$
In the expression for the forcing on the right-hand side, we have used the equation for the vertical group velocity,
$c_{gz}=-Nk_{x}k_{z}/|\boldsymbol{k}|^{3}$
. Next, we obtain the Fourier transform of (3.18), solve the resulting algebraic equations and then obtain the inverse Fourier transform. This gives the following general integral equation for the long-wave contribution to the wave-induced flow:
$$\begin{eqnarray}\boldsymbol{u}_{\mathit{LW}}=\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }\frac{\widehat{u_{\mathit{DF}}}(\unicode[STIX]{x1D705},\unicode[STIX]{x1D706},\unicode[STIX]{x1D707})}{c_{gz}^{2}\unicode[STIX]{x1D707}^{4}-N^{2}\unicode[STIX]{x1D705}_{H}^{2}}\left[\begin{array}{@{}c@{}}c_{gz}^{2}\unicode[STIX]{x1D707}^{4}-N^{2}\unicode[STIX]{x1D706}^{2}\\ N^{2}\unicode[STIX]{x1D705}\unicode[STIX]{x1D706}\\ -c_{gz}^{2}\unicode[STIX]{x1D705}\unicode[STIX]{x1D707}^{3}\end{array}\right]\text{e}^{\text{i}(\unicode[STIX]{x1D705}\tilde{x}+\unicode[STIX]{x1D706}{\tilde{y}}+\unicode[STIX]{x1D707}\tilde{z})}\,\text{d}\unicode[STIX]{x1D705}\,\text{d}\unicode[STIX]{x1D706}\,\text{d}\unicode[STIX]{x1D707}.\end{eqnarray}$$
The singularity in the integrand of (3.19) corresponds to hydrostatic internal waves with wavenumber
$(\unicode[STIX]{x1D705},\unicode[STIX]{x1D706},\unicode[STIX]{x1D707})$
set so that the vertical phase speed of the induced (long internal) waves equals the vertical group velocity of the wavepacket. It is the singularity in (3.19) that results in an induced flow with non-negligible vertical velocity. In order to select only outgoing waves, we make use of Cauchy’s residue theorem to integrate around the singularity in (3.19) and reduce the expression to a formula involving only a double integral (see also van den Bremer & Sutherland Reference van den Bremer and Sutherland2014). In order to apply the correct branch cut that selects only the outgoing waves, we switch to cylindrical coordinates in wavenumber space, defining
$\unicode[STIX]{x1D703}$
so that
$(\unicode[STIX]{x1D705},\unicode[STIX]{x1D706})=\unicode[STIX]{x1D705}_{H}(\cos \unicode[STIX]{x1D703},\sin \unicode[STIX]{x1D703})$
. Equation (3.19) becomes
$$\begin{eqnarray}\boldsymbol{u}_{\mathit{LW}}=\int _{-\infty }^{\infty }\int _{0}^{2\unicode[STIX]{x03C0}}\!\!\int _{0}^{\infty }\frac{\widehat{u_{\mathit{DF}}}(\unicode[STIX]{x1D705}_{H},\unicode[STIX]{x1D703},\unicode[STIX]{x1D707})}{c_{gz}^{2}\unicode[STIX]{x1D707}^{4}-N^{2}\unicode[STIX]{x1D705}_{H}^{2}}\left[\begin{array}{@{}c@{}}c_{gz}^{2}\unicode[STIX]{x1D707}^{4}-N^{2}\unicode[STIX]{x1D705}_{H}^{2}\sin ^{2}(\unicode[STIX]{x1D703})\\ N^{2}\unicode[STIX]{x1D705}_{H}^{2}\cos (\unicode[STIX]{x1D703})\sin (\unicode[STIX]{x1D703})\\ -c_{gz}^{2}\unicode[STIX]{x1D707}^{3}\unicode[STIX]{x1D705}_{H}\cos (\unicode[STIX]{x1D703})\end{array}\right]\text{e}^{\text{i}(\unicode[STIX]{x1D705}_{H}\tilde{r}\cos (\unicode[STIX]{x1D703}-\unicode[STIX]{x1D703}_{p})+\unicode[STIX]{x1D707}\tilde{z})}\unicode[STIX]{x1D705}_{H}\,\text{d}\unicode[STIX]{x1D705}_{H}\,\text{d}\unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D707},\end{eqnarray}$$
where
$\tilde{r}\equiv \sqrt{\tilde{x}^{2}+{\tilde{y}}^{2}}$
and
$\unicode[STIX]{x1D703}_{p}\equiv \tan ^{-1}({\tilde{y}}/\tilde{x})$
. The singularity in (3.20) is split into two terms:
$$\begin{eqnarray}G\equiv \frac{1}{c_{gz}^{2}\unicode[STIX]{x1D707}^{4}-N^{2}\unicode[STIX]{x1D705}_{H}^{2}}=\frac{1}{2c_{gz}\unicode[STIX]{x1D707}^{2}}\left[\frac{1}{N\unicode[STIX]{x1D705}_{H}+c_{gz}\unicode[STIX]{x1D707}^{2}}-\frac{1}{N\unicode[STIX]{x1D705}_{H}-c_{gz}\unicode[STIX]{x1D707}^{2}}\right].\end{eqnarray}$$
This expression in the integral can be decomposed into the sum of the Cauchy principal value and a Dirac delta function (see also Voisin Reference Voisin1991):
$$\begin{eqnarray}G=G_{PV}-\frac{\text{i}2\unicode[STIX]{x03C0}}{2c_{gz}\unicode[STIX]{x1D707}^{2}}\unicode[STIX]{x1D6FF}(N\unicode[STIX]{x1D705}_{H}-|c_{gz}|\unicode[STIX]{x1D707}^{2})\,\text{sgn}(\unicode[STIX]{x1D707}),\end{eqnarray}$$
where the absolute value and sign of
$\unicode[STIX]{x1D707}$
ensure outward-propagating waves. It is the delta function that gives rise to the radiating waves in the solution of the integral with respect to
$\unicode[STIX]{x1D705}_{H}$
in (3.20):
$$\begin{eqnarray}\displaystyle \boldsymbol{u}_{\mathit{LW}} & = & \displaystyle -\frac{\text{i}\unicode[STIX]{x03C0}\,\text{sgn}(c_{gz})}{N^{2}}\int _{-\infty }^{\infty }\int _{0}^{2\unicode[STIX]{x03C0}}\widehat{u_{\mathit{DF}}}\left[\begin{array}{@{}c@{}}c_{gz}^{2}\unicode[STIX]{x1D707}^{4}\cos ^{2}(\unicode[STIX]{x1D703})\\ c_{gz}^{2}\unicode[STIX]{x1D707}^{4}\cos (\unicode[STIX]{x1D703})\sin (\unicode[STIX]{x1D703})\\ -c_{gz}^{2}\unicode[STIX]{x1D707}^{5}|c_{gz}|\cos (\unicode[STIX]{x1D703})/N\end{array}\right]\nonumber\\ \displaystyle & & \displaystyle \times \,\text{e}^{\text{i}[((|c_{gz}|\unicode[STIX]{x1D707}^{2}\tilde{r}\cos (\unicode[STIX]{x1D703}-\unicode[STIX]{x1D703}_{p}))/N)+\unicode[STIX]{x1D707}\tilde{z}]}\,\text{sgn}(\unicode[STIX]{x1D707})\,\text{d}\unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D707},\end{eqnarray}$$
in which
$\unicode[STIX]{x1D705}_{H}=|c_{gz}|\unicode[STIX]{x1D707}^{2}/N$
in
$\widehat{u_{\mathit{DF}}}$
. For
$c_{gz}>0$
, the condition for outward-propagating waves requires
$\unicode[STIX]{x1D707}\cos (\unicode[STIX]{x1D703}-\unicode[STIX]{x1D703}_{p})\geqslant 0$
, so that
$\unicode[STIX]{x1D703}_{p}-\unicode[STIX]{x03C0}/2\geqslant \unicode[STIX]{x1D703}\geqslant \unicode[STIX]{x1D703}_{p}+\unicode[STIX]{x03C0}/2$
for
$\unicode[STIX]{x1D707}\geqslant 0$
and
$\unicode[STIX]{x1D703}_{p}+\unicode[STIX]{x03C0}/2<\unicode[STIX]{x1D703}<\unicode[STIX]{x1D703}_{p}+3\unicode[STIX]{x03C0}/2$
for
$\unicode[STIX]{x1D707}<0$
(and vice versa for
$c_{gz}<0$
). Separately integrating over positive and negative
$\unicode[STIX]{x1D707}$
and combining the results gives
$$\begin{eqnarray}\displaystyle \boldsymbol{u}_{\mathit{LW}} & = & \displaystyle \frac{2\unicode[STIX]{x03C0}}{N^{2}}\int _{0}^{\infty }\!\!\int _{\unicode[STIX]{x1D703}_{p}-\unicode[STIX]{x03C0}/2}^{\unicode[STIX]{x1D703}_{p}+\unicode[STIX]{x03C0}/2}\widehat{u_{\mathit{DF}}}\left[\begin{array}{@{}c@{}}c_{gz}^{2}\unicode[STIX]{x1D707}^{4}\cos ^{2}(\unicode[STIX]{x1D703})\\ c_{gz}^{2}\unicode[STIX]{x1D707}^{4}\cos (\unicode[STIX]{x1D703})\sin (\unicode[STIX]{x1D703})\\ -c_{gz}^{2}|c_{gz}|\unicode[STIX]{x1D707}^{5}\cos (\unicode[STIX]{x1D703})/N\end{array}\right]\nonumber\\ \displaystyle & & \displaystyle \times \sin \left(\frac{|c_{gz}|\unicode[STIX]{x1D707}^{2}\tilde{r}\cos (\unicode[STIX]{x1D703}-\unicode[STIX]{x1D703}_{p})}{N}+\unicode[STIX]{x1D707}\tilde{z}\,\text{sgn}(c_{gz})\right)\text{d}\unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D707}.\end{eqnarray}$$
Having removed the singularity, the double integral in (3.24) is readily solved by standard numerical integration techniques.
Equation (3.24) corresponds to the solution of the differential equation (31) in TA07 derived for flat wavepackets when neglecting rotation by setting the Coriolis parameter to zero. In our consideration of a wide, but not necessarily flat, wavepacket, the total induced mean flow corresponds to the sum of the Bretherton flow
$\boldsymbol{u}_{\mathit{BF}}$
in (3.15) and the induced long waves
$\boldsymbol{u}_{\mathit{LW}}$
in (3.24). Which of the two flows is dominant in the horizontal plane depends on the horizontal and vertical aspect ratios
$R_{y}\equiv \unicode[STIX]{x1D70E}_{y}/\unicode[STIX]{x1D70E}_{x}=\unicode[STIX]{x1D700}_{x}/\unicode[STIX]{x1D700}_{y}$
and
$R_{z}\equiv \unicode[STIX]{x1D70E}_{z}/\unicode[STIX]{x1D70E}_{x}=\unicode[STIX]{x1D700}_{x}/\unicode[STIX]{x1D700}_{z}$
. This is examined below in the specific case of a Gaussian wavepacket.
3.5 Dominant induced flows for different wavepacket aspect ratios
Although the equations above apply to any quasi-monochromatic wavepacket, here we specifically examine the predictions for Gaussian wavepackets with amplitude envelope
$A=A_{0}\exp [-(\tilde{x}^{2}/\unicode[STIX]{x1D70E}_{x}^{2}+{\tilde{y}}^{2}/\unicode[STIX]{x1D70E}_{y}^{2}+\tilde{z}^{2}/\unicode[STIX]{x1D70E}_{z}^{2})/2]$
, in which
$A_{0}$
is constant. Without loss of generality, we will assume that
$A_{0}$
is real and positive,
$k_{x}>0$
and
$k_{z}<0$
, so that the vertical group velocity is positive. Our focus is to compare the largest values of the horizontal velocity associated with the Bretherton flow (3.15) with that associated with induced long waves (3.24): the former should be dominant for relatively round wavepackets (small
$\unicode[STIX]{x1D70E}_{y}$
); the latter dominant for relatively wide wavepackets (large
$\unicode[STIX]{x1D70E}_{y}$
).
From (3.10), we see that the largest magnitude of the divergent-flux-induced flow is
$\Vert u_{\mathit{DF}}\Vert =N|\boldsymbol{k}||A_{0}|^{2}/2$
. Using (3.13), the Fourier-transformed divergent-flux-induced flow is
$\widehat{u_{\mathit{DF}}}=[N|\boldsymbol{k}|\unicode[STIX]{x1D70E}_{x}\unicode[STIX]{x1D70E}_{y}\unicode[STIX]{x1D70E}_{z}/(16\unicode[STIX]{x03C0}^{3/2})]|A_{0}|^{2}\exp [-(\unicode[STIX]{x1D705}^{2}\unicode[STIX]{x1D70E}_{x}^{2}+\unicode[STIX]{x1D706}^{2}\unicode[STIX]{x1D70E}_{y}^{2}+\unicode[STIX]{x1D707}^{2}\unicode[STIX]{x1D70E}_{z}^{2})/4]$
. We substitute this into (3.15) and extract the largest value of the horizontal induced circulation of the Bretherton flow, which occurs at the centre of the wavepacket. Thus, we find the scaled maximum flow,
$U_{\mathit{BF}}\equiv u_{\mathit{BF}}(0,0,0)/\Vert u_{\mathit{DF}}\Vert$
, to be (cf. § 14.3.3 of Bühler Reference Bühler2009)
$$\begin{eqnarray}U_{\mathit{BF}}=\frac{1}{4\unicode[STIX]{x03C0}}\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }\frac{\hat{\unicode[STIX]{x1D706}}^{2}}{R_{y}^{2}\hat{\unicode[STIX]{x1D705}}^{2}+\hat{\unicode[STIX]{x1D706}}^{2}}\exp (-\hat{\unicode[STIX]{x1D705}}^{2}/4-\hat{\unicode[STIX]{x1D706}}^{2}/4)\,\text{d}\hat{\unicode[STIX]{x1D705}}\,\text{d}\hat{\unicode[STIX]{x1D706}}=\frac{1}{R_{y}+1},\end{eqnarray}$$
in which
$\hat{\unicode[STIX]{x1D705}}\equiv \unicode[STIX]{x1D705}\unicode[STIX]{x1D70E}_{x}$
,
$\hat{\unicode[STIX]{x1D706}}\equiv \unicode[STIX]{x1D706}\unicode[STIX]{x1D70E}_{y}$
and
$R_{y}\equiv \unicode[STIX]{x1D70E}_{y}/\unicode[STIX]{x1D70E}_{x}$
denotes the horizontal aspect ratio of the packet. This result is plotted as the thick solid line in figure 1(a).
Figure 1. (a) Relative maximum value of the horizontal velocity and (b) its vertical location predicted as a function of the horizontal aspect ratio
$R_{y}$
of the wavepacket for the Bretherton flow (solid line), induced long waves (blue dashed lines) and their combination (thick red dashed line). For values of
$R_{y}$
for which our long-wave solution is not strictly valid, the blue dashed lines are replaced by blue dotted lines. Values are computed for a wavepacket with
$k_{z}/k_{x}=-1$
and
$\unicode[STIX]{x1D700}_{x}=\unicode[STIX]{x1D700}_{z}=1/5,1/20,1/100$
, as indicated (
$R_{z}=1$
). The circles indicate the results measured from numerical simulations with
$\unicode[STIX]{x1D700}_{x}=\unicode[STIX]{x1D700}_{z}=1/20$
and
$R_{y}=1$
, 10, 40 and 100. The black diamond indicates the solution found for the numerical simulation of a 2D wavepacket by van den Bremer & Sutherland (Reference van den Bremer and Sutherland2014).
The induced long waves, on the contrary, have zero horizontal flow at the centre of the wavepacket, with the maximum in the streamwise direction occurring vertically above the centre. Defining the non-dimensional horizontal flow along the
$\tilde{z}$
axis to be
$U_{\mathit{LW}}(\hat{z})\equiv u_{\mathit{LW}}(0,0,\hat{z})/\Vert u_{\mathit{DF}}\Vert$
, we have, from (3.24),
$$\begin{eqnarray}U_{\mathit{LW}}(\hat{z})=\unicode[STIX]{x1D700}_{x}^{2}\frac{R_{y}}{R_{z}^{4}}\frac{K^{2}}{4\sqrt{\unicode[STIX]{x03C0}}}\underbrace{\int _{0}^{\infty }\!\!\!\int _{-\unicode[STIX]{x03C0}/2}^{\unicode[STIX]{x03C0}/2}\!\exp \!\left[-\frac{\hat{\unicode[STIX]{x1D707}}^{2}}{4}-\frac{\unicode[STIX]{x1D700}_{x}^{2}K^{2}R_{y}^{2}\hat{\unicode[STIX]{x1D707}}^{4}\sin ^{2}(\unicode[STIX]{x1D703})}{4R_{z}^{4}}\right]\!\hat{\unicode[STIX]{x1D707}}^{4}\cos ^{2}(\unicode[STIX]{x1D703})\sin (\hat{\unicode[STIX]{x1D707}}\hat{z})\,\text{d}\unicode[STIX]{x1D703}\,\text{d}\hat{\unicode[STIX]{x1D707}}}_{\equiv \,I(\hat{z},\unicode[STIX]{x1D6E5}\equiv \unicode[STIX]{x1D700}_{x}KR_{y}/R_{z}^{2})},\end{eqnarray}$$
in which
$\hat{z}=\tilde{z}/\unicode[STIX]{x1D70E}_{z}$
and
$K\equiv k_{x}^{2}|k_{z}|/|\boldsymbol{k}|^{3}$
. Consistent with the assumption that the generated waves are long relative to the
$x$
extent of the packet, we can ignore the contribution
$\unicode[STIX]{x1D705}^{2}\unicode[STIX]{x1D70E}_{x}^{2}/4=\unicode[STIX]{x1D705}_{H}^{2}\cos ^{2}(\unicode[STIX]{x1D703})\unicode[STIX]{x1D70E}_{x}^{2}/4$
in the envelope function
$\widehat{u_{\mathit{DF}}}$
. The non-dimensional group
$\unicode[STIX]{x1D6E5}$
plays the role of an effective wavepacket aspect ratio. For large
$R_{y}$
, (3.26) reduces to its 2D equivalent (cf. (52) of van den Bremer & Sutherland Reference van den Bremer and Sutherland2014). For small
$R_{y}$
, (3.26) corresponds to the wave field radiated by a point source in a 3D domain, noting the restriction that the packet remains wide (cf. § 3.4). How large
$R_{y}$
has to be depends on the bandwidth parameter
$\unicode[STIX]{x1D700}_{x}$
, the orientation of the wavevector
$K$
and the aspect ratio
$R_{z}$
, as captured by
$\unicode[STIX]{x1D6E5}$
.
We begin by considering the case for which
$R_{z}=1$
: the wavepacket is round in the
$x$
–
$z$
plane, but of arbitrary width in the spanwise direction. The maximum and its relative vertical location are computed numerically and are plotted as a function of
$R_{y}$
as the blue dashed lines in figure 1 for three different values of
$\unicode[STIX]{x1D700}_{x}$
, as indicated; in both panels, we set
$k_{z}=-k_{x}$
. It is worth emphasizing here that the long-wave contribution (3.26) is derived based on the assumption of a wide wavepacket:
$\unicode[STIX]{x1D700}_{y}\sim \unicode[STIX]{x1D700}_{z}^{2}$
. Taking a practical approach, we thus set
$\unicode[STIX]{x1D700}_{y}\leqslant \unicode[STIX]{x1D700}_{z}^{2}$
as its domain of validity, which corresponds to
$R_{y}\geqslant R_{z}^{2}/\unicode[STIX]{x1D700}_{x}$
. For the values
$\unicode[STIX]{x1D700}_{x}=\unicode[STIX]{x1D700}_{z}=1/5,1/20,1/100$
considered in figure 1, the dashed blue lines are thus only strictly valid for values of
$R_{y}$
greater than 5, 20 and 100, as indicated by the transition of the blue dashed lines into blue dotted lines. It is evident that, for smaller values of
$R_{y}$
, the magnitude of the long waves is sufficiently small compared to the Bretherton flow, so the long-wave contribution can be ignored regardless.
Figure 2. Maximum value of the double integral
$I(\hat{z},\unicode[STIX]{x1D6E5})$
in (3.26) as a function of
$\unicode[STIX]{x1D6E5}\equiv \unicode[STIX]{x1D700}_{x}KR_{y}/R_{z}^{2}$
. Asymptotic limits are indicated by the dashed lines: for
$\unicode[STIX]{x1D6E5}\ll 0.1$
,
$I^{\star }\sim 29.71$
; for
$\unicode[STIX]{x1D6E5}\gg 1$
,
$I^{\star }\sim 10.30\unicode[STIX]{x1D6E5}^{-1}$
.
More generally, the asymptotic behaviour of
$U_{\mathit{LW}}$
is determined first by finding the maximum, as it depends upon
$\tilde{z}$
, of the double integral
$I(\hat{z},\unicode[STIX]{x1D6E5})$
in (3.26), as a function of
$\unicode[STIX]{x1D6E5}\equiv \unicode[STIX]{x1D700}_{x}KR_{y}/R_{z}^{2}$
. The resulting value, denoted by
$I^{\star }(\unicode[STIX]{x1D6E5})$
, is plotted in figure 2. The double integral asymptotically approaches a constant value for small
$\unicode[STIX]{x1D6E5}$
, and varies as
$\unicode[STIX]{x1D6E5}^{-1}$
for large
$\unicode[STIX]{x1D6E5}$
. The non-dimensional height where the maximum occurs varies only little, increasing from
$\tilde{z}\simeq 0.474$
for small
$\unicode[STIX]{x1D6E5}$
to
$\tilde{z}\simeq 0.596$
for large
$\unicode[STIX]{x1D6E5}$
. Taking a practical approach, the crossover from the small- to large-
$\unicode[STIX]{x1D6E5}$
asymptotic regimes is estimated from the intersection of the two asymptotic curves, occurring for
$\unicode[STIX]{x1D6E5}=\unicode[STIX]{x1D6E5}_{c}\simeq 0.35$
. Combining these results with the terms in front of the double integral in (3.26), we arrive at the following asymptotic approximations for the maximum induced flow in the
$x$
-direction associated with induced long waves:
$$\begin{eqnarray}\max _{\tilde{z}}\{U_{\mathit{LW}}\}\simeq \left\{\begin{array}{@{}ll@{}}4.19\,K^{2}\unicode[STIX]{x1D700}_{x}^{2}{\displaystyle \frac{R_{y}}{R_{z}^{4}}} & \text{for}~\unicode[STIX]{x1D6E5}\ll \unicode[STIX]{x1D6E5}_{c},\\ 1.45\,K\unicode[STIX]{x1D700}_{x}{\displaystyle \frac{1}{R_{z}^{2}}} & \text{for}~\unicode[STIX]{x1D6E5}\gg \unicode[STIX]{x1D6E5}_{c}.\end{array}\right.\end{eqnarray}$$
By comparing the maximum horizontal velocity of the Bretherton flow (3.25) and the limits of and the maximum horizontal velocity associated with induced long waves (3.27), we construct the regime diagram, shown in figure 3, which illustrates whether the flow induced by the wavepacket over its extent is best described by that of 1D, 2D or fully 3D wavepackets, as it depends upon the horizontal and vertical aspect ratios
$R_{y}$
and
$R_{z}$
. We will refer to these induced flows as quasi-1D, quasi-2D and quasi-3D and summarize the regimes in table 2. In order to draw boundaries to delineate what is a smooth transition between the different regimes, we make the following assumptions.
Figure 3. Regime diagram indicating for what vertical and horizontal aspect ratios (
$R_{z}$
and
$R_{y}$
, respectively) a Gaussian wavepacket will induce a flow similar to that of a 1D wavepacket (medium grey), a 2D wavepacket (light grey), which induces long waves, and a fully 3D wavepacket (white), which induces a horizontal circulation known as the Bretherton flow. The long-wave transition boundary is computed using (3.28) taking
$|k_{z}|=|k_{x}|$
(so that
$K=2^{-3/2}$
) and
$\unicode[STIX]{x1D700}_{x}=1/20$
. The cross-hatched regions indicate where the quasi-monochromatic wavepacket assumption breaks down for
$\unicode[STIX]{x1D700}_{x}=1/20$
because either
$\unicode[STIX]{x1D700}_{y}>1/5$
or
$\unicode[STIX]{x1D700}_{z}>1/5$
. The region to the right of the dotted line corresponds to the domain of validity of the long-wave solution. The closed circles indicate the values of
$R_{y}$
and
$R_{z}$
for the five simulations with snapshots shown in figure 5.
A wavepacket is considered to be quasi-1D if the wave-induced flow has a magnitude within 10 % of the divergent-flux-induced flow
$u_{\mathit{DF}}$
. From (3.25), this occurs for
$R_{y}\lesssim 0.1$
. This transitional value is plotted as the vertical long-dashed line in figure 3.
A wavepacket is considered to be quasi-2D if the flow induced by long waves exceeds the Bretherton flow. An estimate of the critical vertical aspect ratio of the wavepacket
$R_{z}^{\star }$
(as a function of
$R_{y}$
) when these flows are comparable is given by equating
$U_{\mathit{BF}}$
in (3.25) to
$U_{\mathit{LW}}$
in (3.27), in which the asymptotic approximations are extrapolated outside their strict domains of validity to
$\unicode[STIX]{x1D6E5}=\unicode[STIX]{x1D6E5}_{c}$
. For
$R_{y}\gtrsim 1.0$
, the large-
$\unicode[STIX]{x1D6E5}$
limit dominates in (3.27). Otherwise the small-
$\unicode[STIX]{x1D6E5}$
limit dominates. Together we find the transition boundary is given by
$$\begin{eqnarray}R_{z}^{\star }\simeq \left\{\begin{array}{@{}ll@{}}1.43[K^{2}\unicode[STIX]{x1D700}_{x}^{2}R_{y}(R_{y}+1)]^{1/4}\sim 1.43[K\unicode[STIX]{x1D700}_{x}]^{1/2}R_{y}^{1/4} & \text{for}~R_{y}\lesssim 1.0,\\ 1.20[K\unicode[STIX]{x1D700}_{x}(R_{y}+1)]^{1/2}\sim 1.20[K\unicode[STIX]{x1D700}_{x}]^{1/2}R_{y}^{1/2} & \text{for}~R_{y}\gtrsim 1.0,\end{array}\right.\end{eqnarray}$$
where the second similarity in two respective limits follows from taking the small- and large-
$R_{y}$
limits of the
$(1+R_{y})$
term. Equation (3.28) is plotted as the short-dashed curve in figure 3 for the case of a wavepacket with
$|k_{z}|=|k_{x}|$
(
$K=2^{-3/2}$
) and
$\unicode[STIX]{x1D700}_{x}=1/20$
. In particular, for large
$R_{y}$
and this choice of parameters, the transition curve between quasi-2D and quasi-3D wavepackets is given approximately by
$R_{z}^{\star }\simeq 0.161R_{y}^{1/2}$
. Conversely, this transition written in terms of the horizontal aspect ratio as it depends upon the vertical aspect ratio is
$R_{y}^{\star }\simeq [0.689/(K\unicode[STIX]{x1D700}_{x})]R_{z}^{2}\simeq 39.0R_{z}^{2}$
.
Evidently, it is worth emphasizing again that the long-wave contribution (3.26) is derived based on the assumption of a wide wavepacket and that its validity is thus limited to
$R_{z}\leqslant \sqrt{\unicode[STIX]{x1D700}_{x}R_{y}}$
, denoted by the dotted line in figure 3. By equating the value of
$R_{z}$
corresponding to this boundary to the transition value
$R_{z}^{\star }$
in (3.28), we can obtain an estimate of how wide the wavepacket needs to be for our classification to be valid:
$$\begin{eqnarray}R_{y,min}\simeq \left\{\begin{array}{@{}ll@{}}1.43^{4}K^{2}/(1-1.43^{4}K^{2}) & \text{for}~R_{y}\lesssim 1.0,\\ 1.20^{2}K/(1-1.20^{2}K) & \text{for}~R_{y}\gtrsim 1.0.\end{array}\right.\end{eqnarray}$$
This is plotted in figure 4 as a function of the angle
$\unicode[STIX]{x1D703}=\tan ^{-1}(|k_{z}/k_{x}|)$
that lines of constant phase make with the vertical. This minimum width is not dependent on
$\unicode[STIX]{x1D700}_{x}$
. For the case
$|k_{z}|=|k_{x}|$
considered herein,
$\unicode[STIX]{x1D703}=45^{\circ }$
and
$R_{y,min}=1.10$
. Our estimates in figure 3 are thus only valid for wavepackets that are at least marginally wider than they are long in the
$x$
-direction, i.e. for
$R_{y}\geqslant R_{y,min}=1.10$
.
Figure 4. Minimum width for which the regime diagram is valid as a function of the angle
$\unicode[STIX]{x1D703}$
that lines of constant phase make with the vertical (
$k_{x}=|\boldsymbol{k}|\cos (\unicode[STIX]{x1D703})$
).
Finally, it is also important to keep in mind that these results rely on the starting assumption that the wavepacket is quasi-monochromatic, so that dispersion is a negligible effect. Hence all of
$\unicode[STIX]{x1D700}_{x}$
,
$\unicode[STIX]{x1D700}_{y}$
and
$\unicode[STIX]{x1D700}_{z}$
must be much smaller than unity. For given
$\unicode[STIX]{x1D700}_{x}$
, this puts a limitation on the aspect ratios
$R_{y}=\unicode[STIX]{x1D70E}_{y}/\unicode[STIX]{x1D70E}_{x}=\unicode[STIX]{x1D700}_{x}/\unicode[STIX]{x1D700}_{y}$
and
$R_{z}=\unicode[STIX]{x1D70E}_{z}/\unicode[STIX]{x1D70E}_{x}=\unicode[STIX]{x1D700}_{x}/\unicode[STIX]{x1D700}_{z}$
for which the transition boundaries can be considered reliable. In particular, the cross-hatched regions in figure 3 show the values of
$R_{y}<5\unicode[STIX]{x1D700}_{x}$
and
$R_{z}<5\unicode[STIX]{x1D700}_{x}$
for which the wavepacket is, respectively, too narrow (
$\unicode[STIX]{x1D700}_{y}>1/5$
) or insufficiently tall (
$\unicode[STIX]{x1D700}_{z}>1/5$
) for the case
$\unicode[STIX]{x1D700}_{x}=1/20$
. The quasi-1D wavepacket regime and the long-wave regime in the small
$R_{y}$
limit lie entirely within the cross-hatched regions for the value of
$\unicode[STIX]{x1D700}_{x}$
considered here.
The different asymptotic regimes are summarized in table 2. In theory one could examine cases in which
$\unicode[STIX]{x1D700}_{x}$
is so small that transitions to flows associated with 1D wavepackets and 2D wavepackets at small
$R_{y}$
could be reliably predicted. However, in the numerical simulations that follow, we restrict our test of theory to the examination of wavepackets with
$\unicode[STIX]{x1D700}_{x}=1/20$
. This is because simulations with significantly smaller
$\unicode[STIX]{x1D700}_{x}$
would require much larger domain sizes and consequently more memory and longer computation time. Furthermore, vertically propagating wavepackets with
$\unicode[STIX]{x1D70E}_{x}\gg 20k_{x}^{-1}\simeq 3\unicode[STIX]{x1D706}_{x}$
are unlikely to be generated by a realistic geophysical system on mesoscales and larger.
4 Comparison with numerical simulations
4.1 Code description
The flows induced by a Gaussian wavepacket were examined using fully nonlinear numerical simulations that solved the equations (2.1) with the addition of damping terms to the momentum and internal energy equations for the purpose of numerical stability. Specifically, under the assumption that
$N^{2}$
is constant, the code evolved in time the horizontal velocity and vertical displacement fields according to
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x2202}_{t}u=-\unicode[STIX]{x2202}_{x}\text{(}uu\text{)}-\unicode[STIX]{x2202}_{y}\text{(}vu\text{)}-\unicode[STIX]{x2202}_{z}\text{(}wu\text{)}-\unicode[STIX]{x2202}_{x}P+\unicode[STIX]{x1D708}{\mathcal{D}}\text{(}u\text{)}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x2202}_{t}v=-\unicode[STIX]{x2202}_{x}\text{(}uv\text{)}-\unicode[STIX]{x2202}_{y}\text{(}vv\text{)}-\unicode[STIX]{x2202}_{z}\text{(}wv\text{)}-\unicode[STIX]{x2202}_{y}P+\unicode[STIX]{x1D708}{\mathcal{D}}\text{(}v\text{)}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D709}=-\unicode[STIX]{x2202}_{x}(u\unicode[STIX]{x1D709})-\unicode[STIX]{x2202}_{y}(v\unicode[STIX]{x1D709})-\unicode[STIX]{x2202}_{z}(w\unicode[STIX]{x1D709})+w+\unicode[STIX]{x1D705}{\mathcal{D}}(\unicode[STIX]{x1D709}), & \displaystyle\end{eqnarray}$$
$w$
and pressure divided by density
$P\equiv p/\unicode[STIX]{x1D70C}_{0}$
were found from the respective diagnostic equations: 
$$\begin{eqnarray}\unicode[STIX]{x2202}_{z}w=-\unicode[STIX]{x2202}_{x}u-\unicode[STIX]{x2202}_{y}v\end{eqnarray}$$
and
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FB}^{2}P & = & \displaystyle -\unicode[STIX]{x2202}_{xx}\text{(}uu\text{)}-\unicode[STIX]{x2202}_{yx}\text{(}vu\text{)}-\unicode[STIX]{x2202}_{zx}\text{(}wu\text{)}\nonumber\\ \displaystyle & & \displaystyle -\,\unicode[STIX]{x2202}_{xy}\text{(}uv\text{)}-\unicode[STIX]{x2202}_{yy}\text{(}vv\text{)}-\unicode[STIX]{x2202}_{zy}\text{(}wv\text{)}\nonumber\\ \displaystyle & & \displaystyle -\,\unicode[STIX]{x2202}_{xz}\text{(}uw\text{)}-\unicode[STIX]{x2202}_{yz}\text{(}vw\text{)}-\unicode[STIX]{x2202}_{zz}\text{(}ww\text{)}-N^{2}\unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D709}.\end{eqnarray}$$
The equations were solved in Fourier space (corresponding to
$x$
,
$y$
and
$z$
) in a triply periodic domain. The nonlinear (e.g. advective) terms in the equations were computed by fast Fourier transforming the fields to real space, multiplying and then transforming back to Fourier space. The diffusion operator
${\mathcal{D}}$
was prescribed to act as a Laplacian operator only upon horizontal wavenumbers greater than
$2k_{x}$
. The diffusivities of momentum and substance acting on these high wavenumbers were set to be
$\unicode[STIX]{x1D708}=\unicode[STIX]{x1D705}=0.001N/k_{x}^{2}$
. The fields were advanced in time using a leapfrog scheme with an Euler backstep taken every 20 steps.
4.2 Initial conditions
The initial wavepacket centred at the origin had wavenumber set so that
$k_{z}=-k_{x}$
, and the Gaussian envelope had amplitude
$A_{0}=0.01/k_{x}$
(
$\unicode[STIX]{x1D6FC}=0.01$
) and
$x$
extent
$\unicode[STIX]{x1D70E}_{x}=20/k_{x}$
(
$\unicode[STIX]{x1D700}_{x}=1/20$
). The details of the five simulations presented here are summarized in table 3. We first examined four vertically round (
$\unicode[STIX]{x1D70E}_{z}=\unicode[STIX]{x1D70E}_{x}$
,
$R_{z}=1$
) wavepackets with lateral extents
$k_{x}\unicode[STIX]{x1D70E}_{y}=20$
, 200, 800 and 2000, thus exploring the range between round and wide wavepackets with
$R_{y}=1$
, 10, 40 and 100, bracketing the predicted critical transition value of
$R_{y}^{\star }=39.0$
, as derived below (3.28). Furthermore, we examined a more flat wavepacket chosen to lie on our transition curve in figure 3 (
$R_{y}=10$
,
$R_{z}=0.5$
).
Table 3. Comparison of the scaled maximum horizontal velocity
$U^{(2)}$
from numerical simulations and perturbation theory.
The choice of initial conditions for the mean flow is crucial for the interpretation of results from numerical simulations. Up to the orders we consider (see also § 5), our analytical solution consists of a packet that translates at the group velocity, and the mean flow (Bretherton flow and long waves) it induces is entirely steady. Numerical solutions intended to reproduce this steady solution should therefore in principle have both aspects, the linear packet and the second-order mean flow (Bretherton flow and long waves), superimposed at the outset. If the second-order solutions are not superimposed at the outset, the flow responds as follows. The second-order (amplitude) equations instantaneously demand the generation of a mean flow (Bretherton flow and long waves) that are bound to the packet and travel with it. An ‘error wave’ (or flow) then forms that is equal and of opposite sign to the mean flow, so that the total second-order signal is zero, as prescribed by the initial conditions and by conservation of momentum. The error wave does not generally travel at the same speed as the packet. Here, it stays behind. The error wave can be thought of as unphysical in the sense that it is not part of a (quasi-)steady solution. If one is interested in this (quasi-)steady solution, it is imperative to wait long enough for the solution and the error wave to separate (cf. van den Bremer & Sutherland Reference van den Bremer and Sutherland2014).
Here, we take the following practical approach. Out of the two parts of the second-order mean flow (the Bretherton flow and the long waves), we have only initialized with the Bretherton flow predicted by (3.15). This way, we do not leave an error flow (with opposite sign to the Bretherton flow) behind and can see whether waves are indeed generated, but we do leave behind error waves (with opposite sign to the long waves). We are left with the subtlety that the long waves will consist of both the actual long waves travelling with the packet and the error wave left behind. We will have to wait for long enough for the two to separate, conceptually analogous to what was done by van den Bremer & Sutherland (Reference van den Bremer and Sutherland2014).
4.3 Domain size and post-processing
In all cases, the domain was set to be much larger than the anticipated spatial extent of the induced flow over the duration of the simulation, up to time
$t=200/N$
when the wavepacket had propagated sufficiently far away from its initial position centred at the origin. The simulation resolved disturbances with wavenumbers up to
$4k_{x}$
and
$4|k_{z}|$
in the
$x$
- and
$z$
-directions. In the case
$R_{y}=1$
, the spanwise domain range was
$|y|\leqslant L_{y}$
with
$L_{y}=100k_{x}^{-1}$
(
$k_{z}=-k_{x}$
) and resolved by wavenumbers up to
$128\unicode[STIX]{x03C0}/L_{y}$
. In the wider cases, the domain was substantially larger to accommodate the generation of long waves and so the resolution was reduced in order for the numerical simulations to complete in a reasonable time frame (see table 3), with wavenumbers being resolved up to
$32\unicode[STIX]{x03C0}/L_{y}$
. The time resolution was
$0.010/N$
in the case
$R_{y}=1$
and
$0.025/N$
in the cases
$R_{y}=10$
, 40 and 100. Each run took four days to complete running in serial on a 2.9 GHz Intel Core i5.
The flows induced by the wavepacket were assessed by filtering the Fourier modes with horizontal wavenumbers between
$0.5k_{x}$
and
$2k_{x}$
corresponding to the linear waves that make up the wavepacket. The inverse transform thereby revealed the flow induced by the propagating wavepacket. The result of this analysis at the end of five simulations is shown in figure 5, which plots cross-sections of the
$x$
-component of the induced flow at time
$t=200/N$
, when the wavepacket had propagated sufficiently far away from its initial position at the origin, being centred at
$(x,y,z)\simeq (70.7,0,70.7)$
. (A multimedia version of figure 5 is available online at https://doi.org/10.1017/jfm.2017.745, which shows 3D representations of the induced flows and their cross-sections.)
Figure 5. The wavepacket-filtered
$x$
-component of velocity determined at time
$Nt=200$
from five simulations initialized with a Gaussian wavepacket of maximum amplitude
$A_{0}=0.01k_{x}^{-1}$
centred at the origin, having vertical wavenumber
$k_{z}=-k_{x}$
,
$x$
extents of
$\unicode[STIX]{x1D70E}_{x}=20k_{x}^{-1}$
, and
$y$
and
$z$
extents as reported in the panels. The left (right) panels show a cross-section of the flow in the
$y$
–
$z$
(
$x$
–
$z$
) plane through the centre of the wavepacket. The thick black ovals are drawn around one standard deviation from the centre of the wavepacket predicted to be at
$(c_{gx}t,0,c_{gz}t)$
; these ovals appear as vertical lines in the right panels of (b)–(e) because the extents in
$x$
are small compared to the width of the domain. The colour bars to the lower left in the left panels indicate the magnitude of the induced flows in the left and corresponding right panels.
4.4 Results
Consistent with predictions, the round wavepacket with
$\unicode[STIX]{x1D70E}_{y}=\unicode[STIX]{x1D70E}_{z}=\unicode[STIX]{x1D70E}_{x}$
(
$R_{y}=R_{z}=1$
) induces a horizontal circulation only (the Bretherton flow), as shown in figure 5(a). There is no trace of the Bretherton flow that was initially superimposed on the wavepacket at the origin, as it has translated with the wavepacket at the predicted group velocity. The maximum
$x$
-component of the induced flow at this time indeed occurs at the centre of the wavepacket, with a value consistent with that predicted by theory (see figure 1 and table 3). In simulations with a spanwise wider wavepacket (
$R_{y}$
large), keeping
$R_{z}=1$
fixed, the magnitude of the Bretherton flow becomes smaller and the flow associated with long waves becomes larger. For
$R_{y}=10$
(
$\unicode[STIX]{x1D70E}_{y}=200/k_{x}$
), the largest flow associated with long waves is predicted to be 6 % that of the maximum Bretherton flow. Indeed, a weak wave signature is seen to be superimposed on the Bretherton flow in the snapshot cross-sections taken from the numerical simulation shown in figure 5(b).
Also as predicted, long waves eventually dominate in the simulation with
$R_{y}=100>R_{y}^{\star }$
(
$\unicode[STIX]{x1D70E}_{y}=2000/k_{x}$
), as shown in figure 5(d). In this case the
$x$
-component of the induced flow is close to zero at the centre of the wavepacket, with the maximum shifted up relative to the centre of the packet, as predicted by theory and shown in figure 1. The small discrepancies between the simulated and predicted values of
$\max (U^{(2)})$
and its vertical location, as compared in figure 1(b), are attributed to the necessarily coarse spatial and temporal resolution of the simulations.
In the two simulations for which the Bretherton flow and long waves are predicted to have comparable amplitudes (
$R_{y}=40$
and
$R_{z}=1$
, shown in figure 5
c;
$R_{y}=10$
and
$R_{z}=0.5$
, shown in figure 5
e) at
$t=200/N$
, the superposition is evident of the Bretherton flow, for which a negative horizontal induced flow exists to either side of the wavepacket, and long waves, for which horizontal velocity perturbations occur in the lee and to either side of the wavepacket. The errors associated with our predictions are now somewhat larger, but remain small, as is evident from table 3.
5 Comparison with Tabaei & Akylas (Reference Tabaei and Akylas2007)
Tabaei & Akylas (Reference Tabaei and Akylas2007, TA07) consider two cases of wavepackets having finite extent in all three spatial dimensions: ‘three-dimensional equally strong modulations’ (i.e. round wavepackets for which
$R_{y}\sim R_{z}\sim 1$
) and ‘flat wavetrains’ (i.e. wavepackets for which
$R_{y}\sim 1$
and
$R_{z}\ll 1$
). They also consider the effects of rotation, finite-amplitude effects and they revisit spanwise-infinite wavepackets (‘two-dimensional equally strong modulations’, i.e.
$R_{y}\rightarrow \infty$
and
$R_{z}\sim 1$
).
A significant contribution of TA07 is the examination of modulational stability of a wavetrain through interaction with the wave-induced flow. Specifically, TA07 examine the stability of weakly nonlinear flat wavetrains (their (66)–(68)) and show that perturbations can always be unstable (cf. their (79)), except when the wavepacket remains periodic in
$x$
and there is zero rotation. In the latter case, TA07 recover the stability criterion derived by Sutherland (Reference Sutherland2006b
) for 1D wavepackets (see § 1). More generally, TA07 demonstrate that the instability is strongest when the two sets of roots from their stability analysis coalesce, these representing free hydrostatic gravity-inertial waves (their (77)) and modulations of the horizontal wavetrain travelling with the vertical group velocity (their (78)). Physically, this corresponds to the self-interaction of the packet with the 3D long waves that it generates. For the weakly nonlinear equations in TA07 and to leading order, the equations for the induced mean flow can be solved analytically for a given amplitude envelope by the solutions we propose herein. This explicitly computed mean flow affects the way the wavepacket evolves through the nonlinear term in the nonlinear Schrödinger-type equations in TA07, effectively through Doppler shifting the waves within the wavepacket. If the amplitude of these waves is sufficiently large, the Doppler shifting by the induced mean flow changes the extrinsic frequency of the wave and hence the structure of the amplitude envelope, as prescribed by the advective and dispersive terms in the nonlinear Schrödinger-type equations in TA07. The change in amplitude envelope then requires recalculation of the corresponding induced mean flow at every time given the shape of the packet at that time.
In the work presented here we have assumed that the amplitude of the waves is sufficiently small that the weakly nonlinear feedback between the induced mean flow and the waves is negligible. We also ignore changes of shape of the packet due to linear dispersion. Thus, the induced mean flow is steady in a frame moving with the constant group velocity of the wavepacket. Nonetheless, we are able to compare our predictions and the results of our numerical simulations with those of TA07, who presented simulation results for relatively small-amplitude waves in the particular limits of round wavepackets (§ 5.1) and flat wavetrains (§ 5.2).
5.1 Round wavepackets
If the modulations are equally strong in all three directions and indeed as strong as the effect of rotation, TA07 also conclude that the velocity is essentially horizontal (their § 2.1). Without making restrictive assumptions about the amplitude of the waves, the mean flow is then governed by their (17) and (18), which is in agreement with Grimshaw (Reference Grimshaw1972) in the absence of rotation, as noted by these authors. Without rotation and in the weakly nonlinear limit, their (17) and (18) reduced to their (22), which in turn agrees with Shrira (Reference Shrira1981) (and Bretherton Reference Bretherton1969), as noted by TA07. In the small-amplitude limit, our prediction given by (3.14) is the solution to (22) in TA07.
The right panel of our figure 5(a) can be compared at least qualitatively with the induced flow shown by line contours in the left panel of figure 4 of TA07. Their simulation was somewhat different: it included the effects of rotation, and the vertical wavenumber was relatively smaller in magnitude. Despite these differences, two instructive comparisons can be made. First, because TA07 did not initialize the simulation with the predicted Bretherton flow, they observe a negative, non-translating, induced flow at the origin in addition to the positive induced flow that translates with the wavepacket. Second, because their wavepacket was round in this simulation, TA07 also did not observe the generation of long waves. Our figure 1 indeed confirms that the maximum induced flow associated with long waves is three orders of magnitude smaller than the Bretherton flow in the case
$R_{y}=1$
and
$\unicode[STIX]{x1D700}_{x}=1/20$
(and
$R_{z}=1$
).
5.2 Flat wavetrains
Assuming the horizontal extent of the wavepacket is comparable in the
$x$
- and
$y$
-directions but both are much larger than the vertical extent, TA07 predicted the wavepacket should excite long waves according to their (31). They solved this prognostic differential equation in a simulation with no rotation, setting
$|k_{z}/k_{x}|=0.4$
(so
$K=0.32$
), and setting the horizontal wavepacket extent in both directions to be 100 times larger than the vertical extent indeed shows trailing waves (so
$R_{y}=1$
,
$R_{z}=0.01$
and
$\unicode[STIX]{x1D700}_{x}=10^{-4}$
). This indeed showed the generation of trailing long waves, as revealed in figure 2 of TA07. Yet evidence of an induced flow at the origin is indicative of an error flow resulting from not superimposing the predicted induced flow at the outset and thus considering an inherently unsteady problem.
Generally, we predict that a flat wavepacket should excite long waves of comparable or larger amplitude than the Bretherton flow if
$R_{z}\lesssim R_{z}^{\star }$
, with
$R_{z}^{\star }$
given by (3.28). With
$R_{y}=1$
,
$K=0.32$
and
$\unicode[STIX]{x1D700}_{x}=10^{-4}$
, we predict the transition occurs at
$R_{z}^{\star }\simeq 0.007$
. Thus, the simulation of TA07 shown in their figure 2 corresponds to the regime of induced flows dominated by the Bretherton flow, but which is nonetheless close to the transition in which induced long waves have comparable amplitude. This explains the appearance of the Bretherton-flow-like error flow at the origin as well as the manifestation of induced long waves trailing the propagating wavepacket.
6 Conclusions
Historically, the induction of the Bretherton flow and of long waves by internal wavepackets have been treated as two distinct problems, the former resulting from a round wavepacket and the latter resulting from a flat or infinitely wide wavepacket (Bretherton Reference Bretherton1969; Bühler & McIntyre Reference Bühler and McIntyre1998; Tabaei & Akylas Reference Tabaei and Akylas2007). Using perturbation theory for quasi-monochromatic internal wavepackets, we have rederived the equations in an approach that lends physical insight to the connection and overlap of these distinct flows: the Bretherton flow, which has zero vertical velocity and is forced by vertical vorticity, as has been derived by several authors, including Bühler & McIntyre (Reference Bühler and McIntyre1998); and the induced long waves, which have non-negligible vertical velocity and are forced by horizontal vorticity. We also provide explicit integral solutions for the 3D induced long waves. Combined, these predict the mean flow induced by wavepackets of arbitrary vertical and horizontal aspect ratio, provided they are marginally wider than round in the horizontal. In particular, these explicit integral solutions allow us to assess for which aspect ratios induced long waves dominate over the Bretherton flow and lead us to predict a transitional vertical wavepacket aspect ratio
$R_{z}^{\star }$
, as it depends upon the horizontal aspect ratio
$R_{y}$
through (3.28).
Numerical simulations examined the evolution of small-amplitude wavepackets with a range of aspect ratios designed to explore the transition from induced flows dominated by the Bretherton flow to induced flows dominated by long waves. In accordance with theory, the simulations predict well the peak amplitudes of the induced horizontal flow and their locations. In particular, the results clearly show that the Bretherton flow and induced long waves both contribute to the observed induced flow particularly near the transition boundary.
Acknowledgements
The authors are grateful to the three reviewers whose constructive comments helped greatly to improve this article. This research was performed in part through funding from the Natural Sciences and Engineering Research Council (NSERC) of Canada. T.S.v.d.B. was partially supported by the Faculty of Science at the University of Alberta through its Visiting Fellowship programme.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2017.745.
Appendix A. Explicit forcing equation
Here we explicitly give the three components of the matrix equation (2.3). From substitution of the linear (in
$\unicode[STIX]{x1D6FC}$
) fields in table 1 (both columns), including only those terms that are non-zero, we obtain
$$\begin{eqnarray}\displaystyle -\unicode[STIX]{x2202}_{ttz}v^{\text{(}2\text{)}}+(\unicode[STIX]{x2202}_{tt}+N^{2})\unicode[STIX]{x2202}_{y}\text{}w\text{}^{\text{(}2\text{)}} & = & \displaystyle -\unicode[STIX]{x2202}_{t}[\unicode[STIX]{x2202}_{x}(u_{0}\unicode[STIX]{x1D701}_{x1})+\unicode[STIX]{x2202}_{z}(w_{0}\unicode[STIX]{x1D701}_{x1})]\nonumber\\ \displaystyle & & \displaystyle +\,N^{2}\unicode[STIX]{x2202}_{y}[\unicode[STIX]{x2202}_{x}(u_{1}\unicode[STIX]{x1D709}_{0})+\unicode[STIX]{x2202}_{y}(v_{1}\unicode[STIX]{x1D709}_{0})+\unicode[STIX]{x2202}_{z}(w_{1}\unicode[STIX]{x1D709}_{0})]\nonumber\\ \displaystyle & & \displaystyle +\,\unicode[STIX]{x2202}_{t}[\unicode[STIX]{x1D701}_{x1}\unicode[STIX]{x2202}_{x}u_{0}+\unicode[STIX]{x1D701}_{y0}\unicode[STIX]{x2202}_{y}u_{1}+\unicode[STIX]{x1D701}_{y1}\unicode[STIX]{x2202}_{y}u_{0}],\qquad \quad\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{ttz}u^{\text{(}2\text{)}}-(\unicode[STIX]{x2202}_{tt}+N^{2})\unicode[STIX]{x2202}_{x}w^{\text{(}2\text{)}} & = & \displaystyle -\unicode[STIX]{x2202}_{t}[\unicode[STIX]{x2202}_{x}(u_{0}\unicode[STIX]{x1D701}_{y1}+u_{1}\unicode[STIX]{x1D701}_{y0})+\unicode[STIX]{x2202}_{y}(v_{1}\unicode[STIX]{x1D701}_{y0})+\unicode[STIX]{x2202}_{z}(w_{0}\unicode[STIX]{x1D701}_{y1}+w_{1}\unicode[STIX]{x1D701}_{y0})]\nonumber\\ \displaystyle & & \displaystyle -\,N^{2}\unicode[STIX]{x2202}_{x}[\unicode[STIX]{x2202}_{x}(u_{1}\unicode[STIX]{x1D709}_{0})+\unicode[STIX]{x2202}_{y}(v_{1}\unicode[STIX]{x1D709}_{0})+\unicode[STIX]{x2202}_{z}(w_{1}\unicode[STIX]{x1D709}_{0})]+\unicode[STIX]{x2202}_{t}[\unicode[STIX]{x1D701}_{y0}\unicode[STIX]{x2202}_{y}v_{1}],\qquad \quad\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle -\unicode[STIX]{x2202}_{tty}u^{\text{(}2\text{)}}+\unicode[STIX]{x2202}_{ttx}v^{\text{(}2\text{)}}=\unicode[STIX]{x2202}_{t}[\unicode[STIX]{x1D701}_{x1}\unicode[STIX]{x2202}_{x}w_{0}+\unicode[STIX]{x1D701}_{y0}\unicode[STIX]{x2202}_{y}w_{1}+\unicode[STIX]{x1D701}_{y1}\unicode[STIX]{x2202}_{y}w_{0}], & & \displaystyle\end{eqnarray}$$
where we have left out the superscripts (1) on the linear (in
$\unicode[STIX]{x1D6FC}$
) fields on the right-hand side for convenience. It is also understood that the expressions on the right-hand side are averaged over the fast spatial scale of waves in the wavepacket.
