On physical optics approximation of viscous stratified shear flows

TL;DR

This paper develops a mathematical model for viscous stratified shear flows, deriving a modified Taylor-Goldstein equation and applying geometric and physical optics approximations to analyze wave behavior.

Contribution

It introduces a modified TG equation incorporating vertical and horizontal viscosity effects and applies WKB methods to study asymptotic solutions near singularities.

Findings

Derived a fourth-order modified TG equation with viscosity effects.

Applied WKB method to analyze wave solutions in stratified shear flows.

Investigated wave behavior near turning points and critical levels.

Abstract

We study a mathematical model of a perturbed stratified shear mean flow in the presence of eddy coefficients of turbulent viscosity. We adopt the standard Boussinesq approximation in the natural convection of the buoyancy-driven flow and neglect the influence of the eddy coefficients of turbulent diffusivity. Comprising both the vertical and horizontal viscosity effects, a model for the vertical velocity perturbation corresponds to a fourth-order Taylor-Goldstein (TG) differential equation. Considering only the latter, we obtained a modified TG equation with the same order as the classical, inviscid one. Under an assumption of the slowly varying Brunt-V\"ais\"al\"a frequency and background horizontal velocity, we discuss the corresponding geometrical and physical optics approximations of the modified TG equation using the WKB method. We further investigate the behavior of these asymptotic solutions near singular values of a turning point and critical level.

Equations

$$\frac{\partial \tilde{\mathbf{u}}}{\partial t}+\left(\tilde{\mathbf{u}}\cdot \nabla \right)\tilde{\mathbf{u}}=-\frac{1}{{\rho }_0}\nabla \tilde{p}+\tilde{\theta }{\mathbf{e}}_z+\frac{\partial }{\partial x}\left({\tilde{a}}_H\frac{\partial \tilde{\mathbf{u}}}{\partial x}\right)+\frac{\partial }{\partial z}\left({\tilde{a}}_T\frac{\partial \tilde{\mathbf{u}}}{\partial z}\right).$$

$$\frac{\partial \tilde{\theta }}{\partial t}+\left(\tilde{\mathbf{u}}\cdot \nabla \right)\tilde{\theta }=\frac{\partial }{\partial x}\left({\tilde{k}}_H\frac{\partial \tilde{\theta }}{\partial x}\right)+\frac{\partial }{\partial z}\left({\tilde{k}}_T\frac{\partial \tilde{\theta }}{\partial z}\right),$$

$$\tilde{u}=U(z),\tilde{v}=0,\tilde{w}=0,\tilde{\theta }=\Theta (z),\text{and}\tilde{p}=P(z).$$

$$\left(\frac{\partial }{\partial t}+U\frac{\partial }{\partial x}\right)\theta +w\frac{d\Theta }{dz}$$

$$\left(\frac{\partial }{\partial t}+U\frac{\partial }{\partial x}\right)u+w\frac{dU}{dz}$$

$$\left(\frac{\partial }{\partial t}+U\frac{\partial }{\partial x}\right)w$$

$$\frac{\partial u}{\partial x}+\frac{\partial w}{\partial z}$$

$${\left(\frac{\partial }{\partial t}+U\frac{\partial }{\partial x}\right)}^2{\nabla }^{2}w$$

$$=\left(\frac{\partial }{\partial t}+U\frac{\partial }{\partial x}\right)\left[\frac{\partial }{\partial z}\left(A_N\frac{{\partial }^{3}w}{\partial z\partial x^2}\right)+A_H\frac{{\partial }^{4}w}{\partial x^4}+\frac{{\partial }^2}{\partial z^2}\left(A_T\frac{{\partial }^{2}w}{\partial z^2}\right)\right],$$

$${\hat{w}}^{\prime \prime }+\left[\frac{N^2}{(U-c{)}^2}-\frac{U^{\prime \prime }}{U-c}-k^2\right]\hat{w}=\frac{1}{ik(U-c)}\left[\frac{d^2}{d{z}^2}\left(A_T{\hat{w}}^{\prime \prime }\right)-k^2\frac{d}{dz}\left(A{\hat{w}}^{\prime }\right)+k^4{A}_H\hat{w}\right],$$

$${\hat{w}}^{\prime \prime }+Q_1(z){\hat{w}}^{\prime }+Q_0(z)\hat{w}=0,$$

$$Q_0(z)$$

$$Q_1(z)$$

$$\hat{w}(z)=\text{exp}\left(-\frac{1}{2}{\int }^z{Q}_1(\zeta )d\zeta \right)\tilde{w}(z),$$

$${\tilde{w}}^{\prime \prime }+Q_2(z)\tilde{w}=0,$$

$$Q_2(z)=Q_0-\frac{Q_1^2}{4}-\frac{Q_1^{\prime }}{2}.$$

$${\epsilon }^2\ddot{\tilde{w}}+Q_{\epsilon }(Z)\tilde{w}=0,$$

$$Q_{\epsilon }(Z)=\frac{N^2}{u_1{u}_2}-k^2-\frac{{\epsilon }^2}{u_2}\left(\frac{u_3}{u_2}+u_4\right),$$

$$u_3$$

$$u_4$$

$$Q_{\epsilon }(Z)=A(Z)+iB(Z)-{\epsilon }^2\left[C(Z)+iD(Z)\right]$$

$$A(Z)$$

$$B(Z)$$

$$\tilde{C}(Z)$$

$$+\frac{1}{2}\ddot{A_H}\left[k(U-c_r{)}^2\left(2c_i-k{A}_H\right)+k^2{c}_i^2{A}_H\right]$$

$$C(Z)$$

$$\tilde{D}(Z)$$

$$-\frac{1}{2}\left[(U-c_r{)}^2-c_i^2\right]\left[2\ddot{U}(c_i+k{A}_H)-k\dot{U}\dot{A_H}+k(U-c_r)\ddot{A_H}\right]$$

$$D(Z)$$

$$M(Z)$$

$$=m_0(Z)+i{n}_0(Z)+\epsilon \left[m_1(Z)+i{n}_1(Z)\right]+{\epsilon }^2\left[m_2(Z)+i{n}_2(Z)\right]+\dots .$$

$$\tilde{w}(Z)=\tilde{w}(Z=0)\text{exp}\left(\frac{1}{\epsilon }{\int }_0^{Z}M(\zeta )d\zeta \right).$$

$$m_0^2-n_0^2+A$$

$$2m_0{n}_0+B$$

$$m_0^2(Z)$$

$$n_0^2(Z)$$

$$\text{exp}\left\{\frac{1}{\epsilon }{\int }_0^Z\left[m_0(\zeta )+i{n}_0(\zeta )\right]d\zeta \right\}\text{and}\text{exp}\left\{-\frac{1}{\epsilon }{\int }_0^Z\left[m_0(\zeta )+i{n}_0(\zeta )\right]d\zeta \right\}.$$

$${\tilde{w}}^{\prime \prime }+\left[(A+iB)\mp (m_0^{\prime }+i{n}_0^{\prime })\right]\tilde{w}=0,$$

Full text

\Title

On physical optics approximation of viscous

 stratified shear flows \Author  Natanael Karjanto \orcidA \AuthorNamesNatanael Karjanto

1 Introduction

We consider dynamics of internal gravity waves in the atmosphere for continuously stratified parallel flows in the presence of horizontal eddy coefficients of turbulent viscosity. Under the Boussinesq approximation, we derive a vertical velocity perturbation equation from the two-dimensional Euler equations as the governing model. By including both the horizontal and vertical eddy viscosity, an equation for the vertical velocity perturbation is given by a fourth-order ordinary differential equation (ODE). On the other hand, in the absence of vertical eddy viscosity, it is modeled by a second-order ODE, an extended version of the classical Taylor-Goldstein (TG) equation for the inviscid case Taylor (1931); Goldstein (1931); Haurwitz (1931). The latter governs the dynamics of Kelvin-Helmholtz (KH) instability, which in general refers to the instability for continuous distribution of velocity and density Kelvin (1871); Helmholtz (1868).

The study of flow with viscosity has a long history, starting from the Orr-Sommerfeld equation, an eigenvalue equation that models two-dimensional linear stability problem for nearly parallel viscous flows in a straight channel or a boundary layer Orr (1907a, b); Sommerfeld (1909). An extended version of the TG equation to include the effects of constant viscosity and diffusivity is modeled by Koppel using a sixth-order ODE Koppel (1964). Hazel solved numerically the higher-order equation for the vertical velocity of internal gravity waves with kinematic viscosity and thermal conductivity Hazel (1967). A couple of years later, the infinitesimal stability of inviscid, parallel, stratified shear flows described by the TG equation was studied numerically Hazel (1972). For a mathematical perspective of instability in viscous fluids, see Renardy and Renardy (2003).

The TG equation and its variants have been used as mathematical models in various physical settings where stratified fluids or density stratification and shear flow occur. The interest lies in hydrodynamic stability since it sheds light on the problem of transition from the laminar to turbulent flow. Some examples where the TG equation aids in understanding this stability include surface gravity waves in the atmosphere Moninger and Gossard (1976); Janssen and Komen (1985a, b), internal gravity waves in a stably stratified boundary layer Anne and Driedonks (1985); Chimonas and Nappo (1989), and internal waves in the ocean Sun et al. (1998); Watson (1994); Thorpe (1999), cf. Hall (1993); Lawrence et al. (1991); Farrell and Ioannou (1993); Sutherland (2005); Zhou and Morton (2007); Tedford et al. (2009); Camassa and Viotti (2013); Pütz (2018).

There exists, however, theoretical studies of the TG equation. For instance, Baldwin and Roberts studied the critical layer in a stratified shear flow where they introduced small viscosity and heat conductivity to remove the singularity Baldwin and Roberts (1970). Banerjee presented a geometry of the TG equation and some generalized results related to the stability Banerjee (2005). Using a variational approach, Hirota and Morrison developed an efficient method for finding marginally stable (neutral mode) solutions of the TG equation Hirota and Morrison (2015, 2016).

The interest in studying the TG model has been growing for the past two decades. Pratt and collaborators proposed an extended TG equation for modeling the stratified shear flow of arbitrary cross-section Pratt et al. (2000); Deng et al. (2003). While these authors also formulated the corresponding semicircle theorem, Subbiah and Ganesh sharpened it by deriving a semi-ellipse theorem for unstable modes. In addition to estimating their growth rate, they also obtained the upper and lower bounds for the neutral modes phase velocity Subbiah and Ganesh (2008); Ganesh and Subbiah (2009).

After the classical works on the TG equation with viscosity in the 1970s Hazel (1972); Baldwin and Roberts (1970), it was not until the 1990s that the work starts to blossom again. Hall discussed the nonlinear instability of a boundary layer on a heated flat plate positioned in an incoming flow Hall (1992). Dando examined the nonlinear evolution of inviscid Görtler vortices in three-dimensional boundary layers Dando (1996). Ng and Reid considered the stability problem in plane Poiseuille flow with the combined effects of stratification and viscosity Ng and Reid (1997). Mureithi et al. investigated the effect of buoyancy on the upper-branch linear stability characteristics of an accelerating boundary layer flow in Tollmien-Schlichting waves Mureithi et al. (1997). See also Watson and Otto (2001); Hogg et al. (2001, 2002); Thorpe (2012); Lawrence et al. (2013); Xiao et al. (2015).

Fast forward to the past decade, Smyth et al. investigated properties of shear instabilities with small-scale turbulence represented by vertical eddy viscosity and diffusivity coefficients Smyth et al. (2011). Liu et al. examined an extended TG equation by including the effects of small-scale turbulence on the stability of stratified shear flows Liu et al. (2012). Thorpe et al. investigated the marginal stability of a stably stratified shear flow with ambient turbulence represented by vertical small eddy viscosity and diffusivity Thorpe et al. (2013). In some circumstances, they found that viscosity amplifies instability. Li et al. further examined the effects of vertically varying eddy viscosity and diffusivity on the KH instability of a stratified shear flow Li et al. (2015). Additionally, Khani and Waite measured horizontal and vertical eddy viscosity in stratified turbulence using high-resolution simulations Khani and Waite (2013). Recently, Lian et al. performed numerical computations for the stability analysis of stratified and parallel shear flows with the effects of small-scale turbulence described by eddy viscosity and diffusivity Lian et al. (2020).

Our focus in this article is a modified TG equation with the presence of horizontal eddy viscosity but the absence of vertical eddy viscosity. Under the assumption of slowly varying background horizontal velocity and Brunt-Väisälä frequency, we propose an asymptotic solution for this modified TG equation using the WKB method. The technique is named after and developed by physicists Gregor Wentzel, Hendrik Anthony Kramers, and Léon Brillouin in 1926. They refined the ideas of George Green and Joseph Liouville Gill (1982); Steele (1976). Although mathematician Harold Jeffreys had developed a general method for approximating solutions of second-order linear ODEs in 1923, his contribution is often neglected since the trios WKB were not really aware of this earlier work. Indeed, a better name for this method should be coined as the JWKB or WKBJ approximation Bender and Orszag (2010); Baines (1995). There are in fact at least 12 people who contributed to the development of WKB theory, including Liouville, Green, Horn, Rayleigh, Gans, Jeffrey, Wentzel, Kramers, Brillouin, Langer, Olver, and Meyer Hinch (1991).

The WKB method has some limitations. The so-called physical optics WKB approximation turns singular and vanishes at a turning point and critical level, respectively. In the vicinity of the critical level, i.e., a level where the mean fluid velocity is equal to the horizontal phase velocity, an internal gravity wave flow is unstable Booker and Bretherton (1967); Geller et al. (1975); Fritts and Geller (1976). The behavior of both 2D and 3D internal gravity wave packets approaching a critical level has been examined both analytically and numerically Winters and D’Asaro (1989); Winters and Riley (1992); Winters and D’Asaro (1994). Shear instability is responsible for maintaining deep-cycle turbulence due to the extraction of energy from the mean flow near the critical level through the interaction of mean shear and subgrid-scale stress Sun et al. (1998). Note that although this reference refers to the latter as Reynolds stress, mean shear generally interacts with non-resolving fluxes, which are not necessarily Reynolds stress.

We organize this article as follows. The first part of Section 2 presents a derivation of mathematical modeling for vertical velocity perturbation with eddy coefficients of turbulent viscosity. Subsection 2.2 focuses on a modified TG equation in the presence of horizontal coefficient eddy viscosity only. Section 3 derives both the geometrical and physical optics approximation using the WKB method. Section 4 presents the solution near a turning point and critical level for different values of wind shears. Section 5 concludes from our discussion.

2 Mathematical modeling

2.1 An extended Taylor-Goldstein equation

Consider the Navier-Stokes momentum equation for an incompressible flow of the Newtonian fluid Batchelor (2000); White (2006); Kundu and Cohen (2015):

[TABLE]

The fluid motion satisfies this equation and the continuity equation $\nabla\cdot\tilde{\mathbf{u}}=0$. For each quantity, we write $q=Q+\delta q^{\prime}$, where $\delta\ll 1$ is a small parameter, $Q$ denotes a background state, undisturbed flow, taken to be steady of slowly varying and horizontally uniform, but varying in the vertical direction, of the corresponding boundary layer; and $\delta q^{\prime}$ is the first-order perturbation. Assume that the background flow is in the hydrostatic balance.

For another small parameter $\eta\ll 1$, $\tilde{\theta}=\eta g\tilde{\theta}^{\prime}/\Theta$ is the potential temperature. Furthermore, $g$ is the gravitational acceleration, $\tilde{p}$ is the pressure, $\rho_{0}$ is the reference density, $\tilde{a}_{H}(x,z,t)$ and $\tilde{a}_{T}(x,z,t)$ are the horizontal and the vertical eddy coefficients of turbulent viscosity, respectively. The continuity equation for potential temperature is given by

[TABLE]

where $\tilde{k}_{H}(x,z,t)$ and $\tilde{k}_{T}(x,z,t)$ are the horizontal and the vertical eddy coefficients of turbulent diffusivity, respectively. In what follows, we will assume that $\tilde{a}_{H}$, $\tilde{a}_{T}$, $\tilde{k}_{H}$, and $\tilde{k}_{T}$ are independent of $x$, following the argument in Liu et al. (2012); Thorpe et al. (2013); Li et al. (2015), cf. Kundu and Cohen (2015); Garratt (1992); Kantha and Clayson (2000). In general, however, they should depend on $(x,z,t)$, as shown recently that even at anisotropic grid points with $x\gg z$, both eddies in horizontal and vertical directions are in fact coupled Khani and Waite (2020).

The Navier-Stokes equation (1) and the continuity equation for potential temperature (2) admit the initial and boundary conditions $\tilde{\mathbf{u}}(x,z,0)=\tilde{\mathbf{u}}_{0}(x,z)$, and $\tilde{\mathbf{u}}(0,z,t)=\mathbf{f}_{1}(z,t)$, $\tilde{\mathbf{u}}(x,0,t)=\mathbf{f}_{2}(x,t)$, respectively. They also possess an exact stationary solution for the flow velocity $\tilde{\mathbf{u}}=(\tilde{u},\tilde{v},\tilde{w})$, potential temperature, and pressure, where the former is written in the component form as follows:

[TABLE]

Here, $U(z)$, $\Theta(z)$, and $P(z)$ denote the horizontal mean flow, potential temperature, and pressure at the basic state, respectively. This indicates that the equilibrium of a stratified fluid in gravity field can be achieved when these quantities depend only on the vertical height $z$. Hence, the flow velocity vector with first-order perturbation can be written as $\tilde{\mathbf{u}}=(U(z)+\delta u,0,\delta w)$.

The linearized equations in the presence of eddy viscosity and diffusivity terms read

[TABLE]

We can exclude the pressure terms by subtracting $\partial_{z}$ of (4) from $\partial_{x}$ of (5). Applying a derivative operator $\left(\frac{\partial}{\partial t}+U\frac{\partial}{\partial x}\right)\frac{\partial}{\partial x}$ to this equation, and employing (3) and (6), we can further eliminate $u$ and $\theta$ to obtain a single equation for $w$. In the absence of diffusivity, it is given as follows:

[TABLE]

where $\nabla^{2}=\partial_{x}^{2}+\partial_{z}^{2}$ is a Laplacian operator, $N^{2}=g\frac{d}{dz}\ln\Theta$ is the square of the Brunt-Väisälä (buoyancy) frequency, and $A_{N}=A_{H}+A_{T}$ Holton (2004); Nappo (2013). Substituting $w(x,z,t)=\hat{w}(z)e^{ik(x-ct)}$ into (7), we obtain a fourth-order ODE in $\hat{w}$, where the prime denotes differentiation with respect to $z$:

[TABLE]

with boundary conditions $\hat{w}(0)=\hat{w}_{0}$ and $\hat{w}^{\prime}(0)=\hat{w}_{1}$. Note that the absence of the viscosity terms on the right-hand side of (8) leads to the classical Taylor-Goldstein equation, which governs the behavior of a perturbed vertical velocity in continuously stratified parallel flows Taylor (1931); Goldstein (1931). For unidirectional flow, the Miles-Howard semicircle theorem is established Miles (1961); Howard (1961). It states that the complex-valued wave velocity $c$ of any unstable mode of a perturbed inviscid fluid parallel flows must lie inside the semicircle in the upper half of the $c$-plane with a diameter the range of $U$ Huppert (1973); Engevik et al. (1985); Miles (1986); Bayly et al. (1988).

In the absence of stratification, Rayleigh’s inflection-point theorem for inviscid, bounded, parallel, homogeneous shear flow is established. This necessary but not sufficient condition states that the flow may be linearly unstable to perturbations only if $U$ has an inflection point $z_{\ast}$ Rayleigh (1880). Furthermore, a stronger condition for instability was obtained by Fjørtoft, who showed that $U^{\prime\prime}(U-U_{s})<0$ somewhere in the field flow, where $z_{\ast}$ is an inflection point at which $U^{\prime\prime}(z_{\ast})=0$ and $U_{s}=U(z_{\ast})$ Fjørtoft (1950), see also Squire (1933); Høiland (1953); Drazin (2002); Drazin and Reid (2004), cf. Gubarev (2007, 2013) for further discussion on necessary and sufficient conditions. In what follows, we will consider the special case where $A_{T}=0$ and nonzero $A_{H}$.

2.2 Modified Taylor-Goldstein equation with horizontal eddy viscosity

Let $U(z)$ be the background velocity, $A_{H}(z)$ be the horizontal eddy viscosity, where both are real-valued functions of the vertical coordinate $z$. Let $c=c_{r}+ic_{i}\in\mathbb{C}$ with $c_{i}>0$ for unstable situation. We write $U(z)-c=u_{1}(z)=\rho_{1}e^{i\phi_{1}}$ and $U(z)-c-ikA_{H}=u_{2}(z)=\rho_{2}e^{i\phi_{2}}$, where both the amplitudes $\rho_{1,2}(z)$ and the phases $\phi_{1,2}(z)$ are real-valued functions. The square of the Brunt-Väisälä frequency is defined as $N^{2}=g\frac{d}{dz}\ln\Theta(z)$, with $\Theta$ as the basic state of potential temperature. Consider a modified TG equation for the vertical velocity profile $\hat{w}$ in the absence of $A_{T}$ but only the presence of $A_{H}$. It reads

[TABLE]

where both $Q_{0}$, $Q_{1}\in\mathbb{C}$ and are given as follows:

[TABLE]

Applying the following Liouville transformation by introducing a new variable $\tilde{w}$ and an integration factor

[TABLE]

we can simplify the modified TG equation (9) by removing the first-order derivative term. It yields

[TABLE]

where

[TABLE]

Rescaling the vertical coordinate $z$ as $Z=\varepsilon z$, where $\varepsilon\ll 1$ is a small dimensionless parameter, then the modified TG equation (13) reads

[TABLE]

where the complex-valued potential term $Q_{\varepsilon}$ is given by

[TABLE]

with

[TABLE]

Note that the differentiation with respect to the scaled vertical variable $Z$ has now been replaced by a dot.

3 WKB method and matrix exponential

WKB theory provides relatively simple approximate solutions for the Fourier component of the perturbed vertical velocity $\tilde{w}$, which is expressed as a linear second-order ODE given by the modified TG equation (13) or (15). When the frequency-like term $Q_{\varepsilon}$ (16) varies slowly vertically, an approximate solution to (15) can be sought by implementing the WKB method.

WKB theory is a particular case of multiple-scale analysis and an early application of this method comes from quantum mechanics when the conditions of the particle motion are nearly classical. This method does not only enjoy applications but also has been successfully implemented in the lower troposphere and boundary layer regions, including internal gravity Gill (1982); Bretherton (1966); Pedlosky (2003), atmospheric acoustic-gravity Pitteway and Hines (1965); Einaudi and Hines (1970), radio Budden (1988), and mountain waves Laprise (1993); Grisogono (1994); Teixeira et al. (2004); Teixeira and Miranda (2004, 2006); Miranda et al. (2009); Teixeira and Miranda (2009).

3.1 WKB approximation

Consider the modified TG equation (15) and write its potential term $Q_{\varepsilon}$ by separating the real and imaginary parts, given as follows:

[TABLE]

where

[TABLE]

Let also $M(Z)=m(Z)+in(Z)$, $m,n\in\mathbb{R}$, where all three quantities $M$, $m$, and $n$ are expanded in the order of $\varepsilon$:

[TABLE]

Note that two conditions must be satisfied for the WKB approach to be valid and useful. First, the series $\varepsilon^{-1}\int M(Z)$ must be an asymptotic series in $\varepsilon$ as $\varepsilon\rightarrow 0$ uniformly for all $Z$ in the interval of approximation. Second, for the WKB solution truncated at the term $\varepsilon^{j-1}M_{j}$, the next term must be small compared to one for all $Z$ in the interval of approximation, i.e.,  $\varepsilon^{j}|M_{j+1}(Z)|\ll 1$ as $\varepsilon\rightarrow 0$ Bender and Orszag (2010).

The WKB anzats for (15) is given by

[TABLE]

Substituting this anzats to (15), collecting the terms in the order of $\varepsilon$, we obtain coupled equations by separating the real and imaginary parts. The lowest order equations are given by

[TABLE]

Solving these equations, we obtain

[TABLE]

By noting $m_{0}$ and $n_{0}$ as the positive roots of the expressions above, there are two fundamental solutions for the geometrical optics WKB approximation of (15):

[TABLE]

The corresponding solution represents a plane-wave with amplitude independent of the vertical variable $z$. This geometrical optics approximation satisfies

[TABLE]

where $A+iB+M_{0}^{2}=0$, which contributes from the lowest-order term of $Q_{\varepsilon}$. From here, we observe that for the WKB approximation to be valid, then we require that the condition $\left|M_{0}^{\prime}(z)\right|\ll\left|M_{0}^{2}(z)\right|$ needs to be satisfied. Hence, we can obtain a reasonable definition for $\varepsilon>0$ Teixeira et al. (2004):

[TABLE]

As an example, for $U(z)=1+2\tanh(z)$, $N^{2}(z)=1+4\tanh^{2}(z)$, $c=5+0.1i$, $k=0.1$, and $A_{H}(z)=2e^{-z^{2}}$, any value of $\varepsilon$ satisfying $0<\varepsilon\leq 0.3$ would be sufficient.

The first-order of $\varepsilon$ yields the following system of equations:

[TABLE]

Solving these equations and using the fact that $|M_{0}|^{2}=m_{0}^{2}+n_{0}^{2}=\sqrt{A^{2}+B^{2}}$, we obtain $m_{1}$ and $n_{1}$ and their integrals:

[TABLE]

The physical optics approximation of the WKB method reads

[TABLE]

where we have defined

[TABLE]

In the literature, what is called as the “WKB solution” usually refers to the physical optics approximation. The geometric approximation is rarely used when we discuss the WKB method, as it shortfalls to provide the correction for the spatially varying coefficients.

The system of equations for the second-order $\varepsilon$ is given by

[TABLE]

These solve as

[TABLE]

Although in principal we can proceed further to obtain an expression for the WKB approximation at higher-order of $\varepsilon$, integrating $M_{2}(Z)$ is already cumbersome at this stage due to its dependence on $C(Z)$ and $D(Z)$. Hence, the WKB approximation (18) suffices for our analysis. The left panel of Figure 1 displays a physical optics approximation of the WKB solution (18) for a Gaussian profile of horizontal eddy viscosity $A_{H}$ at the onset of instability. The horizontal background velocity and the square of the Brunt-Väisälä frequency are taken as $U(z)=1+2\tanh(z)$ and $N^{2}(z)=1+4\tanh^{2}(z)$, respectively, which are shown in the right panel of Figure 1. Both quantities are symmetric with respect to the $z$-axis.

3.2 Matrix exponential and fundamental matrix physical optics approximations

Consider the following well-known theorem on matrix exponential and fundamental matrix solutions from ODE theory Braun (1993); Boyce and DiPrima (2012); Edwards et al. (2018); Strang (2014). Let $\boldsymbol{G}$ be an $n\times n$ matrix and let $\boldsymbol{\Phi}(z)$ be a fundamental matrix for the homogeneous linear system $\boldsymbol{x}^{\prime}=\boldsymbol{G}\boldsymbol{x}$, then the unique solution of the initial (in this case, boundary) value problem (IVP or BVP)

[TABLE]

is given by

[TABLE]

Consider again the modified TG equation (13) or (15). Defining $\tilde{v}^{\prime}=\tilde{w}$, we can write it as a system of ODEs in the form of (19):

[TABLE]

where the prime and dot signs represent the derivative with respect to $z$ and $Z$, respectively and

[TABLE]

The eigenvalues of $\boldsymbol{G}_{\varepsilon}$ are complex-valued functions that depends on the variable $Z$:

[TABLE]

Figure 2 depicts the real and imaginary parts of eigenvalues $\lambda$ as well as their trajectories in the complex plane for two different values of $\varepsilon$.

We need to seek the corresponding eigenvectors from these eigenvalues to obtain exact solutions using the fundamental matrix and matrix exponential, which turns out to be nontrivial. Instead, we will use the corresponding eigenfunctions from the WKB solution that we have obtained in the previous subsection to construct the physical optics approximation for the matrix exponential. Pursuing this course is particularly useful once we follow the path in solving the BVP numerically. One notable novel numerical integration technique combines the standard Runge-Kutta approach and the WKB method, known as the Runge-Kutta-Wentzel-Kramers-Brillouin (RKWKB) method Handley et al. (2016); Agocs et al. (2020); Bamber and Handley (2020).

For the geometrical optics approximation of the WKB method, the complex-valued frequency can be approximated as

[TABLE]

The physical optics approximation provides a correction to the phase and amplitude of the solution. Using the fact that there are two linearly independent eigenfunction solutions, the matrix exponential approximation can be expressed as follows:

[TABLE]

where the Wronskian $W(z)=\text{det}\,\boldsymbol{\Phi}(z)=-2q(z)/\sqrt{M_{0}(z)}$ and

[TABLE]

We observe that $\sqrt{\left|\frac{M_{0}(0)}{M_{0}(z)}\right|}\cosh\phi(z)$ and $\sqrt{\left|\frac{M_{0}(0)}{M_{0}(z)}\right|}\sinh\phi(z)$ play a role as two linearly independent eigenfunctions of the WKB solution. This is guaranteed from the fact that if the Wronskian starts from $W\neq 0$ at $z=0$, then $W(z)\neq 0$ for all $z>0$ Strang (2014).

4 Behavior near turning point and critical level

There are at least two possibilities where the WKB approximation becomes invalid, i.e., at a turning point or at a critical level. In this section, we will discuss the behavior of the perturbed vertical velocity profile near these points.

4.1 Behavior near a turning point

A turning point is the point $z=z_{0}$ where $M_{0}$ vanishes, i.e., when it satisfies $N^{2}(z_{0})=k^{2}\left[U(z_{0})-c\right]\left[U(z_{0})-c-ikA_{H}(z_{0})\right]$. At this incidence, the WKB solution (18) becomes singular. Let $f(z)=N^{2}-k^{2}u_{1}u_{2}$, then its Taylor expansion about the turning point $z_{0}$ gives

[TABLE]

Multiplying the TG equation (13) with $u_{1}u_{2}$, taking the first-order of Taylor-expansion of $f$, and dividing again with $u_{1}u_{2}$ after evaluating them at the turning point $z_{0}$, we obtain

[TABLE]

By defining a new variable $\zeta$:

[TABLE]

we arrive at the Airy ODE

[TABLE]

The Airy equation (20) has two linearly independent solutions, called the Airy function of the first and second kinds, denoted by Ai$(\zeta)$ and Bi$(\zeta)$, respectively. Since Bi$(\zeta)$ blows up for $\zeta>0$, this will not correspond to the physical setting of the stratified shear flows. We are more interested in Ai$(\zeta)$ because it decays for $\zeta>0$ or for $z<z_{0}$. The function Ai$(\zeta)$ is defined by the following improper Riemann integral:

[TABLE]

The convergence of this function can be shown by Dirichlet’s test using the integration by parts Hall (2013).

4.2 Behavior near a critical level

A critical level (or steering level) is a height $z=z_{c}$, if it exists, at which the background state of the velocity flow $U(z)$ is equal to the real part of the horizontal phase speed $c$, i.e., $U(z_{c})=c_{r}$. When this height is reached and $c\in\mathbb{R}$ ($c_{i}=0$), the denominator of $Q_{2}$ or $Q_{\varepsilon}$ vanishes. Hence, the modified TG equation (13) or (15) becomes singular, and consequently, the WKB solution (18) vanishes too. Assume that $c_{i}$ is small, in particular, $|U^{\prime\prime}(z)|c_{i}/[U^{\prime}(z)]^{2}\ll 1$. For stable flow, the function $\tilde{w}(z)$ oscillate rapidly in $z$ near the critical level and the oscillations become infinitely rapid as the wave approaches it Haynes (2051). Following an analysis presented in Booker and Bretherton (1967); Nappo (2013), we examine the behavior of the vertical velocity profile near the critical level.

Denoting $\zeta=z-z_{c}$, we proceed with Taylor-expanding the background wind speed $U$ as well as $U+ikA_{H}$ up to second-order terms

[TABLE]

Here, $\alpha_{1}=kA_{H}^{\prime}(z_{c})$ and $\alpha_{2}=kA_{H}^{\prime\prime}(z_{c})$. Upon substitution to the modified TG equation (13), we obtain another TG equation with a regular singular point at $\zeta=0$:

[TABLE]

where

[TABLE]

We seek a series solution near the critical level by employing the method of Frobenius by expanding $\tilde{w}$ as follows Braun (1993); Boyce and DiPrima (2012); Edwards et al. (2018); Frobenius (1873):

[TABLE]

Substituting (28) to (24) and ordering the terms according to the powers of $\zeta$ yields the indicial equation $\lambda(\lambda-1)+Q_{22}=0$, with indicial roots or characteristic exponents at the singularity $\lambda$ at the regular singular point $\zeta=0$ are given by

[TABLE]

where

[TABLE]

We can express the solution of the modified TG equation (24) as the following linear combination:

[TABLE]

where $\tilde{w}_{1}$ and $\tilde{w}_{2}$ are two linearly independent Frobenius solutions, given as follows:

[TABLE]

The real-valued coefficients are given by $A_{n}=\text{Re}\left\{C_{n}\right\}$ and $B_{n}=\text{Im}\left\{C_{n}\right\}$, $n\in\mathbb{N}_{0}$, and the complex-valued coefficients $C_{n}$, $n\in\mathbb{N}_{0}$ satisfy the following recursive relationship:

[TABLE]

The behavior of Frobenius solutions is rather nontrivial. Although both trigonometric terms are bounded, since they oscillate infinitely many times between their maximum and minimum values, they do not approach any single value as $z\rightarrow z_{c}$, cf. Haynes (2051). Additionally, we also need to ensure the positive power for the small term $|z-z_{c}|$ as the limit of $|z-z_{c}|$ for negative power as $z\rightarrow z_{c}$ does not exist.

In what follows, we will verify that the corresponding power series (28) without the term $\zeta^{\lambda}$ converges for all $|\zeta|<R$, where $R\geq 0$ is the radius of convergence of the term $(Q_{22}-Q_{21}\zeta+Q_{20}\zeta^{2})$. In this particular case, $R=\infty$. Let $I(n+\lambda)=n(2\lambda+n-1)$ be a second-degree polynomial function in $n$ and $N_{0}\in\mathbb{N}$ such that $I(n+\lambda)\neq 0$ for $n=N_{0}$, $N_{0}+1$, $N_{0}+2$, $\cdots$. The coefficients of the power series (28) in the absence of $\zeta^{\lambda}$ satisfy the recursive relationship (30).

Let us express the modified TG equation (24) as $\tilde{w}^{\prime\prime}+\zeta^{-2}\widetilde{Q}(\zeta)\tilde{w}=0$, where

[TABLE]

with

[TABLE]

Note that this particular representation extends easily to the case where $\widetilde{Q}(\zeta)$ is any power series convergent for $|\zeta|<R$, where $R>0$ is the corresponding radius of convergence. For any $\zeta\in(-R,R)$, choose $\Lambda_{1}>0$ and $\Lambda_{2}>0$ such that for $k\in\mathbb{N}$:

[TABLE]

Furthermore, since $\tilde{q}_{k}=0$ for $k\geq 3$, we have

[TABLE]

We now consider a modified power series

[TABLE]

with

[TABLE]

One could verify immediately that

[TABLE]

As a consequence, the convergence of the series (28) without the factor $\zeta^{\lambda}$ can be confirmed by showing the convergence of the series (31) by means of the Ratio Test. Moreover, for $n>N_{0}$, we have

[TABLE]

Since $I$ is a second-degree polynomial in $n$, we can straightforwardly verify that

[TABLE]

Since $|\zeta|<\Lambda_{1}$, we obtain

[TABLE]

Hence, the Frobenius series solution (29) converges for $|\zeta|<R$. A similar proof for a more general second-order linear homogeneous ODEs and an exposition of different possible solutions obtained using the Frobenius method can be found in advanced texts on differential equations, e.g., Agarwal and O’Regan (2009); Birkhoff and Rota (1989); Coddington (1989); Howell (2019); Teschl (2012). Note that the Frobenius solution presented in Booker and Bretherton (1967); Nappo (2013) only includes two terms of the lowest-order and discards the higher-order terms in the series. This suggests that those higher-order terms might not be essential and can be neglected from the physical point of view.

Finally, by employing the continuity equation (6), applying the Liouville transformation (12), and restoring the normal mode, one may obtain the corresponding horizontal and vertical velocity fields.

5 Conclusion

We have considered a modified TG equation describing a model for complex-valued vertical velocity perturbation in stratified shear flows with horizontal eddy coefficient of turbulent viscosity. Under an assumption of slowly varying background horizontal velocity and Brunt-Väisälä frequency, we implemented an asymptotic approach using the WKB method to obtain the geometrical and physical optics solutions for the modified TG equation.

Furthermore, we also investigate the behavior of the WKB solution near a turning point and seek a series solution near a critical level by employing Frobenius’ method. The former can be described by the Airy function while the latter exhibits an oscillating behavior but vanishes as the height approaches a critical level.

For the future work, it would be illuminating to discuss and analyze rigorously the relationship between the WKB method and a transformation at complex geometry. Although there exists an abundant literature on complex WKB method, e.g., Voros (1983); Maslov (1994); Fedotov and Klopp (2004), most of them deals with the complex variable. In our case, however, the variable remains real-valued even though the potential-like function is complex-valued. Furthermore, we are interested in investigating the effect of vertical eddy coefficients of turbulent viscosity where its contribution is more essential in the study of atmospheric boundary layer in comparison to the horizontal one. Additionally, we would like to investigate the case when both horizontal and vertical eddy coefficients of turbulent diffusivity are present in the model equation.

\funding

“This research received no external funding.”

Acknowledgements.

The author would like to thank Kenneth Howell for illuminating on the convergence proof of the Frobenius series, Vera Liu, Andy Chan, Ardhasena Sopaheluwakan, Francisco Marcelo Fernández, Lars Pesch, Mason Porter, Ramoshweu Solomon Lebelo, Shijun Liao, Sivasankaran Sivanandam, Charles Steele for fruitful discussions, and four anonymous referees for the suggestions in improving this article, whom two of them bringing the references Khani and Waite (2013, 2020); Gubarev (2007, 2013) to our attention. \conflictsofinterest“The authors declare no conflict of interest.” \abbreviationsThe following abbreviations are used in this manuscript:

[TABLE]

\reftitleReferences