2.1. Basic equations and physical quantities

In this work, we consider a Boussinesq incompressible flow with constant buoyancy frequency $N_{0}$ and neglect the effects of thermohaline diffusion on the density of fluid parcels, in which the density perturbation from hydrostatic equilibrium is much smaller than the background density. Additionally, the neglect of Coriolis effects is justified for high-frequency internal waves in the ocean. In the current context, we have the balance of inertial force, buoyancy force, pressure gradient force and viscous force in the momentum equation. We follow the same scalings as in Tabaei & Akylas (Reference Tabaei and Akylas2003) and use a characteristic wavelength $L$ as length scale and $1/N_{0}$ as time scale. Thus, the non-dimensional buoyancy equation, continuity equation and momentum conservation equation are

(2.1) \begin{align} \frac {\partial \rho }{\partial {t}}+\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }\!\rho -\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{k}&=0, \end{align}

(2.2) \begin{align} \frac {\partial \boldsymbol{u}}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}&=-\boldsymbol{\nabla }{\!p}-\rho \boldsymbol{k}+\boldsymbol{\nabla }\boldsymbol{\cdot }\left(\frac {1}{ {\textit{Re}}}\boldsymbol{\nabla }\boldsymbol{u}\right)\!, \end{align}

(2.3) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}&=0, \end{align}

where $\boldsymbol{u}=u\boldsymbol{i}+v\boldsymbol{\!j}+w\boldsymbol{k}$ is the three-dimensional velocity field in the $(x, y, z)$ Cartesian coordinate system. As shown in figure 1, the $y$ axis aligns with the isobaths and the $z$ direction points in the opposite direction of gravity. Here, $\rho$ and $p$ are the reduced density and reduced pressure from hydrostatic equilibrium, respectively. Equation (2.3) indicates that the fluid is incompressible. This assumption is typically used in equations describing oceanic fluid dynamics. Notably, two parameters are used in the subsequent equation derivations. The amplitude parameter $\varepsilon =U/N_0 L$ measures nonlinear effects, where $U$ is a characteristic velocity in the wave beam. Another parameter in (2.2), $({1}/{\textit{Re}})=\nu _*/N_0 L^2$ , is defined as the inverse Reynolds number ( $\nu _*$ denotes the kinematic viscosity).

According to Buhler (Reference Bühler2014), the concept of PV needs to be introduced as

(2.4) \begin{equation} \mathit{\varOmega }=\left (\boldsymbol{\nabla }\times \boldsymbol{u}\right )\boldsymbol{\cdot }\left (\boldsymbol{k}-\boldsymbol{\nabla }\!{\rho }\right )\!, \end{equation}

where $\varOmega$ denotes the non-dimensional PV. Under inviscid fluid conditions, the PV is a material conserved quantity that follows

(2.5) \begin{equation} \frac {\partial \varOmega }{\partial {t}}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\!\varOmega =0. \end{equation}

This property helps us seek the Lagrangian mean flow solution in the perfect fluid (i.e. inviscid fluid) region and establish the relationship between the mean flow and PV theory.

Within the framework of Lagrangian kinematics, the definition of particle displacement needs to be known. Assuming that the position of the marked particle is released at an arbitrary time $t_0$ (called the initial time) is $\boldsymbol{x_0}=(x_0,y_0,z_0)$ , the net displacement of the marked parcel is expressed as $\boldsymbol{\xi }=\boldsymbol{x}-\boldsymbol{x_0}=\xi \!\boldsymbol{i}+\eta \!\boldsymbol{j}+\zeta \!\boldsymbol{k}$ in the ( $x$ , $y$ , $z$ ) Cartesian coordinate system ( $\boldsymbol{x}$ is the position of the fluid parcel at time $t$ ). It follows that the trajectory equation is

(2.6) \begin{equation} \frac {{\textrm d}\xi }{u}=\frac {{\textrm d}\eta }{v}=\frac {{\textrm d}\zeta }{w}={\textrm d}t. \end{equation}

Nonlinearity causes great difficulties in solving the physical quantities of the internal wave and the wave-induced Lagrangian mean flow equations. Assuming weakly nonlinear disturbances, that is, an amplitude parameter $\varepsilon \ll 1$ , we can perform a perturbation expansion for each physical quantity in the governing equation

(2.7) \begin{equation} (\boldsymbol{u},\rho ,p,\varOmega ,\boldsymbol{\xi })=\varepsilon ^1\!\big(\boldsymbol{u^{(0)}},\rho ^{(0)},p^{(0)},\varOmega ^{(0)},\boldsymbol{\xi ^{(0)}}\big)+\varepsilon ^2\big(\boldsymbol{u^{(1)}},\rho ^{(1)},p^{(1)},\varOmega ^{(1)},\boldsymbol{\xi ^{(1)}}\big)+\ldots \end{equation}

Viscous dissipation affects both the structure of the internal wave and the evolution of the Lagrangian mean flow. We assume that the inverse Reynolds number is small $ ({1}/{\textit{Re}} \ll 1 )$ in the exterior region. However, this does not mean that we omit the viscosity terms in the equations. In fact, in the boundary layer, the viscosity has a significant effect on the internal wave and the Lagrangian mean flow, which is shown in our subsequent derivation.

2.2. Linear zeroth-order wave solutions

In this subsection, on the basis of expansion (2.7), we give the governing equations for zeroth-order internal waves in three-dimensional space with weak viscosity and derive the theoretical solutions. Our current research focuses on an internal wave beam incident on a uniform slope of constant angle $\alpha$ to the horizontal at an arbitrary direction. The conceptual framework for defining the geometry of internal wave beams and the $x'y'z'$ -axis is derived from the work of Kataoka & Akylas (Reference Kataoka and Akylas2019), as illustrated in figure 1. Unlike the case of internal waves generated in finite-sized tanks where physical quantities vary in all directions, our study assumes that all derivatives of physical quantities along the constant phase lines of the wave beam plane vanish. Hence, while the fluid is distributed in three-dimensional space, the solution of the equations is fundamentally a two-dimensional problem.

Figure 1.

The rotated coordinate system configuration refers to the set-up formed by Kataoka & Akylas (Reference Kataoka and Akylas2019). Here, $C_{i}$ ( $C_{r}$ ) represents the phase velocity of the incident (reflected) internal wave and $C_{gi}$ ( $C_{gr}$ ) represents the velocity of the incident (reflected) internal wave beam. Also, $\theta$ is the angle between the phase velocity and the $z$ -axis (parallel to the gravity). We establish a right-hand Cartesian coordinate system by setting the intersection of the beam planes of constant phase (blue shaded) and the bottom plane of constant angle $\alpha$ as the $y'$ -axis, defining the line perpendicular to the bottom plane and intersecting with the $y'$ -axis as the $z'$ -axis. The $x'$ -axis lies in the bottom plane and defined as the line perpendicular to the $y'$ - $z'$ plane. Here, $\beta$ represents the angle between the $y'$ -axis and the isobaths of the bottom slope.

Owing to the incompressibility of the fluid, we commonly introduce a streamfunction $\psi (x',z',t )$ to simplify the governing equations

(2.8) \begin{equation} u=\frac {\partial \psi }{\partial z'},\quad w=-\frac {\partial \psi }{\partial x'}. \end{equation}

It can be observed that $\psi (x',z',t )$ inherently satisfies (2.3). Unlike Tabaei, Akylas & Lamb (Reference Tabaei, Akylas and Lamb2005) and Kataoka & Akylas (Reference Kataoka and Akylas2019), the momentum equations we have constructed here consider the effects of viscous dissipation. In fact, the issue of internal wave reflection under viscous conditions was discussed by Kataoka & Akylas (Reference Kataoka and Akylas2020). However, their research focused more on the conditions of the intersection line of the wave beam and the slope parallel to the isobath (two-dimensional, $\beta =0$ ) and paid little attention to the evolution of the mean flow. In a viscous environment, internal waves propagate in an arbitrary direction, reflecting and generating Lagrangian mean flow, which is the focus of our research. Thus, we collect all the $O (\varepsilon )$ terms based on (2.7) and temporarily retain the viscous terms in the governing equations. In the new coordinate system, the governing equations become

(2.9) \begin{align} \frac {\partial {\nabla} ^{2} \psi ^{(0)}}{\partial {t}}&+\frac {\partial \rho ^{(0)}}{\partial {z}'}\sin \alpha \cos \beta -\frac {\partial \rho ^{(0)}}{\partial {x}'}\cos \alpha =\frac {1}{\textit{Re}} {\nabla} ^{4} \psi ^{(0)}, \end{align}

(2.10) \begin{align}& \frac {\partial {v}^{(0)}}{\partial {t}}-\rho ^{(0)}\sin \alpha \sin \beta =\frac {1}{\textit{Re}}{\nabla} ^{2}v^{(0)}, \end{align}

(2.11) \begin{align} \frac {\partial \rho ^{(0)}}{\partial {t}}+v^{(0)}&\sin \alpha \sin \beta +\frac {\partial \psi ^{(0)}}{\partial {x}'}\cos \alpha -\frac {\partial \psi ^{(0)}}{\partial {z}'}\sin \alpha \cos \beta =0, \end{align}

where the gradient operator $\boldsymbol{\nabla }=({\partial }/{\partial x'})\boldsymbol{i}\boldsymbol{'}+( {\partial }/{\partial z'})\boldsymbol{k}\boldsymbol{'}$ in this paper and where $ (\boldsymbol{i}\boldsymbol{'},\boldsymbol{j}\boldsymbol{'},\boldsymbol{k}\boldsymbol{'} )$ denotes the coordinate system composed of unit vectors along the $x'$ , $y'$ and $z'$ -axis directions. The viscosity parameter $ {1}/{\textit{Re}}$ is assumed to be constant and therefore can be extracted outside the operator.

Phillips (Reference Phillips1966) provided the form of solutions for linearised time-harmonic wave reflection on a slope as

(2.12a,b) \begin{align} \psi ^{(0)}=\varPsi _{0} (x',z')e^{-{\textrm{i}}\omega t}&+\textrm{c.c.},\quad v^{(0)}=V_{0}(x',z')e^{-{\textrm{i}}\omega t}+\textrm{c.c.}, \end{align}

(2.13) \begin{align} \rho ^{(0)}&=R_{0}(x',z')e^{-{\textrm{i}}\omega t}+\textrm{c.c.}, \end{align}

where the subscript $0$ represents the zeroth-order first-harmonic term and c.c. denotes complex conjugate. Assuming that the propagation direction of the incident wave beam is oblique to horizontal by $\theta$ (defined as the incidence angle), the dimensionless frequency $\omega$ satisfies

(2.14) \begin{equation} \omega =\sin \theta . \end{equation}

For the inviscid wave reflection problem, the solution consists of two components: incident and reflected waves (Tabaei et al. Reference Tabaei, Akylas and Lamb2005; Kataoka & Akylas Reference Kataoka and Akylas2019). After taking into account the influence of viscosity, the geometric relationship between the incident wave and the reflected wave will change. Therefore, it is necessary to establish the dispersion relation for linear internal waves based on the context of our study. Following the approach of Kataoka & Akylas (Reference Kataoka and Akylas2020), we first give the form of a linear sinusoidal plane wave

(2.15) \begin{equation} (\varPsi _{0}(x',z'), V_{0}(x',z'), R_{0}(x',z'))\sim e^{{\textrm{i}}l \left (x'+mz'\right )}, \end{equation}

where $l$ and $m$ represent the upslope wavenumber and the ratio of the normal to the tangential wavenumber, respectively. By substituting expressions (2.12), (2.13) and (2.15) into the governing equations (2.9)–(2.11), we thus obtain the dispersion relation

(2.16) \begin{equation} \textrm{i}\omega (m^{2}+1)+\frac {(m\sin \alpha \cos \beta -\cos \alpha )^{2}}{\textrm{i}\omega -\dfrac {\sin ^{2} \alpha \sin ^{2} \beta }{\dfrac {1}{\textit{Re}} (m^{2}+1 )l^{2}-\textrm{i}\omega }}=\frac {1}{\textit{Re}} (m^{2}+1)^{2}l^{2}. \end{equation}

Under perfect fluid conditions, the terms on the right-hand side of expression (2.16) tend to zero. In this case, we obtain the following two solutions for $m$ :

(2.17a,b) \begin{align}& m_{\textit{inc}}=\frac {\cos \beta \!\left (\omega r\cos \theta -\sin \alpha \cos \alpha \right )}{\omega ^{2}-\sin ^{2} \alpha },\quad m_{\textit{refl}}=-\frac {\cos \beta \!\left (\omega r\cos \theta +\sin \alpha \cos \alpha \right )}{\omega ^{2}-\sin ^{2} \alpha }, \end{align}

(2.18) \begin{align}&\qquad\qquad\qquad\qquad\qquad r=\frac {\sqrt {\omega ^{2}-\sin ^{2} \alpha \sin ^{2} \beta }}{\omega \cos \beta }. \end{align}

These findings are consistent with those of Kataoka & Akylas (Reference Kataoka and Akylas2019).

When ${1}/{\textit{Re}} \neq 0$ , (2.16) is the sixth power in $m$ and has six roots. Assuming that the inverse Reynolds number is small $ ( {1}/{\textit{Re}} \ll 1 )$ , we find two of the roots in (2.16) by perturbing the inviscid solutions $m_{\textit{inc}}$ and $m_{\textit{refl}}$

(2.19a,b) \begin{align} m_{1}=m_{\textit{inc}}-\textrm{i}\frac {s_{1}}{\textit{Re}},\quad m_{2}=m_{\textit{refl}}+\textrm{i}\frac {s_{2}}{\textit{Re}}. \end{align}

The sign of the imaginary part in (2.19) ensures that the viscosity effect attenuates when the incident wave approaches the slope or when the reflected wave moves away from the slope. Substituting (2.19) into (2.16), $s_{1}$ and $s_{2}$ satisfy

(2.20) \begin{align} s_{1}&=\frac {\omega ({m}_{\textit{inc}}^{2}+1)^{2}l^{2}}{2 (\omega ^{2}m_{\textit{inc}}-\sin ^{2} \alpha m_{\textit{inc}}+\sin \alpha \cos \alpha \cos \beta )}, \end{align}

(2.21) \begin{align} s_{2}&=\frac {-\omega ({m}_{\textit{refl}}^{2}+1)^{2}l^{2}}{2 (\omega ^{2}m_{\textit{refl}}-\sin ^{2} \alpha m_{\textit{refl}}+\sin \alpha \cos \alpha \cos \beta)}. \end{align}

Under weakly viscous conditions, the real parts of $m_{1}$ and $m_{2}$ still represent the ratios of the normal to tangential wavenumbers for incident and reflected waves, respectively, consistent with the results obtained by Kataoka & Akylas (Reference Kataoka and Akylas2019) for inviscid flow. However, in (2.16), there are two other physically meaningful roots that represent disturbances in the boundary layer, which take the form

(2.22a,b) \begin{align} m_{3}={\textit{Re}}^{1/2}\mathit{s_{3}},\quad m_{4}={\textit{Re}}^{1/2}\mathit{s_{4}}, \end{align}

where

(2.23a,b) \begin{align} s_{3}=\frac {\sqrt {2} (1+\textrm{i})}{2l}\sqrt {\omega } ,\quad s_{4}=\frac {\sqrt {2} (1+\textrm{i})}{2l}\sqrt {\frac {\omega ^{2}-\sin ^{2} \alpha }{\omega }}. \end{align}

Two other parameters with opposite signs to (2.23) lead to the streamfunction growing exponentially with increasing $z'$ and thus being ignored. Notably, Kistovich & Chashechkin (Reference Kistovich and Chashechkin1994), Kataoka & Akylas (Reference Kataoka and Akylas2020) and Chashechkin (Reference Chashechkin2021) discussed approximate roots analogous to (2.19) and (2.22) for non-critical reflection cases in the study of the dispersion relation for two-dimensional internal waves under weakly viscous conditions.

Considering that the wave beam is a superposition of sinusoidal waves, by means of the Fourier transformation, we have

(2.24) \begin{align} \varPsi _{0} (x^{\prime},z^{\prime})&=\int _{0}^{\infty }A(l)e^{{\textrm{i}}l \big (x'+m_{1}z'\big )} {\textrm d}l+\int _{0}^{\infty }B(l)e^{{\textrm{i}}l \big (x'+m_{2}z'\big )} {\textrm d}l\nonumber \\ &\quad +\int _{0}^{\infty }C(l)e^{{\textrm{i}}l \big (x'+m_{3}z'\big )} {\textrm d}l+\int _{0}^{\infty }D(l)e^{{\textrm{i}}l \big (x'+m_{4}z'\big )} {\textrm d}l, \end{align}

where $A (l )$ is a prescribed function determined by the specific structure of the incident wave beam. Here, $B (l )$ , $C (l )$ and $D (l )$ are solved by a series of boundary conditions. The first two terms on the right-hand side of (2.24) represent the incident and reflected wave beam affected by viscosity, respectively. Notably, the last two terms in (2.24) represent that the internal waves dissipate exponentially in two ways within the boundary layer. It can be seen that there are three undetermined coefficients in (2.24), which is different from the two-dimensional problem studied by Kataoka & Akylas (Reference Kataoka and Akylas2020). This implies that relying solely on no-slip boundary conditions is insufficient to determine all the coefficients.

We defined the sum of the first two terms and the sum of the last two terms on the right-hand side of (2.24) as ${\varPsi }_{0}^{\textit{out}}$ and ${\varPsi }_{0}^{b}$ , respectively. Similar to expression (2.24), the solutions for the reduced density and transverse velocity also consist of both outer region and boundary layer components in their forms. By performing a unified decomposition of these physical quantities, we obtain

(2.25) \begin{equation} ({\varPsi }_{0}, R_{0},{V }_{0})=\big({\varPsi }_{0}^{\textit{out}}+{\varPsi }_{0}^{b}, {R}_{0}^{\textit{out}}+{R}_{0}^{b}, {V }_{0}^{\textit{out}}+{V }_{0}^{b}\big), \end{equation}

where the superscript $out$ denotes the part of the solution that corresponds to the ideal flow within the outer region, and the superscript $b$ denotes the portion of the wave solution that accounts for boundary layer disturbances. Therefore, we utilise the boundary layer theory in conjunction with a new boundary condition to solve for the unknown coefficients in the linear solution. First, we need to introduce the coordinates within the boundary layer, which can be achieved by stretching the normal coordinate ( $z'$ )

(2.26) \begin{equation} Z=\varepsilon ^{-1/2}z^{\prime}. \end{equation}

In this study, we assume that both nonlinearity and viscosity are weak and equally significant. Consequently, we obtain the following balance:

(2.27) \begin{equation} \frac {1}{\textit{Re}} =\varepsilon . \end{equation}

By substituting (2.26) and (2.27) into the governing equations (2.9)–(2.11) and collecting the leading order, we obtain

(2.28) \begin{align} \qquad -\!\textrm{i}\omega {V}_{0}^{b}-{R}_{0}^{b}\sin \alpha \sin \beta &= \frac {\partial ^{2}{V}_{0}^{b}}{\partial\! Z^{2}}, \end{align}

(2.29) \begin{align} \quad -\!\textrm{i}\omega \frac {\partial ^{2}{\varPhi }_{0}^{b}}{\partial\! Z^{2}}+\frac {\partial {R}_{0}^{b}}{\partial\! Z}\sin \alpha \cos \beta &= \frac {\partial ^{4}{\varPhi }_{0}^{b}}{\partial\! Z^{4}}, \end{align}

(2.30) \begin{align} -\!\textrm{i}\omega {R}_{0}^{b}+{V}_{0}^{b}\sin \alpha \sin \beta -\frac {\partial {\varPhi }_{0}^{b}}{\partial\! Z}\sin \alpha \cos \beta &=0, \end{align}

where ${\varPhi }_{0}^{b}$ is a redefined streamfunction satisfying ${\varPhi }_{0}^{b}=\varepsilon ^{-1/2}{\varPsi }_{0}^{b}$ . Notably, the remaining part of (2.25) can be derived from the ideal fluid governing equations in the outer region, similar to Kataoka & Akylas (Reference Kataoka and Akylas2019).

After a series of algebraic calculations, (2.28) to (2.30) are simplified to a differential equation that contains only one unknown physical quantity ${R}_{0}^{b}$

(2.31) \begin{equation} \frac {\partial ^{3}{R}_{0}^{b}}{\partial\! Z^{3}}+\frac {\textrm{i}\!\left (\omega ^{2}-\sin ^{2} \alpha \right )}{\omega }\frac {\partial {R}_{0}^{b}}{\partial\! Z}=0. \end{equation}

It is significant to note that the boundary layer disturbances are confined within the boundary layer region, which means that, as the boundary layer coordinates $Z\rightarrow \infty$ , ${R}_{0}^{b}\rightarrow 0$ . Under this condition, the form of the solution to (2.31) is as follows:

(2.32) \begin{equation} {R}_{0}^{b}=\int _{0}^{{\infty }}k(l)e^{{\textrm{i}}lx^{\prime}}e^{{\textrm{i}}l{s_{4}}Z} {\textrm d}l. \end{equation}

Inspired by the solution methodology for shallow-water internal solitary waves (Lamb & Yan Reference Lamb and Yan1996), we find a boundary condition to determine $k (l )$ , where the amplitude function of the reduced density satisfies

(2.33) \begin{equation} R_{0}\rightarrow 0\quad {\textrm{as}}\quad z^{\prime}\rightarrow 0. \end{equation}

Since the density perturbation from hydrostatic equilibrium varies with the motion of fluid parcels, boundary condition (2.33) is consistent with the no-slip boundary condition. In the outer region, without considering the viscous terms, the relationship between the reduced density and streamfunction is given by

(2.34) \begin{align} {R}_{0}^{\textit{out}}&=\frac {-\textrm{i}\omega \left (\cos \alpha {\varPsi }_{0x'}^{\textit{out}}-\sin \alpha \cos \beta {\varPsi }_{0z'}^{\textit{out}}\right )}{\omega ^{2}-\sin ^{2} \alpha \sin ^{2} \beta }\nonumber \\ &=\int _{0}^{{\infty }}\omega l\!\left (\cos \alpha -\sin \alpha \cos \beta m_{\textit{inc}}\right )\!A(l)e^{\textrm{i}l \left (x^{\prime}+m_{1}z'\right )} {\textrm d}l\nonumber \\ &\quad +\int _{0}^{{\infty }}\omega l\!\left (\cos \alpha -\sin \alpha \cos \beta m_{\textit{refl}}\right )\!B(l)e^{\textrm{i}l \left (x^{\prime}+m_{2}z'\right )} {\textrm d}l, \end{align}

where the subscripts $x^{\prime}$ and $z^{\prime}$ denote the partial derivatives of the physical quantity in the $x^{\prime}$ and $z^{\prime}$ directions, respectively. Substituting the spectral amplitude of expression (2.32) and (2.34) into the boundary condition (2.33), we obtain the following result:

(2.35) \begin{equation} k(l)=\frac {\omega l\sin \alpha \cos \beta \!\left (m_{\textit{inc}}A(l)+m_{\textit{refl}}B(l)\right )-\omega l\cos \alpha \left (A(l)+B(l)\right )}{\omega ^{2}-\sin ^{2} \alpha \sin ^{2} \beta }. \end{equation}

Having established the boundary layer disturbance part for the reduced density, we substitute the solution (2.32) into (2.29), and the new equation satisfies

(2.36) \begin{equation} \frac {\partial ^{2}{\varPhi }_{0}^{b}}{\partial\! Z^{2}}+\textrm{i}\omega {\varPhi }_{0}^{b} =\frac {\sin \alpha \cos \beta }{\textrm{i}ls_{4} ^{1/2}}{R}_{0}^{b}+c_{1} (x^{\prime})Z+c_{2} (x^{\prime}). \end{equation}

Matching the solution between the boundary layer and the outer region, we can obtain the spectral amplitude of the last term on the right-hand side of (2.24) as

(2.37) \begin{equation} D(l)=-\frac {\sqrt {2}}{2}\left (1-\textrm{i}\right )\frac {\omega ^{3/2}\cos \beta }{\sqrt {{\textit{Re}}(\omega ^{2}-\sin ^{2} \alpha ) }\sin \alpha }k(l). \end{equation}

Therefore, $D (l )$ is ultimately determined by the asymptotic trend in the reduced density near the boundary. Owing to the presence of viscous effect and the rigidity of the slope boundary, restrictions are imposed on the tangential and normal velocities. In this context, we can use the no-slip boundary condition to determine other parameters

(2.38) \begin{equation} \varPsi _{0x'}\quad {\textrm{and}}\quad \varPsi _{0z'}\rightarrow 0\quad {\textrm{as}} \quad z^{\prime}\rightarrow 0. \end{equation}

To establish the relationship with the given function $A (l )$ , we change the boundary condition (2.38) into the following form via the Fourier transform:

(2.39) \begin{align} \qquad A(l)+B(l)+C(l)+D(l)&=0, \end{align}

(2.40) \begin{align} m_{1}A(l)+m_{2}B(l)+m_{3}C(l)+m_{4}D(l)&=0. \end{align}

After a series of calculations and simplifications, we finally obtain the values for each coefficient

(2.41) \begin{align} B(l)&=\left (-1+\frac {m_{\textit{inc}}-m_{\textit{refl}}}{m_{3}}+\frac {m_{4}-m_{3}}{m_{3}}\left (p_{1}-p_{2}\right )\right )\!A(l), \end{align}

(2.42) \begin{align} C(l)&=-A(l)-B(l)-D(l)\quad {\textrm{and}}\quad D(l)=p_{1}A(l)-p_{2}B(l), \end{align}

where

(2.43) \begin{align} p_{1}&=G^{\prime}\!\big (\sin \alpha \cos \beta m_{\textit{inc}}-\cos \alpha \big )\quad {\textrm{and}}\quad p_{2}=G^{\prime}\!\big (\sin \alpha \cos \beta m_{\textit{refl}}-\cos \alpha \big ), \end{align}

(2.44) \begin{align} G^{\prime}&=-\frac {\sqrt {2}}{2} (1-\textrm{i})\frac {l\omega ^{5/2}\cos \beta }{\sin \alpha \sqrt {{\textit{Re}}(\omega ^{2}-\sin ^{2} \alpha ) } (\omega ^{2}-\sin ^{2} \alpha \sin ^{2} \beta)}. \end{align}

The spectral amplitude $B (l )$ indicates that the viscous effect imposes a dual correction on the reflected wave, making it more complex than the two-dimensional reflection (Kataoka & Akylas Reference Kataoka and Akylas2020). Importantly, all coefficients have omitted terms beyond $O ({\textit{Re}} ^{-1/2} )$ . All relationships from (2.41) to (2.43) are substituted into (2.24) to obtain the amplitude function. Using (2.12a ), the analytical expression for the zeroth-order streamfunction can be derived.

Moreover, for the problem of internal wave reflection at an arbitrary angle, there is an extra transverse velocity $v$ , which is different from the case of normal incidence ( $\beta =0$ ). By performing calculations on (2.28) and (2.30) to eliminate the terms involving reduced density, we can obtain the following differential equation:

(2.45) \begin{equation} \big (\omega ^{2}-\sin ^{2} \alpha \sin ^{2} \beta \big ){V}_{0}^{b}+\frac {\partial {\varPhi }_{0}^{b}}{\partial\! Z}\sin ^{2} \alpha \sin \beta \cos \beta =\textrm{i}\omega \frac {\partial ^{2}{V}_{0}^{b}}{\partial\! Z^{2}}. \end{equation}

On the basis of the asymptotic matching condition, we discard the term in the homogeneous solution in which the real part of the exponential function is positive. Since the analytical solution for the streamfunction has already been determined, the form of the solution for (2.45) can be expressed as

(2.46) \begin{align} {V}_{0}^{b}&=\int _{0}^{{\infty }}{\hat {V}}_{0}^{b}(l)e^{{\textrm{i}}lx'}e^{-\frac {\sqrt {2} (1-{\textrm{i}})}{2}\sqrt {\frac {{ (\omega ^{2}-\sin ^{2} \alpha \sin ^{2} \beta )}}{\omega }}Z} {\textrm d}l\nonumber \\ &\quad +\textrm{i}{\textit{lm}}_{3}\cot \beta \int _{0}^{{\infty }}C(l)e^{{\textrm{i}}l (x^{\prime}+s_{3}Z )} {\textrm d}l+{\textrm{i}}{\textit{lm}}_{4}\tan \beta \int _{0}^{{\infty }}D(l)e^{{\textrm{i}}l (x^{\prime}+s_{4}Z)} {\textrm d}l. \end{align}

Considering the influence of viscosity, the transverse velocity $V_{0}$ also satisfies the no-slip boundary condition

(2.47) \begin{equation} V_{0}\rightarrow 0\quad {\textrm{as}}\quad z^{\prime}\rightarrow 0. \end{equation}

Similar to the method used to solve for the reduced density, the relationship between the transverse velocity in the outer region and the streamfunction satisfies

(2.48) \begin{align} {V}_{0}^{\textit{out}}&=\frac {\left (\cos \alpha \frac {\partial {\varPsi }_{0}^{\textit{out}}}{\partial x^{\prime}}-\sin \alpha \cos \beta \frac {\partial {\varPsi }_{0}^{\textit{out}}}{\partial z^{\prime}}\right )\sin \alpha \sin \beta }{\omega ^{2}-\sin ^{2} \alpha \sin ^{2} \beta }\nonumber \\ &=\int _{0}^{\infty }\hat {V}_{A}^{\textit{out}}(l)e^{{\textrm{i}}l (x'+m_{1}z')} {\textrm d}l+\int _{0}^{\infty }\hat {V}_{B}^{\textit{out}}(l)e^{{\textrm{i}}l (x'+m_{2}z')} {\textrm d}l, \end{align}

(2.49) \begin{align} {\hat {V}}_{A}^{\textit{out}}(l)&=\frac {\big(\textrm{i}l\cos \alpha -\textrm{i}{\textit{lm}}_{\textit{inc}}\sin \alpha \cos \beta \big )\sin \alpha \sin \beta }{\omega ^{2}-\sin ^{2} \alpha \sin ^{2} \beta }A(l), \end{align}

(2.50) \begin{align} {\hat {V}}_{B}^{\textit{out}}(l)&=\frac {\big (\textrm{i}l\cos \alpha -\textrm{i}{\textit{lm}}_{\textit{refl}}\sin \alpha \cos \beta \big )\sin \alpha \sin \beta }{\omega ^{2}-\sin ^{2} \alpha \sin ^{2} \beta }B(l). \end{align}

Combining (2.47) to (2.50), we eventually determine the viscous decay part of the solution for the transverse velocity

(2.51) \begin{equation} {\hat {V}}_{0}^{b}(l)=-\textrm{i}{\textit{lm}}_{3}C(l)\cot \beta +\textrm{i}{\textit{lm}}_{4}D(l)\tan \beta -{\hat {V}}_{A}^{\textit{out}}(l)-{\hat {V}}_{B}^{\textit{out}}(l). \end{equation}

Since (2.41) to (2.43) have already given a detailed expression for $B (l )$ , $C (l )$ and $D (l )$ , the complete analytical solutions for the reduced density (2.13) and transverse velocity (2.12b ) are solved.

In this section, on the basis of previous research, we solved the zeroth-order analytical solutions of each physical quantity in the three-dimensional internal wave governing equations in the context of considering viscosity effects. In the next section, we combine the results of the zeroth-order solution to study the higher-order Lagrangian mean flow and discuss its dynamic mechanism through Lagrangian-averaged governing equations.

3.1. Eulerian mean flow equations

For two reasons, we need to introduce the Eulerian mean flow equations in the first subsection. On the one hand, the construction of the Lagrangian-averaged equations relies on the Eulerian-averaged equations because the governing equations (2.1) to (2.3) are established in the background of the Eulerian field. Moreover, in weakly nonlinear systems, a well-defined relationship between Eulerian-averaged physical quantities and the $O (\varepsilon ^{2} )$ approximation of Lagrangian-averaged physical quantities was obtained by Longuet-Higgins (Reference Longuet-Higgins1953). Taking velocity $\boldsymbol{u}$ as an example, it satisfies

(3.1) \begin{equation} \boldsymbol{u}_{\boldsymbol{L}}=\big \langle \big (\boldsymbol{u}^{(\textbf{1})}+\boldsymbol{\xi }^{(\textbf{0})}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}^{(\textbf{0})}\big )\big \rangle , \end{equation}

where the subscript $L$ represents the Lagrangian-averaged term. Regarding a specific physical quantity $f$ , the Eulerian-averaged operator $\langle f \rangle$ means that

(3.2) \begin{equation} \big \langle f(t)\big \rangle =\frac {\omega }{2n\pi }\int _{t_{0}}^{t_{0}+\frac {2n\pi }{\omega }}f(t)\, {\textrm d}t, \end{equation}

where $t_{0}$ is the initial time and the positive integer $n$ denotes the number of periods under study. The time derivative of an arbitrary time-harmonic physical quantity results in zero after time averaging $\langle {\partial f (t )}/{\partial t}\rangle =0$ . The second term on the right-hand side of (3.1) is referred to as `the Stokes drift’ $\boldsymbol{u}_{\boldsymbol{s}}$ , which is entirely determined by the zeroth-order governing equations. Notably, we need to assume $O(n)\lt \varepsilon ^{-2}$ to ensure that the particle displacement induced by the Lagrangian mean flow remains much smaller than the wavelength in our study. If this constraint is violated, all the physical quantities in expression (3.1) would no longer be functions of a fixed spatial field under Eulerian description (Feng Reference Feng1987). On the other hand, we compare the physical differences between the two types of mean flow equations later. Therefore, it is necessary to first clarify the form of the Eulerian mean flow equations.

We collect all the first-order terms and substitute them into the governing (2.1) to (2.3). According to expression (2.7) and the form of linear solution (2.12) and (2.13), all $O (\varepsilon ^{2} )$ physical quantities consist of two parts

(3.3) \begin{equation} f^{(1)}=f_{1}(x',z')+f_{1h} (x',z',t)=F_{1}(x',z')+\big (F_{1h}(x',z')e^{-2{\textrm{i}}\omega t}+\textrm{c.c.}\big ), \end{equation}

where the subscripts $1$ and $1h$ represent the Eulerian-averaged term and the second-harmonic term, respectively. Since the harmonic components are filtered out after periodic averaging, the first-order Eulerian-averaged physical quantities are equivalent to their own amplitude. Therefore, we apply the averaging operator (3.2) to the governing equations to filter out the second-harmonic terms, resulting in the following Eulerian-averaged equations:

(3.4) \begin{align}&\qquad\qquad\qquad \big \langle \boldsymbol{u}^{(\textbf{0})}\boldsymbol{\cdot }\boldsymbol{\nabla }\!\rho ^{(0)}\big \rangle -\big \langle \boldsymbol{u}^{(\textbf{1})}\boldsymbol{\cdot }\boldsymbol{k}\big \rangle =0, \end{align}

(3.5) \begin{align}& \big \langle \boldsymbol{u}^{(\textbf{0})}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}^{(\textbf{0})}\big \rangle =-\big \langle \boldsymbol{\nabla }\!p^{(1)}\big \rangle -\big \langle \rho ^{(1)}\big \rangle \boldsymbol{k}+\left \langle \boldsymbol{\nabla }\boldsymbol{\cdot }\left (\frac {1}{\textit{Re}} \boldsymbol{\nabla }\boldsymbol{u}^{(\textbf{1})}\right )\right \rangle\! , \end{align}

(3.6) \begin{align}&\qquad\qquad\qquad\qquad\quad \big \langle \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}^{(\textbf{1})}\big \rangle =0, \end{align}

where the left-hand side of (3.5) represents the nonlinear stress terms in the Eulerian mean flow equation. Part of nonlinear stress terms can be converted into internal wave-induced body force terms in the Lagrangian-averaged equations, and its physical significance will be explained in the next section. Notably, (3.4) represents the dimensionless internal energy conservation equation after Eulerian averaging. However, the form of this equation will undergo significant changes in the Lagrangian-averaged equations and finally degenerate into a purely kinematic relation with well-defined physical meanings.

3.2. Three-dimensional Lagrangian mean flow equations

Using the Lagrangian mean velocity equation in the $x'$ direction as a prototype, we systematically demonstrate the construction procedure of the Lagrangian mean flow equation, which can be applied similarly to other directions. In § 2.2 of the paper, we obtain the zeroth-order momentum equation in a rotated coordinate system (2.9). Similarly, substituting expansions (2.7) into (2.2), the governing equations for zeroth-order velocity in a rotated coordinate system are

(3.7) \begin{align} \frac {\partial u^{(0)}}{\partial t}+\rho ^{(0)}\sin \alpha \cos \beta &=-\frac {\partial p^{(0)}}{\partial x'}+\boldsymbol{\nabla }\boldsymbol{\cdot }\left (\frac {1}{\textit{Re}} \boldsymbol{\nabla }\!u^{(0)}\right )\!, \end{align}

(3.8) \begin{align} \frac {\partial w^{(0)}}{\partial t}+\rho ^{(0)}\cos \alpha &=-\frac {\partial p^{(0)}}{\partial z'}+\boldsymbol{\nabla }\boldsymbol{\cdot }\left (\frac {1}{\textit{Re}} \boldsymbol{\nabla }w^{(0)}\right )\!. \end{align}

The rotated coordinate system is shown in figure 1. Performing cross-differentiation on (3.7) and (3.8) in combination with expression (2.8), we can virtually obtain (2.9). Starting from the Eulerian-averaged momentum equation (3.5), the Eulerian mean flow equations can be derived for the rotated coordinate system, which are

(3.9) \begin{align} \big \langle \rho ^{(1)}\sin \alpha \cos \beta \big \rangle +\big \langle \boldsymbol{u}^{(\textbf{0})}\boldsymbol{\cdot }\boldsymbol{\nabla }\!u^{(0)}\big \rangle &=-\frac {\partial \big \langle p^{(1)}\big \rangle }{\partial x'}+\left \langle \boldsymbol{\nabla }\boldsymbol{\cdot }\left (\frac {1}{\textit{Re}} \boldsymbol{\nabla }\!u^{(1)}\right )\right \rangle \!, \end{align}

(3.10) \begin{align} \big \langle \rho ^{(1)}\cos \alpha \big \rangle +\big \langle \boldsymbol{u}^{(\textbf{0})}\boldsymbol{\cdot }\boldsymbol{\nabla }w^{(0)}\big \rangle &=-\frac {\partial \big \langle p^{(1)}\big \rangle }{\partial z'}+\left \langle \boldsymbol{\nabla }\boldsymbol{\cdot }\left (\frac {1}{\textit{Re}} \boldsymbol{\nabla }w^{(1)}\right )\right \rangle \!, \end{align}

(3.11) \begin{align} -\big \langle \rho ^{(1)}\sin \alpha \sin \beta \big \rangle +\big \langle \boldsymbol{u}^{(\textbf{0})}\boldsymbol{\cdot }\boldsymbol{\nabla }v^{(0)}\big \rangle &=\left \langle \boldsymbol{\nabla }\boldsymbol{\cdot }\left (\frac {1}{\textit{Re}} \boldsymbol{\nabla }v^{(1)}\right )\right \rangle \!, \end{align}

where the angle bracket notation $\langle \boldsymbol{\cdot }\rangle$ denotes the Eulerian (i.e. periodic) averaging operator. For convenience, since all physical quantities in the study vary in only two directions, in this subsection and Appendices A and B, we define that the displacement and velocity vector contain components solely in the $x'$ and $z'$ directions (e.g. $\boldsymbol{\xi }^{(\textbf{0})}=(\xi ^{ (0 )}$ , $\zeta ^{(0)})$ ).

To eliminate the nonlinear stress term in the Eulerian-averaged equation (3.9), we should obtain its equivalent expression from (3.7). The detailed derivation can be found in Appendix A. According to expression (A4), the Eulerian-averaged nonlinear stress term can be clearly transformed into a sum of three Stokes correction terms: the buoyancy force, the pressure gradient force and the viscous force. By substituting (A4) into (3.9) to eliminate the Eulerian-averaged nonlinear stress term, we thus derive

(3.12) \begin{align} \rho _{L}\sin \alpha \cos \beta &=-\frac {\partial \big \langle p^{(1)}\big \rangle }{\partial x'}-\left \langle \boldsymbol{\xi }^{(\textbf{0})}\boldsymbol{\cdot }\boldsymbol{\nabla }\!\left (\frac {\partial p^{(0)}}{\partial x'}\right )\right \rangle \nonumber \\&\quad +\left \langle \boldsymbol{\nabla }\boldsymbol{\cdot }\left (\frac {1}{\textit{Re}} \boldsymbol{\nabla }\!u^{(1)}\right )\right \rangle +\left \langle \boldsymbol{\xi }^{(\textbf{0})}\boldsymbol{\cdot }\boldsymbol{\nabla }\!\left (\boldsymbol{\nabla }\boldsymbol{\cdot }\left (\frac {1}{\textit{Re}} \boldsymbol{\nabla }\!u^{(0)}\right )\right )\right \rangle \nonumber \\&=\left(-\frac {\partial {p}}{\partial {x}}\right) _{L}+F_{vL}, \end{align}

where $\rho _{L}=\rho ^{(1)}+\boldsymbol{\xi }^{ (\textbf{0} )}\boldsymbol{\cdot }\boldsymbol{\nabla }\!\rho ^{ (0 )}$ is based on relationship (3.1). Here, $\rho _{L}$ also represents the Lagrangian mean reduced density, which reflects the density perturbation from hydrostatic equilibrium that varies with the mean motion of fluid parcels in a weakly nonlinear system. From a mechanical equilibrium perspective, the three terms $-\rho _{L}\sin \alpha \cos \beta$ , $(-({\partial {p}}/{\partial {x}})) _{L}$ and $F_{vL}$ represent the Lagrangian-averaged buoyancy force, the Lagrangian-averaged reduced pressure gradient force and the Lagrangian-averaged viscous force in the $x'$ direction, respectively. Notwithstanding distinct averaging methods, Dong, Bühler & Smith (Reference Dong, Bühler and Smith2020) yields an analogous mechanical balance in the wave-induced GLM mean flow equations.

Equation (3.12) eliminates the Eulerian-averaged nonlinear stress term, but the mechanisms that influence the Lagrangian mean velocity have not been elucidated. To acquire a profound understanding of the internal wave-induced body force that affects the Lagrangian mean velocity, we need to further construct the governing equations. Following the derivation in Appendix A, we obtain the equivalent wave-induced body force terms for two Stokes correction terms $\langle \boldsymbol{\xi }^{ (\textbf{0} )}\boldsymbol{\cdot }\boldsymbol{\nabla }( {\partial p^{ (0 )}}\!/{\partial x'} )\rangle$ and $\langle \boldsymbol{\xi }^{(\textbf{0})}\boldsymbol{\cdot }\boldsymbol{\nabla }(\boldsymbol{\nabla }\boldsymbol{\cdot }(({1}/{\textit{Re}}) \boldsymbol{\nabla }\!u^{ (0 )} ) )\rangle$ on the right-hand side of (3.12). Subsequently, by substituting the Lagrangian mean velocity for the Eulerian mean velocity in (3.12) and rebalancing the equation, we have (A14). Equation (A14) represents the basic form of the internal wave-induced Lagrangian mean flow equation under conditions where the viscosity coefficient is a function of space. If we follow Feng’s (Reference Feng1987) derivation of the tide-induced body force in the Lagrangian mean flow equation, we can further simplify and merge several terms. However, in our study, we assume that the Reynolds number ${\textit{Re}}$ is a constant, which greatly simplifies our equation. The final result is as follows:

(3.13) \begin{align} \rho _{L}\sin \alpha \cos \beta &=-\frac {\partial \big \langle p^{(1)}\big \rangle }{\partial x^{\prime}}+\frac {1}{\textit{Re}} {\nabla} ^{2}u_{L}+\mathit \varPi _{1}, \end{align}

(3.14) \begin{align} \mathit \varPi _{1}&=\mathit \varPi _{N1}+\mathit \varPi _{R1}, \end{align}

(3.15) \begin{align} \mathit \varPi _{N1}&=-\frac {1}{2}\left \langle \frac {\partial }{\partial x^{\prime}}\big (\boldsymbol{\xi }^{(\textbf{0})}\boldsymbol{\cdot }\boldsymbol{\nabla }\!p^{(0)}\big )\right \rangle \!, \end{align}

(3.16) \begin{align} \mathit \varPi _{R1}&=-\frac {1}{2}\left \langle \frac {\partial \xi ^{(0)}}{{\partial x^{\prime} }}\rho ^{(0)}\sin \alpha \cos \beta \right \rangle -\frac {1}{2}\left \langle \frac {\partial \zeta ^{(0)}}{\partial x^{\prime} }\rho ^{(0)}\cos \alpha \right \rangle +\frac {1}{2}\left \langle \frac {\partial \rho ^{(0)}}{\partial x^{\prime} }\xi ^{(0)}\sin \alpha \cos \beta \right \rangle \nonumber \\ &\quad +\frac {1}{2}\left \langle \frac {\partial \rho ^{(0)}}{{\partial x^{\prime} }}\zeta ^{(0)}\cos \alpha \right \rangle -\frac {1}{2{\textit{Re}}}\left \langle \boldsymbol{\xi }^{(\textbf{0})}\boldsymbol{\cdot }\frac {\partial ^{2}}{\partial x^{\prime 2}}\left (\frac {\partial \boldsymbol{u}^{(\textbf{0})}}{\partial x^{\prime} }\right )+\boldsymbol{\xi }^{(\textbf{0})}\boldsymbol{\cdot }\frac {\partial ^{2}}{\partial x^{\prime} \partial z^{\prime} }\left (\frac {\partial \boldsymbol{u}^{(\textbf{0})}}{\partial z^{\prime} }\right )\right \rangle \nonumber \\ &\quad -\frac {1}{2{\textit{Re}}}\left \langle \frac {\partial ^{2}u^{(0)}}{\partial x^{\prime 2}}\frac {\partial \xi ^{(0)}}{\partial x^{\prime} }\right \rangle -\frac {5}{2{\textit{Re}}}\left \langle \frac {\partial \zeta ^{(0)}}{\partial x^{\prime} }\frac {\partial ^{2}u^{(0)}}{\partial x^{\prime} \partial z^{\prime} }\right \rangle -\frac {3}{\textit{Re}} \left \langle \frac {\partial ^{2}u^{(0)}}{\partial x^{\prime} \partial z^{\prime} }\frac {\partial \xi ^{(0)}}{\partial z^{\prime} }\right \rangle \nonumber \\ &\quad -\frac {7}{2{\textit{Re}}}\left \langle \frac {\partial \zeta ^{(0)}}{\partial z^{\prime} }\frac {\partial ^{2}u^{(0)}}{\partial z^{\prime 2}}\right \rangle +\frac {1}{\textit{Re}} \left \langle \frac {\partial ^{2}w^{(0)}}{\partial x^{\prime 2}}\frac {\partial \xi ^{(0)}}{\partial z^{\prime} }\right \rangle +\frac {1}{2{\textit{Re}}}\left \langle \frac {\partial ^{2}w^{(0)}}{\partial x^{\prime 2}}\frac {\partial \zeta ^{(0)}}{\partial x^{\prime} }\right \rangle \!, \end{align}

where $\mathit \varPi _{1}$ is referred to as the internal wave-induced body force in the $x'$ direction, and each term in $\mathit \varPi _{1}$ is related to the displacement of fluid parcels. According to expression (A4), $\langle \boldsymbol{\xi }^{ (\textbf{0} )}\boldsymbol{\cdot }\boldsymbol{\nabla }( {\partial p^{ (0 )}}\!/{\partial x'} )\rangle$ and $-\langle \boldsymbol{\xi }^{(\textbf{0})}\boldsymbol{\cdot }\boldsymbol{\nabla }(\boldsymbol{\nabla }\boldsymbol{\cdot }(( {1}/{\textit{Re}}) \boldsymbol{\nabla }\!u^{ (0 )} ) )\rangle$ constitute part of the nonlinear stress terms. Therefore, following the derivation in Appendix A, $\mathit \varPi _{1}$ is opposite to the sum of a part of the nonlinear stress and the momentum diffusion of the Stokes drift given by (A10); $\mathit \varPi _{1}$ can be divided into two parts according to its rotationality, defined as $\mathit \varPi _{N1}$ and $\mathit \varPi _{R1}$ , respectively. The part $\mathit \varPi _{N1}$ is coupled to the irrotational reduced pressure gradient force, while $\mathit \varPi _{R1}$ is influenced by viscosity and the reduced density field. In (3.16), the first four terms on the right-hand side actually have another expression. The derivation involves integrating (2.11) with respect to time, which establishes the relationship between the zeroth-order reduced density and the displacement of fluid parcels

(3.17) \begin{equation} \rho ^{(0)}=\xi ^{(0)}\sin \alpha \cos \beta +\zeta ^{(0)}\cos \alpha -\eta ^{(0)}\sin \alpha \sin \beta +C_{\rho }, \end{equation}

where $C_{\rho } (x',z' )$ does not vary with time and represents the reduced density of the fluid parcel at the chosen initial moment $t_{0}$ . Expression (3.17) clearly shows that the vertical displacement of the fluid parcel from the initial position affects the value of the dimensionless density perturbation. By substituting this relationship into expression (3.16), the terms related to the reduced density in the internal wave-induced body force can be changed into the product of the particle displacement.

Following the path that led to (3.13), we can similarly derive the Lagrangian mean flow equation in the $z'$ direction via (3.7), (3.8) and (3.10)

(3.18) \begin{align} \rho _{L}\cos \alpha &=-\frac {\partial \big \langle p^{(1)}\big \rangle }{\partial z^{\prime}}+\frac {1}{\textit{Re}} {\nabla} ^{2}w_{L}+\mathit \varPi _{2}, \end{align}

(3.19) \begin{align} \mathit \varPi _{2}&=\mathit \varPi _{N2}+\mathit \varPi _{R2}, \end{align}

(3.20) \begin{align} \mathit \varPi _{N2}&=-\frac {1}{2}\left \langle \frac {\partial }{\partial z^{\prime}}\big (\boldsymbol{\xi }^{(\textbf{0})}\boldsymbol{\cdot }\boldsymbol{\nabla }\!p^{(\textrm{0})}\big )\right \rangle \!, \end{align}

(3.21) \begin{align} \mathit \varPi _{R2}&=-\frac {1}{2}\left \langle \frac {\partial \xi ^{(0)}}{\partial z'}\rho ^{(0)}\sin \alpha \cos \beta \right \rangle -\frac {1}{2}\left \langle \frac {\partial \zeta ^{(0)}}{\partial z'}\rho ^{(0)}\cos \alpha \right \rangle +\frac {1}{2}\left \langle \frac {\partial \rho ^{(0)}}{\partial z'}\xi ^{(0)}\sin \alpha \cos \beta \right \rangle \nonumber \\ &\quad +\frac {1}{2}\left \langle \frac {\partial \rho ^{(0)}}{\partial z'}\zeta ^{(0)}\cos \alpha \right \rangle -\frac {1}{2{\textit{Re}}}\left \langle \boldsymbol{\xi }^{(\textbf{0})}\boldsymbol{\cdot }\frac {\partial ^{2}}{\partial x'\partial z'}\left (\frac {\partial \boldsymbol{u}^{(\textbf{0})}}{\partial x'}\right )+\boldsymbol{\xi }^{(\textbf{0})}\boldsymbol{\cdot }\frac {\partial ^{2}}{\partial z^{\prime 2}}\left (\frac {\partial \boldsymbol{u}^{(\textbf{0})}}{\partial z'}\right )\right \rangle \nonumber \\ &\quad -\frac {1}{2{\textit{Re}}}\left \langle \frac {\partial ^{2}w^{(0)}}{\partial z^{\prime 2}}\frac {\partial \zeta ^{(0)}}{\partial z'}\right \rangle -\frac {5}{2{\textit{Re}}}\left \langle \frac {\partial \xi ^{(0)}}{\partial z'}\frac {\partial ^{2}w^{(0)}}{\partial x'\partial z'}\right \rangle -\frac {3}{\textit{Re}} \left \langle \frac {\partial ^{2}w^{(0)}}{\partial x'\partial z'}\frac {\partial \zeta ^{(0)}}{\partial x'}\right \rangle \nonumber \\ &\quad -\frac {7}{2{\textit{Re}}}\left \langle \frac {\partial \xi ^{(0)}}{\partial x'}\frac {\partial ^{2}w^{(0)}}{\partial x^{\prime 2}}\right \rangle +\frac {1}{\textit{Re}} \left \langle \frac {\partial ^{2}u^{(0)}}{\partial z^{\prime 2}}\frac {\partial \zeta ^{(0)}}{\partial x'}\right \rangle +\frac {1}{2 {\textit{Re}}}\left \langle \frac {\partial ^{2}u^{(0)}}{\partial z^{\prime 2}}\frac {\partial \xi ^{(0)}}{\partial z'}\right \rangle \!, \end{align}

where $\mathit \varPi _{2}$ is the internal wave-induced body force in the $z'$ direction. Analogous to $\mathit \varPi _{N1}$ and $\mathit \varPi _{R1}$ , $\mathit \varPi _{N2}$ and $\mathit \varPi _{R2}$ represent the irrotational and rotational components of $\mathit \varPi _{2}$ , respectively. Through cross-differentiation of momentum equations (3.13) and (3.18) in combination with expression (2.8), we derive the governing equation in terms of the Lagrangian-averaged streamfunction

(3.22) \begin{equation} \frac {\partial \rho _{L}}{\partial z'}\sin \alpha \cos \beta -\frac {\partial \rho _{L}}{\partial x'}\cos \alpha =\frac {\partial \mathit \varPi _{1}}{\partial z'}-\frac {\partial \mathit \varPi _{2}}{\partial x'}+\frac {1}{\textit{Re}} {\nabla} ^{4}\psi _{L}. \end{equation}

In (3.22), only the rotational part of the internal wave-induced body force is involved. On the basis of (2.10) and (3.11), we can derive the Lagrangian mean flow equations in the $y'$ direction. However, there is a certain difference in the derivation processes of the previous two sets of equations. Since the derivative of each physical quantity along the $y'$ direction is zero, a similar simplification of the $-\boldsymbol{\xi }^{ (\textbf{0} )}\boldsymbol{\cdot }\boldsymbol{\nabla }( {\partial p^{ (0 )}}\!/{\partial x'} )$ term on the right side of (3.12) is not necessary. Following the approach outlined above and skipping this step, we thus obtain

(3.23) \begin{align}&\qquad\qquad\qquad -\!\rho _{L}\sin \alpha \sin \beta =\frac {1}{\textit{Re}} {\nabla} ^{2}v_{L}+\mathit \varPi _{3}, \end{align}

(3.24) \begin{align} \mathit \varPi _{3}&=-\left \langle \frac {2}{\textit{Re}} \frac {\partial \xi ^{(0)}}{\partial x'}\frac {\partial ^{2}v^{(0)}}{\partial x^{\prime 2}}\right \rangle -\left \langle \frac {1}{\textit{Re}} \frac {\partial ^{2}\xi ^{(0)}}{\partial x^{\prime 2}}\frac {\partial v^{(0)}}{\partial x^{\prime} }\right \rangle -\left \langle \frac {2}{\textit{Re}} \frac{\partial \zeta ^{(0)}}{\partial x'}\frac {\partial ^{2}v^{(0)}}{\partial x'\partial z^{\prime} }\right \rangle \nonumber \\ &\quad -\left \langle \frac {1}{\textit{Re}} \frac {\partial^{2} \zeta ^{(0)}}{\partial x^{\prime 2}}\frac {\partial v^{(0)}}{\partial z'}\right \rangle -\left \langle \frac {2}{\textit{Re}} \frac {\partial \xi ^{(0)}}{\partial z'}\frac {\partial ^{2}v^{(0)}}{\partial x'\partial z'}\right \rangle -\left \langle \frac {1}{\textit{Re}} \frac {\partial ^{2}\xi ^{(0)}}{\partial z^{\prime 2}}\frac {\partial v^{(0)}}{\partial x'}\right \rangle \nonumber \\ &\quad -\left \langle \frac {2}{\textit{Re}} \frac {\partial \zeta ^{(0)}}{\partial z'}\frac {\partial ^{2}v^{(0)}}{\partial z^{\prime 2}}\right \rangle -\left \langle \frac {1}{\textit{Re}} \frac {\partial ^{2}\zeta ^{(0)}}{\partial z^{\prime 2}}\frac {\partial v^{(0)}}{\partial z'}\right \rangle \!. \end{align}

Through observation, the irrotational internal wave-induced body force, similar to (3.15) and (3.20), disappears. This irrotational force is associated with the Stokes correction to the original pressure gradient force. The assumption that the physical quantities do not vary along the $y'$ direction implies that the pressure gradient force in this direction equals zero. Consequently, only the rotational part of the internal wave-induced body force remains.

Equations (3.13), (3.18) and (3.23) together provide a comprehensive description of the momentum conservation relationships governing internal wave-induced flow in the Lagrangian-averaged framework within the context of the current research. Furthermore, the governing equations also encompass the conservation of internal energy, as expressed in (3.4). In the rotated coordinate system defined in figure 1, (3.4) becomes

(3.25) \begin{equation} \big \langle \boldsymbol{u}^{(\textbf{0})}\boldsymbol{\cdot }\boldsymbol{\nabla }\!\rho ^{(0)}\big \rangle -\big \langle u^{(1)}\sin \alpha \cos \beta \big \rangle -\big \langle w^{(1)}\cos \alpha \big \rangle +\big \langle v^{(1)}\sin \alpha \sin \beta \big \rangle =0. \end{equation}

With relation (3.1), we replace the Eulerian mean velocity in (3.25) with the Lagrangian mean velocity. By substituting the relationship between the reduced density and displacement from (3.17) into (3.25), we derive an expression involving only the streamfunction $\psi$ and the transverse velocity $v$ , which satisfies

(3.26) \begin{equation} \frac {\partial \psi _{L}}{\partial x'}\cos \alpha -\frac {\partial \psi _{L}}{\partial z'}\sin \alpha \cos \beta +v_{L}\sin \alpha \sin \beta =0. \end{equation}

This Lagrangian-averaged expression degenerates from the internal energy conservation equation that includes dimensionless buoyancy into a completely kinematic relationship. It has a clear physical meaning; that is, the Lagrangian mean velocity along the gravitational direction (i.e. $z$ direction) remains constant at zero. This conclusion is also mentioned in the research of three-dimensional internal wave-induced mean flow in an inviscid environment by Kataoka & Akylas (Reference Kataoka and Akylas2019). In the presence of viscosity, this conclusion remains unchanged because the viscosity terms do not affect the original internal energy conservation relationship (2.11). With respect to (2.11), we essentially make a significant assumption that the buoyancy frequency $N$ is a constant and is treated as a fixed time scale. If we assume that background stratification varies in the gravitational direction, it thus satisfies

(3.27) \begin{equation} N^{2}=\kappa (z){N}_{0}^{2}, \end{equation}

where $1/N_{0}$ remains a dimension of time and where $\kappa$ is a coefficient depending on $z$ . In the new configuration, (2.11) and (3.25) undergo changes in their form, resulting in a Lagrangian-averaged equation distinct from (3.26)

(3.28) \begin{equation} \left \langle w _{g} ^{(0)} \zeta _{g} ^{(0)}\frac {\partial {\kappa }}{\partial {z}}\right \rangle + \kappa \left (\frac {\partial \psi _{L}}{\partial x'}\cos \alpha -\frac {\partial \psi _{L}}{\partial z'}\sin \alpha \cos \beta +v_{L}\sin \alpha \sin \beta \right )=0, \end{equation}

where $w _{g} ^{(0)}$ and $\zeta _{g} ^{(0)}$ represent the Eulerian mean velocity and displacement components, respectively, along the gravitational direction. Although $ {\partial {\kappa }}/{\partial {z}} \neq 0$ , we cannot conclude that the vertical component of the Lagrangian mean velocity exists. Conversely, utilising the relationship (2.6) between velocity and displacement, we deduce that the term $\langle w_{g} ^{(0)} \zeta _{g} ^{(0)}\rangle$ equals zero; thereby, (3.28) degenerates into (3.26). Equation (3.28) yields results consistent with the Lagrangian-averaged vertical flow under the Boussinesq approximation as derived by Ou & Maas (Reference Ou and Maas1986). In weakly nonlinear systems, the lowest-order approximation of Lagrangian mean velocity derived in (3.1) describes the time-averaged motion of fluid parcels in the first order equation (3.28). In an adiabatic fluid without energy exchange between external sources and fluid parcels, the Lagrangian mean flow cannot have components normal to isopycnals (Wunsch Reference Wunsch1971). For an equilibrium fluid with vertically stratified (along a $z$ layer) background density, the Lagrangian mean velocity vanishes in the $z$ direction regardless of whether the buoyancy frequency is constant. Under linear stratification conditions, if the intersection line between the internal wave beam plane and the inclined plane is a contour line $ (\beta =0 )$ , the Lagrangian mean flow velocity will be zero in any direction. The differential equation is as follows:

(3.29) \begin{equation} \frac {\partial \psi _{L}}{\partial x'}\cos \alpha -\frac {\partial \psi _{L}}{\partial z'}\sin \alpha =0, \end{equation}

where the boundary condition

(3.30) \begin{equation} \psi _{L}\rightarrow 0\quad {\textrm{as}}\quad z'\rightarrow 0. \end{equation}

Equation (3.29) shows that, in the two-dimensional case, the incompressibility condition requires the Lagrangian mean velocity to satisfy volume conservation. Additionally, constraints on the velocity in the vertical direction lead to the disappearance of the Lagrangian mean flow throughout the entire domain.

To date, we have acquired three equations, (3.22), (3.23) and (3.26), which involve unknown physical quantities $\psi _{L}$ , $\rho _{L}$ and $v_{L}$ . To derive an analytical solution, we combine these three equations through a series of algebraic operations into a single equation that contains only the Lagrangian mean streamfunction, which satisfies

(3.31) \begin{align} \frac {1}{\textit{Re}} {\nabla} ^{4}\psi _{L}&-\frac {2k_{1l}}{\textit{Re}} \frac {\partial ^{4}\psi _{L}}{\partial x^{\prime 3}\partial z^{\prime} }-\frac {2k_{1l}}{\textit{Re}} \frac {\partial ^{4}\psi _{L}}{\partial x'\partial z^{\prime 3}}+\frac {k_{2l}+k_{3l}}{\textit{Re}} \frac {\partial ^{4}\psi _{L}}{\partial x^{\prime 2}\partial z^{\prime 2}}+\frac {k_{2l}}{\textit{Re}} \frac {\partial ^{4}\psi _{L}}{\partial x^{\prime 4}}\nonumber \\&+\frac {k_{3l}}{\textit{Re}} \frac {\partial ^{4}\psi _{L}}{\partial z^{\prime 4}}=\frac {\partial \mathit{\varPi } _{3}}{\partial x^{\prime} }\frac {\cot \alpha }{\sin \beta }-\frac {\partial \mathit{\varPi } _{3}}{\partial z^{\prime} }\cot \beta -\frac {\partial \mathit{\varPi }_{1}}{\partial z'}+\frac {\partial \mathit{\varPi } _{2}}{\partial x'}, \end{align}

where

(3.32) \begin{equation} k_{1l}=\frac {\cos \alpha \cos \beta }{\sin \alpha \sin ^{2} \beta },\quad k_{2l}=\frac {\cos ^{2} \alpha }{\sin ^{2} \alpha \sin ^{2} \beta },\quad k_{3l}=\frac {\cos ^{2} \beta }{\sin ^{2} \beta }. \end{equation}

Importantly, the Lagrangian mean streamfunction for the outer region cannot be acquired through (3.31). In (3.31), all terms containing the undetermined Lagrangian mean streamfunction are viscosity-dependent, which is neglected in the inviscid flow condition. In the present context, the Lagrangian mean streamfunction equation cannot be constructed from the momentum equation without considering the viscous effect. Therefore, we can only use (3.31) to solve for the damping term in the Lagrangian mean streamfunction. Introducing relationships (2.26) and (2.27), we transform (3.31) into a governing equation for the Lagrangian mean flow damping within the boundary layer

(3.33) \begin{equation} \frac {1}{\sin ^{2} \beta } \frac {\partial ^{4}{\psi }_{L}^{b}}{\partial\! Z^{4}}=f_{k} (x',Z), \end{equation}

where the right-hand side forcing term

(3.34) \begin{align} f_{k}&=\varepsilon ^{1/2}\left (\left \langle \frac {\partial \xi ^{(0)}}{\partial x'}\frac {\partial \rho ^{(0)b}}{\partial\! Z}\sin \alpha \cos \beta \right \rangle +\left \langle \frac {\partial \zeta ^{(0)out}}{\partial x'}\frac {\partial \rho ^{(0)b}}{\partial\! Z}\cos \alpha \right \rangle \right )\nonumber \\ &\quad -\varepsilon ^{1/2}\left \langle \frac {\partial \xi ^{(0)b}}{\partial\! Z}\frac {\partial \rho ^{(0)}}{\partial x'}\sin \alpha \cos \beta \right \rangle \nonumber \\ &\quad + \varepsilon ^{1/2}\cot \beta \!\left (2\frac {\partial }{\partial\! Z}\left \langle \frac {\partial \xi ^{(0)b}}{\partial\! Z}\frac {\partial ^{2}v^{(0)b}}{\partial x'\partial\! Z}\right \rangle +\left \langle \frac {\partial ^{3}\xi ^{(0)b}}{\partial\! Z^{3}}\frac {\partial v^{(0)}}{\partial x'}\right \rangle +\left \langle \frac {\partial ^{2}\xi ^{(0)b}}{\partial\! Z^{2}}\frac {\partial ^{2}v^{(0)b}}{\partial x'\partial\! Z}\right \rangle \right )\nonumber \\ &\quad + \cot \beta \!\left (\frac {\partial }{\partial\! Z}\left \langle \frac {\partial ^{2}\zeta ^{(0)b}}{\partial\! Z^{2}}\frac {\partial v^{(0)b}}{\partial\! Z}\right \rangle +2\left \langle \frac {\partial ^{2}\zeta ^{(0)b}}{\partial\! Z^{2}}\frac {\partial ^{2}v^{(0)b}}{\partial\! Z^{2}}\right \rangle +2\left \langle \frac {\partial \zeta ^{(0)}}{\partial\! Z}\frac {\partial ^{3}v^{(0)b}}{\partial\! Z^{3}}\right \rangle \right )\nonumber \\ &\quad + \varepsilon ^{1/2}\left (3\frac {\partial }{\partial\! Z}\left \langle \frac {\partial ^{2}u^{(0)b}}{\partial x'\partial\! Z}\frac {\partial \xi ^{(0)b}}{\partial\! Z}\right \rangle +\frac {1}{2}\left \langle \frac {\partial \xi ^{(0)b}}{\partial\! Z}\frac {\partial ^{3}u^{(0)b}}{\partial x'\partial\! Z^{2}}\right \rangle +\frac {1}{2}\frac {\partial }{\partial x'}\left \langle \frac {\partial ^{2}u^{(0)b}}{\partial\! Z^{2}}\frac {\partial \xi ^{(0)b}}{\partial\! Z}\right \rangle \right )\nonumber \\&\quad -\left (\frac {1}{2}\varepsilon ^{1/2}\left \langle \frac {\partial \xi ^{(0)}}{\partial x'}\frac {\partial ^{3}u^{(0)b}}{\partial\! Z^{3}}\right \rangle +\frac {7}{2}\left \langle \frac {\partial ^{2}\zeta ^{(0)b}}{\partial\! Z^{2}}\frac {\partial ^{2}u^{(0)b}}{\partial\! Z^{2}}\right \rangle +\frac {7}{2}\left \langle \frac {\partial \zeta ^{(0)}}{\partial\! Z}\frac {\partial ^{3}u^{(0)b}}{\partial\! Z^{3}}\right \rangle \right )\!. \end{align}

We aim to focus on the principal balance relationship that the Lagrangian mean flow satisfies in the boundary layer under weakly viscous conditions. Therefore, unlike the body force term on the right-hand side of (3.31), $f_{k}$ in expression (3.34) is the sum of the lowest-order terms selected from the right-hand side of (3.31). It can be inferred from the body force term that the magnitude of ${\psi }_{L}^{b}$ is $O (\varepsilon ^{5/2} )$ , which is theoretically smaller than the Lagrangian mean flow in the outer region ${\psi }_{L}^{\textit{out}}\sim O (\varepsilon ^{2} )$ . This exhibits the same order of magnitude difference between the outer region and boundary layer components of the zeroth-order internal wave streamfunction. The mean flow we study is actually the result of the interaction of zeroth-order internal waves, and since ${\psi }_{L}^{b}$ rapidly decays within the boundary layer, each body force term on the right-hand side of expression (3.34) involves the boundary layer portion of zeroth-order internal waves ( $v^{ (0 )b}$ , $u^{ (0 )b}$ , etc.). Using matching conditions between the boundary layer and the outer region, we can obtain the solution to the differential equation (3.33) as

(3.35) \begin{equation} {\psi }_{L}^{b}(x',Z)=\int _{Z }^{\infty }\int _{\tau _{3} }^{\infty }\int _{\tau _{2} }^{\infty }\int _{\tau _{1} }^{\infty }\sin ^{2} \beta f_{k} (x',\tau _{0}){\textrm d}\tau _{0}\,{\textrm d}\tau _{1}\,{\textrm d}\tau _{2}\,{\textrm d}\tau _{3} . \end{equation}

Expression (3.35) ultimately provides the analytical solution for the viscous decay part of the Lagrangian mean streamfunction, similar to the third and fourth terms on the right-hand side of solution (2.24).

In the context of the mean flow equation for the outer region, where viscous forces and the Coriolis force are not considered, it is impossible to construct an equation governed by the Lagrangian mean velocity (streamfunction). Under this circumstance, to solve for the Lagrangian mean flow in the inviscid region, we must revert to the Eulerian-averaged equations. Indeed, the determination of the inviscid solution for the Lagrangian mean streamfunction induced by internal waves relies on the conservation of PV, as shown in (2.5).

3.3. Potential vorticity and Lagrangian mean flow in the outer region

In the previous section, we provided the exponential decay component in the Lagrangian mean flow solution within the boundary layer. We currently shift our focus back to the outer ideal fluid region. On the basis of the coordinate system established in figure 1, (3.4) to (3.6) are transformed via the following relationship (2.8):

(3.36) \begin{align} \left \langle -\frac {\partial \psi ^{out(1)}}{\partial\! \chi }-v^{out(1)}\sin \alpha \sin \beta \right \rangle &=\big \langle J \big ({R}_{0}^{\textit{out}},{{\varPsi }_{0}^{\textit{out}}}^{*}\big )+\textrm{c.c.}\big \rangle, \end{align}

(3.37) \begin{align} \big \langle \rho ^{out(1)}\sin \alpha \sin \beta \big \rangle &=\big \langle J \big ({V}_{0}^{\textit{out}},{{\varPsi }_{0}^{\textit{out}}}^{*}\big )+\textrm{c.c.}\big \rangle, \end{align}

where ${\partial }/{\partial\! \chi }=\cos \alpha ( {\partial }/{\partial x'})-\sin \alpha \cos \beta ( {\partial }/{\partial z'})$ as mentioned by Kataoka & Akylas (Reference Kataoka and Akylas2019), an additional governing equation for the streamfunction should be obtained through the expressions of the $x'$ and $z'$ directions from (3.5). Nevertheless, the Eulerian-averaged operator filters out terms that vary periodically with time; this governing equation eventually degenerates into the same form as (3.37). Now, the number of governing equations is fewer than the number of unknown physical quantities, making it impossible to determine the Eulerian-averaged streamfunction. Kataoka & Akylas (Reference Kataoka and Akylas2019) suggested that the relationship of PV conservation should be utilised to determine the time-averaged streamfunction. Combining the expansion (2.7) with the form of time-harmonic solutions, we can obtain the series of the PV as

(3.38) \begin{equation} \mathit \varOmega =\varepsilon \big (Q_{0}e^{-{\textrm{i}}\omega t}+\textrm{c.c.}\big )+\varepsilon ^{2}\big (Q_{1h}e^{-2\textrm{i}\omega t}+\textrm{c.c.}\big )+\varepsilon ^{2}Q_{1}+\cdots .\end{equation}

In the absence of viscous effect, (2.9) and (2.10) can be used to determine that the zeroth-order PV is zero. Similarly, from the momentum equations of the first-order second-harmonic waves, we can also obtain that the first-order second-harmonic PV equals zero. Therefore, the lowest-order expression satisfying (2.5) becomes

(3.39) \begin{equation} \frac {\partial\! Q_{1}}{\partial T}+J \big(Q_{1},{\varPsi }_{1}^{\textit{out}}\big )+\big (J \big (Q_{2c},{{\varPsi }_{0}^{\textit{out}}}^{*}\big )+\textrm{c.c.}\big )=0, \end{equation}

where $Q_{1}$ and ${\varPsi }_{1}^{\textit{out}}$ represent the amplitudes of the Eulerian-averaged $O(\varepsilon ^{2})$ PV and streamfunction in the outer region, respectively. Here, $Q_{2c}$ denotes the amplitude of the second-order first-harmonic term and $J (Q_{2c},{{\varPsi }_{0}^{\textit{out}}}^{*} )+\textrm{c.c.}$ is defined as the internal wave dispersion term similar to the tidal dispersion term in Feng (Reference Feng1986). As discussed in § 3, the time derivatives of any Eulerian or Lagrangian-averaged terms are zero. We now introduce a slow time scale $T=\varepsilon ^{-2}t$ , which allows the averaged PV to vary over a longer time scale. This assumption does not change all the $O(\varepsilon )$ wave equations and $O(\varepsilon ^{2})$ mean flow equations if the incident and reflected beams also depend on $T$ . Notably, the Eulerian-averaged operator $\langle \boldsymbol{\cdot }\rangle$ is still defined as the temporal average of physical quantities over $n$ periods, not over the time scale $T$ . In fact, the slow time scale $T$ is much longer than $n$ wave periods, as we assumed $O(n)\lt \varepsilon ^{-2}$ in § 3.

According to expression (2.4), we obtain $Q_{1}$ by collecting all the first-order terms

(3.40) \begin{equation} Q_{1}=-\!\sin \alpha \sin \beta \boldsymbol{\nabla} ^{2}{\varPsi }_{1}^{\textit{out}}+\frac {\partial {V}_{1}^{\textit{out}}}{\partial\! \chi }-\big (J \big ({V}_{0}^{\textit{out}},{{R}_{0}^{\textit{out}}}^{*}\big )+\textrm{c.c.}\big ), \end{equation}

where ${V}_{1}^{\textit{out}}$ is the Eulerian-averaged first-order transverse velocity. The internal wave dispersion term in (3.39) can be transformed into a convective term governed by Stokes drift. This conclusion is derived in Kataoka & Akylas (Reference Kataoka and Akylas2019) through a series of algebraic manipulations of the PV conservation equation. In fact, for any time-averaged physical quantity that meets certain conditions, its convective transport velocity is the Lagrangian mean velocity, with the relevant derivation provided in Appendix B. Therefore, (3.39) becomes

(3.41) \begin{equation} \frac {\partial\! Q_{1}}{\partial T}+J \big (Q_{1},{\varPsi }_{L}^{\textit{out}}\big )=0. \end{equation}

We cannot directly obtain a general solution by solving this fully nonlinear equation evolving with time $T$ . However, we can assign it a suitable initial value that satisfies the conservation relations, resulting in a streamfunction solution that is non-trivial, unlike the two-dimensional case where the entire field’s streamfunction is zero. Considering that the lowest-order PV is zero, that is, the Stokes correction term in the Lagrangian mean PV equals zero, we postulate $Q_{1}=0$ under which the Lagrangian mean flow exists while satisfying (3.41). Combining expression (3.40), we finally obtain a partial differential equation

(3.42) \begin{equation} a_{0}{\varPsi }_{1x^{\prime}x^{\prime}}^{\textit{out}}+2b_{0}{\varPsi }_{1x^{\prime}z^{\prime}}^{\textit{out}}+c_{0}{\varPsi }_{1z^{\prime}z^{\prime}}^{\textit{out}}=f_{0}, \end{equation}

where

(3.43) \begin{equation} a_{0}=\big (\cos ^{2} \alpha +\sin ^{2} \alpha \sin ^{2} \beta \big ),\quad b_{0}=-\!\sin \alpha \cos \alpha \cos \beta ,\quad c_{0}=\sin ^{2} \alpha , \end{equation}

and the forcing term $f_{0}$ on the right side of (3.42) represents

(3.44) \begin{equation} f_{0}=\frac {\textrm{i}\omega }{\omega ^{2}-\sin ^{2} \alpha \sin ^{2} \beta }\frac {\partial ^{2}J \big ({\varPsi }_{0}^{\textit{out}},{{\varPsi }_{0}^{\textit{out}}}^{*}\big )}{\partial\! \chi ^{2}}-\frac {2\textrm{i}\omega \sin ^{2} \alpha \sin ^{2} \beta }{\big (\omega ^{2}-\sin ^{2} \alpha \sin ^{2} \beta \big )^{2}}J \big ({\varPsi }_{0\chi }^{\textit{out}},{{\varPsi }_{0\chi }^{\textit{out}}}^{*}\big ). \end{equation}

Equation (3.42) is formally an inhomogeneous elliptic equation. This can be easily reduced to an ordinary differential equation through the Fourier transform in $x'$ , and then the general solution of the inhomogeneous differential equation can be derived by taking the Laplace transform in $z'$ . The boundary conditions here still utilise the overall tangential and normal Lagrangian mean velocities, which are equal to zero at $z'=0$ . The general solution for the viscous decay part has already been derived in expression (3.35) from the previous section. Therefore, we need to replace the inviscid Eulerian-averaged streamfunction ${\psi }_{1}^{\textit{out}}$ at $z'=0$ with the Lagrangian mean streamfunction through the Euler–Lagrange transformation. The theoretical solution for the streamfunction controlling the outer region is given by

(3.45) \begin{equation} {\psi }_{1}^{\textit{out}}={\varPsi }_{1}^{\textit{out}}=\int _{0}^{{\infty }}{\hat {\varPsi }}_{1}^{\textit{out}} (l,z^{\prime})e^{\textrm{i}lx^{\prime}}{\textrm d}l, \end{equation}

where the spectral amplitude ${\hat {\varPsi }}_{1}^{\textit{out}}$ can be expressed as

(3.46) \begin{equation} {\hat {\varPsi }}_{1}^{\textit{out}}=\left ({\hat {\varPsi }}_{1z^{\prime}}^{\textit{out}}(l,0)+\frac {2\textrm{i}lb_{0}}{c_{0}}{\hat {\varPsi }}_{1}^{\textit{out}}(l,0)\right )\hat {f}_{p1}+{\hat {\varPsi }}_{1}^{\textit{out}}(l,0)\hat {f}_{p2}+\frac {\hat {f}_{0}}{c_{0}}\otimes \hat {f}_{p1}. \end{equation}

The symbol $\otimes$ denotes the convolution operation. More calculation details are provided in Appendix B. Combining expression (2.8) and (3.1), the streamfunction of the Stokes drift takes the form

(3.47) \begin{equation} {\psi }_{s}^{\textit{out}}=-\frac {\textrm{i}}{\omega }J \big ({\varPsi }_{0}^{\textit{out}},{{\varPsi }_{0}^{\textit{out}}}^{\textrm{*}}\big ). \end{equation}

Notably, the streamfunction’s Stokes correction term is derived through treating the streamfunction ${\psi }_{s}^{\textit{out}}$ as an independent physical quantity and, substituting it into expression (3.1), it equals $-({2\textrm{i}}/{\omega })J ({\varPsi }_{0}^{\textit{out}},{{\varPsi }_{0}^{\textit{out}}}^{\textrm{*}} )$ . The desired result is given by expression (3.47) since the boundary conditions essentially impose constraints on the velocity. By combining expressions (3.45) with (3.47), the Lagrangian mean streamfunction in the outer inviscid region is thus obtained

(3.48) \begin{equation} {\psi }_{L}^{\textit{out}}={\psi }_{1}^{\textit{out}}+{\psi }_{s}^{\textit{out}}. \end{equation}

Therefore, adding the terms on the left-hand sides of (3.35) and (3.48), we find the analytical solution for the Lagrangian mean streamfunction over the entire domain. By utilising the Lagrangian-averaged governing equations (3.23) and (3.26), along with the first-order period-averaged equations (3.36) and (3.37) derived in this section, we can readily establish the polarisation relations for the other Lagrangian-averaged physical quantities. Consequently, analytical solutions for all Lagrangian-averaged physical quantities can be derived.

In this section, starting from the concept of Lagrangian-averaged expansion under the weakly nonlinear assumption, we transformed the original first-order Eulerian-averaged equations into Lagrangian-averaged governing equations. By analysing the equilibrium relationship of the three-dimensional Lagrangian-averaged governing equations, in addition to the steady viscous force, buoyancy force and pressure gradient force exerted on the fluid parcel, we determined that the equilibrium relationship is also affected by the internal wave-induced body force. This distinguishes the mean motion of the fluid parcel from the mean flow obtained by performing periodic averaging at each spatial position. On the basis of the boundary layer assumption, we derived the two-part solution of the Lagrangian mean streamfunction. One part was obtained by directly solving the Lagrangian-averaged equations, and the other part needs to be derived with the concept of PV conservation. In the next section, we use numerical methods to observe the spatio-temporal propagation of three-dimensional primary incident and reflected internal waves under the influence of viscosity, and we also present the spatial structure of the Lagrangian mean flow resulting from the self-interaction of internal waves.

4.1. Primary internal waves

We computed the solution of expression (2.24) for different Reynolds numbers as well as the incidence angle and obtained a series of visualisation results. Figure 2 specifically illustrates the evolution of the dimensionless zeroth-order internal wave streamfunction with respect to phase $\phi$ . Since the oblique incidence angle $\beta$ does not affect the geometric relationship between incident and reflected waves in the $x'$ - $z'$ plane, and according to expressions (2.41)–(2.42), $\beta$ has a negligible influence on the magnitude of $\psi ^{(0)}$ . All zeroth-order internal wave results in this subsection are computed at a fixed angle of $\beta = {\pi }/{4}$ . Unlike the numerical results of the internal wave beam velocity field as seen in Rodenborn et al. (Reference Rodenborn, Kiefer, Zhang and Swinney2011), what we present here is a zeroth-order analytical solution with a no-slip bottom boundary. Therefore, it allows us to observe only the incident and reflected wave beams, without the higher-order wave beams. Since the internal wave beam is a superposition of sinusoidal plane waves in the form of expression (2.15), its phase propagation is constrained to be non-transverse (i.e. the phase velocity possesses no component normal to the $x'$ - $z'$ plane). Our study considers non-critical reflection, which results in figures 2(a) and 2(b) showing no significant differences under various Reynolds numbers ${\textit{Re}}$ , especially in the outer region where the viscous effect is weak. This is well reflected in expressions (2.24) and (2.41), where the spectral amplitude of the reflected wave is generally opposite to that of the incident wave, and the viscous correction coefficient $(({m_{\textit{inc}}-m_{\textit{refl}}})/{m_{3}})+ (({m_{4}-m_{3}})/{m_{3}}) (p_{1}-p_{2} )$ is only $O ({\textit{Re}} ^{-1/2} )$ . Even if the incidence angle is altered, as long as critical reflection is avoided, noticeable changes only in the beam geometry can be observed in figures 2(a) and 2(c). As it is a time-harmonic wave, $\phi =0$ and $\phi =0.5\pi$ in figure 2 represent the real and imaginary parts of the streamfunction amplitude $\mathit \varPsi _{0}$ , respectively. From the results in figure 3, it can be found that the imaginary part of $\mathit \varPsi _{0}$ contributes more significantly to the magnitude of the streamfunction.

Figure 3.

Plots of the magnitude of $\psi ^{(0)}$ (the red solid line), the magnitude of the real part of $\mathit \varPsi _{0}$ (the black solid line) and the magnitude of the imaginary part of $\mathit \varPsi _{0}$ (the blue solid line) at each corresponding position of the marked points in figure 2(a). Each point is labelled with its respective $x'$ -coordinate. The settings are the same as those in figure 2(a).

However, the situation within the boundary layer is different. The last two terms in expression (2.24) cause the streamfunction to undergo abrupt changes within the thin boundary layer. Therefore, we will also focus on the plots of the zeroth-order velocity in the boundary layer (figure 4). These results play a decisive role in the solution of the Lagrangian governing equations, as shown in (3.33). One part of the Lagrangian mean streamfunction ${\psi }_{L}^{b}$ essentially balances with the body force term, which rapidly decays in the boundary layer, as expressed in (3.34). The variation patterns of the Lagrangian mean streamfunction ${\psi }_{L}^{b}$ with different parameters are described in detail in the next section.

Figure 4.

The normalised zeroth-order velocity solution $u^{(0)}$ in the boundary layer region: (a) ${\textit{Re}} =900$ , $\theta = {\pi }/{4}$ ; (b) ${\textit{Re}} =900$ , $\theta ={2\pi }/{9}$ ; (c) ${\textit{Re}} =900$ , $\theta = {\pi }/{5}$ . The remaining settings are the same as those in figure 2.

4.2. The Lagrangian mean flow decays within the boundary layer

In this subsection, the image results of the solution to the Lagrangian-averaged equation (3.33) are demonstrated. In contrast to the previous subsection, we analyse the variation in the mean flow within the boundary layer from three aspects: the Reynolds number (figure 5), the incidence angle (figure 6) and the oblique incidence angle (figure 7).

As demonstrated in expression (2.22), the boundary layer thickness that causes the rapid decay of the zeroth-order internal wave depends on the Reynolds number ${\textit{Re}}$ . This is the reason why the normalised Lagrangian mean streamfunction ${\psi }_{L}^{b}$ exhibits different vertical extension ranges under different Reynolds numbers, as shown in figure 5. As the viscous effect intensifies, the boundary layer of the mean flow gradually thickens, ultimately showing a tendency to extend throughout the entire region. Under strong viscous conditions, it is necessary to re-derive the solutions of the governing equations for internal waves and Lagrangian mean flow. The strong viscous effect was discussed in the study of tide-induced Lagrangian mean flow (Feng Reference Feng1987). In the study of Feng (Reference Feng1987), the solution of the tide-induced Lagrangian-averaged equation is applicable to the entire region, while in our study, (3.33) can provide only the solution for mean flow damping within the thin boundary layer.

Figure 5.

The normalised viscous decay streamfunction ${\psi }_{L}^{b}$ (the solution of (3.33)) in the boundary layer: (a) ${\textit{Re}} =900$ , $\theta = {\pi }/{4}$ ; (b) ${\textit{Re}} =225$ , $\theta = {\pi }/{4}$ . The remaining settings are the same as those in figure 2.

Figure 5 revealed that the difference between the incidence angle of the wave beam and the bottom slope angle significantly affects the Lagrangian mean flow in the boundary layer. Consistent with figure 4, three sets of experiments in figure 6 were conducted with different incidence angles. As $\left | \theta -\alpha \right |$ decreases, it can be observed that the Lagrangian mean flow significantly strengthens within the boundary layer. Similar to the mechanism of near-critical reflection in internal waves, a decrease in $\left | \theta -\alpha \right |$ causes both the real and imaginary parts of the two coefficients in (2.19a , b ) to increase, which has a particularly significant impact on $m_{2}$ . This variation can also be clearly seen in figure 4, where the wave velocity $u^{ (0 )}= {\partial \psi ^{ (0 )}}\!/{\partial z'}$ within the boundary layer shows a noticeable increasing trend as the incidence angle decreases. This result inevitably leads to an increase in the internal wave-induced body force in (3.33), thereby causing changes in the Lagrangian mean streamfunction on the left side of the equation. The enhancement of the mean flow is not without limits; near-critical reflection occurs when $ | \sin ^{2} \theta -\sin ^{2} \alpha | \sim O (\varepsilon )$ , necessitating a reconsideration of the balance relationship in the zeroth-order internal wave equations. This provides inspiration, suggesting that future research could focus on the Lagrangian mean flow induced by internal waves under near-critical reflection circumstances. It is foreseeable that ${\psi }_{L}^{b}$ will play a dominant role within the boundary layer.

Figure 6.

The normalised viscous decay streamfunction ${\psi }_{L}^{b}$ (the solution of (3.33)) in the boundary layer: (a) ${\textit{Re}} =900$ , $\theta = {\pi }/{4}$ ; (b) ${\textit{Re}} =900$ , $\theta = {2\pi }/{9}$ ; (c) ${\textit{Re}} =900$ , $\theta = {\pi }/{5}$ . The remaining settings are the same as those in figure 2.

Figure 7.

The normalised viscous decay streamfunction ${\psi }_{L}^{b}$ (the solution of (3.33)) in the boundary layer: (a) ${\textit{Re}} =900$ , $\theta = {\pi }/{4}$ , $\beta = {\pi }/{12}$ ; (b) ${\textit{Re}} =900$ , $\theta = {\pi }/{4}$ , $\beta = {\pi }/{4}$ ; (c) ${\textit{Re}} =900$ , $\theta = {\pi }/{4}$ , $\beta = {5\pi }/{12}$ . The remaining settings are the same as those in figure 2.

Since the oblique incidence angle $\beta$ only appears in $O (\varepsilon ^{{1}/{2}} )$ terms of the spectral amplitudes of the zeroth-order internal wave streamfunction, it has no significant effect on the magnitude of $\psi ^{(0)}$ . However, according to expression (3.35), there is a strong correlation between reflection obliqueness and the intensity of Lagrangian mean flow, which is also shown in figure 7. In figures 7(a) and 7(b), it can be found that the magnitude of the Lagrangian mean streamfunction maintains consistency with variations in $\beta$ , which can be explained by solution (3.35). Notably, as shown in figures 7(b) and 7(c), when $\beta \geqslant {\pi }/{4}$ , a limited number of positions exhibit a significant reduction in the magnitude of $\psi _{L}^{b}$ . This is because the product of $\textrm{sin}^{2}\beta$ and part of the terms in expression (3.34) begin to decrease when $\beta$ varies within $ {\pi }/{4}\sim {\pi }/{2}$ . Generally, the Lagrangian mean flow within the boundary layer strengthens with increasing $\beta$ . In fact, § 3.2 derives that $\psi _{L}$ vanishes in two-dimensional space, the wave-induced Lagrangian mean flow fundamentally depends on three-dimensional wave reflection. This conclusion has crucial implications for the spatial configuration of the internal wave-induced Lagrangian mean flow.