# Disk wakes in nonlinear stratification

Divyanshu Gola 1  S. Nidhan1  J. L. Ortiz-Tarin1  Hieu T. Pham1  S. Sarkar1
<EMAIL_ADDRESS>1Department of Mechanical and Aerospace Engineering,
University of California San Diego, CA 92093, USA

###### Abstract

Nonlinearity of density stratification modulates buoyancy effects. We report
results from a body-inclusive large eddy simulation (LES) of a wake in
nonlinear stratification, specifically for a circular disk at diameter-based
Reynolds number ($\Rey$) of $5000$. Five density profiles are considered; the
benchmark has linear stratification and the other four have hyperbolic tangent
profiles to model a pycnocline. The disk moves inside the central core of the
pycnocline in two of those four cases and, in the other two cases with a
shifted density profile, the disk moves partially/completely outside the
pycnocline. The maximum buoyancy frequency ($N_{max}$) for all the profiles is
the same. The first part of the study investigates the centered cases.
Nonuniform stratification results in increasing wake turbulence relative to
the benchmark owing to reduced suppression of turbulence production as well as
the wave trapping in the pycnocline. Steady lee waves are also quantified to
understand limitations of linear theory. The second part pays attention to the
effect of a relative shift between the pycnocline and the disk. The wake
defect velocity decays faster in the cases with a shift. The effect of disk
location on the Kelvin wake waves (a family of steady waves within the
pycnocline) and its modal form is obtained and explained by solving the
Taylor-Goldstein equation. The family of unsteady internal gravity waves that
are generated by the wake is also studied and the effect of disk shift is
quantified.

##  1 Introduction

A pycnocline is a layer of fluid whose density changes rapidly with depth and
so does the stratification described by the buoyancy frequency ($N$), where
$N^{2}=-(g/\rho_{0})\partial\rho/\partial z$. Such layers in the ocean or
atmosphere affect environmental turbulence and waves, the transport of tracers
(nutrients, pollutants, etc.) and the interaction of the environment with
engineered structures. The upper-ocean pycnocline is flanked by a mixed layer
on the top and a deeper layer at the bottom. The mixed layer and the deep
layer can also be stratified but their stratification levels are typically
lower than that in the pycnocline. Since much of our knowledge of turbulence
in stratified wakes is derived from canonical wakes in uniform stratification,
it is useful to study the wake of a canonical body inside and near a
pycnocline where the change in stratification is rapid and nonlinear.

Most of the past work on the flow features of bodies moving through a
pycnocline has involved the internal wave structure. Robey (1997) studied
internal waves generated by a sphere moving below a pycnocline using
experimental and numerical techniques for a wide range of body-based $\Rey$
and Fr. Nicolaou et al. (1995) experimentally studied the waves of an
accelerating sphere in a thermocline. The phase configuration for two-
dimensional trapped internal waves was theoretically and experimentally
studied by Stevenson et al. (1986) for a cylinder moving in a thermocline
generated using salt brine. Often, these experimental observations of wave
patterns are compared to the theoretical results of earlier studies, e.g.,
Barber (1993) and Keller & Munk (1970), where analytical results involving the
dispersion relation of internal waves and their modal wave forms for nonlinear
stratification profiles are provided.

A hyperbolic tangent profile bridging two regions with different values of $N$
has been often used to model nonlinear density variation within and density
jump across a pycnocline. Ermanyuk & Gavrilov (2002, 2003) used a hyperbolic
tangent profile to study forces on an oscillating cylinder and sphere, while
Grisouard et al. (2011) used a hyperbolic tangent profile with a constant
bottom $N$ to study internal solitary waves in a pycnocline using direct
numerical simulation. A hyperbolic profile of $N$ has been used for studying a
weakly stratified shear layer adjacent to a uniformly stratified region (Pham
et al., 2009) as well as an asymmetrically stratified jet (Pham & Sarkar,
2010).

There have been numerous studies on evolution of turbulent wakes under
stratification. An early study of Froude number $O(10-1000)$ wakes by Lin &
Pao (1979) showed that stratification starts to affect the wake at buoyancy
time scale ($Nt$) of $O(1)$ resulting in a non-axisymmetric evolution. The
multistage wake decay was later characterized by Spedding (1997) into three
separate regimes, namely three-dimensional (3-D), non-equilibrium (NEQ), and
quasi-two-dimensional (Q2D) and quantitative turbulence measurements were
obtained in this and following laboratory studies.

Temporal simulations have enabled longer simulations in terms of evolved time
($T$) or streamwise distance ($x$) using the transformation $T=x/U_{\infty}$,
where $U_{\infty}$ is the free-stream velocity. Gourlay et al. (2001) found
the appearance of pancake vortices in the Q2D regime in their temporal DNS of
a wake at $\Rey=10^{4}$ and $\mbox{{Fr}}=10$. Brucker & Sarkar (2010) in their
DNS of towed and self-propelled wakes at $\Rey=10^{4}$ and $5\times 10^{4}$
quantified the wake turbulence statistics and found that the higher defect
velocity in stratified wakes is a result of the buoyancy induced reduction of
the Reynolds shear stress and ultimately the turbulent production term. This
results in a longer lifetime of the stratified wake, also verified later by
Redford et al. (2015) in their DNS of a weakly stratified turbulent wake. They
also found that the horizontal nature of the wake in Q2D regime resulted in
the dominance of the lateral Reynolds shear stress so that the decay of the
mean wake velocity was faster (compared to NEQ regime) and was accompanied by
an increase in the turbulent kinetic energy. Other studies that use temporal
wake models include Dommermuth et al. (2002), Diamessis et al. (2005), de
Stadler & Sarkar (2012) and Abdilghanie & Diamessis (2013).

Temporal simulations have been shown to be sensitive to the choice of initial
conditions (Redford et al., 2012). Besides, inclusion of vortex shedding at
the body and lee wave generation presents a challenge to temporal models. In
recent work on stratified wakes, these limitations have been overcome by body-
inclusive simulations. Orr et al. (2015) conducted a numerical study of a
sphere wake at $\Rey=200$ and $1000$ where they identified vortex shedding and
lee waves. Pal et al. (2016) simulated the wake of a sphere at $\Rey=3700$
using DNS and found that there is a resurgence of turbulent fluctuations below
a critical Fr instead of a monotone suppression. This observation was further
examined by Chongsiripinyo et al. (2017) who found that, despite the low Fr,
vortex stretching was dominant in the near wake, resulting in small-scale,
three dimensional turbulence. The decay of a disk wake at $\Rey=5\times
10^{4}$ was examined by Chongsiripinyo & Sarkar (2020) who decomposed its
evolution into three regimes based on buoyancy time scale and horizontal
Froude number of the fluctuations (instead of the mean wake): weakly
stratified, intermediately stratified and strongly stratified turbulence.
Turbulence scales were also used to characterize wake transitions by Zhou &
Diamessis (2019). Ortiz-Tarin et al. (2019) performed LES of flow past a
prolate spheroid with laminar boundary layer separation into the near and
intermediate wake and found that at $\mbox{{Fr}}\sim O(1)$, buoyancy effects
are stronger for slender bodies as compared to bluff bodies.

Interaction of the turbulent wake with the stratified background leads to the
generation of unsteady internal waves. Abdilghanie & Diamessis (2013) found
that the NEQ regime is prolonged at high $\Rey$, resulting in internal wave
radiation that persists up to $Nt\approx 100$. Brucker & Sarkar (2010) showed
that internal wave flux dominates the turbulent dissipation during $20<Nt<75$
for a wake at $\Rey=5\times 10^{4}$ and $\mbox{{Fr}}=4$. In contrast, Redford
et al. (2015) showed that for very high Fr, the internal wave activity,
although still more pronounced at the start of the NEQ regime, makes a
negligible contribution to the turbulent energy budget. Rowe et al. (2020)
characterized the angles for the strongest internal waves and found that after
$Nt=10$, internal wave radiation is an important sink for wake kinetic energy.
Nidhan et al. (2022) used spectral proper orthogonal decomposition and
connected the leading wake eigenmodes to the wave field. Besides numerical
techniques, many experimental (Gilreath & Brandt, 1985; Hopfinger et al.,
1991; Chomaz et al., 1993; Bonneton et al., 1993; Brandt & Rottier, 2015;
Meunier et al., 2018) as well as theoretical (Sturova, 1974; Voisin, 1991,
1994, 2003, 2007) studies have been done to analyze the internal wave field in
a stratified flow.

Study of the wake in the vicinity of a pycnocline with a characterization of
turbulence and its interaction with trapped internal waves is limited. Pham &
Sarkar (2010) studied the interaction between internal waves from a shear
layer and an adjacent stratified jet, and classified waves that are
trapped/reflected by the jet and waves that transmit through the jet.
Sutherland & Yewchuk (2004) and Sutherland (2016) calculated analytical
expressions for internal wave transmission through density staircases for a
stationary, two-dimensional Boussinesq fluid. Voropayeva & Chernykh (1997)
simulated a temporally evolving wake in a nonlinear stratification using a
Reynolds-averaged turbulence model. To the best of the authors’ knowledge, the
present study is the first body-inclusive wake LES for a nonlinearly
stratified fluid. We aim to assess the effect of nonlinear stratification on
wake turbulence, its interaction with trapped waves as well as characterize
the far-field lee waves. The disk is chosen as a canonical generator of a
bluff-body wake to avoid the computational expense incurred in resolving the
boundary layer that develops on a long body.

Section 2 describes the numerical setup and the simulation parameters. The
results are divided into two sections. The basic differences between linear
and nonlinear stratification when the disk is centered in the pycnocline layer
are studied in section 3. The effect of relative shift between the density
profile and the disk, so that the disk is partially/completely outside the
pycnocline layer, is discussed in section 4. We summarise and conclude the
study in section 5.

##  2 Methodology

###  2.1 Governing equations and numerical scheme

The wake of a disk is simulated by solving the three dimensional,
incompressible, unsteady form of the conservation equations for mass, momentum
and density. A high-resolution large eddy simulation (LES) with the Boussinesq
approximation for density effects is used. The disk is immersed perpendicular
to a flow with velocity $U_{\infty}$. The equations are numerically solved in
cylindrical coordinates but both Cartesian $(x,y,z)$ and cylindrical
$(r,\theta,x)$ coordinates are appropriately used in the discussion. Here, $x$
is streamwise, $y$ is spanwise and positive along $\theta=0^{\circ}$, and $z$
is vertical and positive along $\theta=90^{\circ}$ (Fig. 1). The density field
$(\rho(\mathbf{x},t))$ is split into a constant reference density
$(\rho_{0})$, the variation of the background $(\Delta\rho_{b}(z))$, and the
flow induced deviation, $(\rho_{d}(\mathbf{x},t))$ so that
$\rho(\mathbf{x},t)=\rho_{0}+\Delta\rho_{b}(z)+\rho_{d}(\mathbf{x},t)$. The
Reynolds number of the flow, defined as $\Rey=U_{\infty}D/\nu$ with $D$ being
the diameter of the disk, is $5000$. Since the stratification is non-uniform,
the Froude number defined by $Fr(z)=U_{\infty}/N(z)D$ is variable and its
case-dependent behavior is described in section 2.2. The minimum value of the
Froude number ($Fr_{min}$) is set to $1$ for all the cases.

The filtered nondimensional equations are as follows

Figure 1: Schematic showing the simulation setup and domain.

$\frac{\partial u_{i}}{\partial x_{i}}=0,$ (1)

$\frac{\partial u_{i}}{\partial t}+u_{j}\frac{\partial u_{i}}{\partial
x_{j}}=-\frac{\partial p}{\partial
x_{i}}+\frac{1}{\Rey}\frac{\partial}{\partial
x_{j}}\Big{[}(1+\frac{\nu_{sgs}}{\nu})\frac{\partial u_{i}}{\partial
x_{j}}\Big{]}-\frac{\rho_{d}}{Fr_{min}^{2}}\delta_{i3},$ (2)

$\frac{\partial\rho}{\partial t}+u_{j}\frac{\partial\rho}{\partial
x_{j}}=\frac{1}{\Rey Pr}\frac{\partial}{\partial
x_{j}}\Big{[}(1+\frac{\kappa_{sgs}}{\kappa})\frac{\partial\rho}{\partial
x_{j}}\Big{]},$ (3)

where $u_{i}$ refers to the filtered velocities in $x$, $y$, and $z$
directions for $i=1,2$, and $3$ respectively. $\nu_{sgs}$ and $\nu$ in Eq. (2)
are the subgrid-scale kinematic viscosity and the kinematic viscosity,
respectively, while $\kappa_{sgs}$ and $\kappa$ in Eq. (3) are the subgrid
scale diffusion coefficient and the diffusion coefficient, respectively. The
Prandtl number defined by $Pr=\nu/\kappa$ is set to one. No major qualitative
influence of $Pr$ was found in the wake study by de Stadler et al. (2010) who
varied $Pr$ between 0.2 and 7. The dynamic eddy viscosity model (Germano et
al., 1991) is employed to obtain $\nu_{sgs}$, following the implementation of
Chongsiripinyo & Sarkar (2020), and the subgrid Prandtl number is also set to
unity to obtain $\kappa_{sgs}$. The parameters used for nondimensionalization
are as follows: free-stream velocity ($U_{\infty}$) for velocity, disk
diameter ($D$) for length, advection time ($D/U_{\infty}$) for time, dynamic
pressure ($\rho U_{\infty}^{2}$) for pressure, and characteristic change in
background density ($-D(\partial\Delta\rho_{b}/\partial z)|_{max}$) for the
flow induced density deviation.

The periodicity in the azimuthal direction is leveraged in solving the
discretised pressure Poisson equation in the predictor step, reducing it to a
pentadiagonal system of linear equations, which is then solved using a direct
solver (Rossi & Toivanen, 1999). The disk is represented by the immersed body
method of Balaras (2004) and Yang & Balaras (2006). At the inlet boundary, a
uniform stream of velocity ($U_{\infty}$) is imposed while an Orlanski-type
convective boundary condition is used for the outflow (Orlanski, 1976).
Neumann boundary condition is imposed at the radial boundary of the domain for
density as well as the three velocity components. To prevent spurious
reflection of waves back into the domain, sponge layers are used at the inlet
(streamwise length $10D$), outlet (streamwise length $10D$) and the
cylindrical walls (radial length $15D$) of the domain to gradually relax the
velocities and the density to their respective background values at the
boundaries.

###  2.2 Density profiles

The profiles chosen for the variation of background density in the five cases
are shown in Fig. 2. Except for the benchmark 111 case, the profiles have
nonuniform $N(z)$ and $Fr(z)=U/N(z)D$. Profile 111 with linear stratification
has a constant $Fr(z)$ which is the conventional body-based Froude number of
$\mbox{{Fr}}=U_{\infty}/ND=1$. The maximum value of the density gradient is
the same for all four nonlinear profiles and corresponds to $Fr_{min}=1$. In
Profile 614 (Fig. 2b), a hyperbolic tangent function is used to bridge two
linearly stratified regions (6 and 4 represent the value, rounded to an
integer, of the local $Fr$ in the linearly stratified regions). Profile-I1I (I
standing for “Infinity”) is a complete hyperbolic tangent profile between two
constant-density regions, therefore having infinite Froude number far above
and below the body.

For the I1I and 614 profiles, the disk center coincides with the center of the
density profile and, furthermore, the disk lies entirely within the
nonlinearly stratified region. Profiles $614-$ and $614+$ are obtained after
shifting profile 614 (and not the disk) vertically by $-2.5D$ and $3D$,
respectively. The disk center is located at $z=0$ for all five profiles and
the negative and positive signs in the names $614-$ and $614+$ indicate
whether the profile is shifted downwards or upwards, respectively. The
vertical shift in $614-$ (Fig. 2d) is chosen so that the upper half of the
disk lies in the linearly stratified region with $\mbox{{Fr}}=6$ and the lower
half in the pycnocline. The vertical shift in $614+$ (Fig. 2e) is chosen so
that the disk lies entirely in the linearly stratified region with
$\mbox{{Fr}}=4$ while being very close to the pycnocline, specifically, its
upper edge is $0.5D$ below the pycnocline.

The profile I1I is given by the following function, which was also used by
Ermanyuk & Gavrilov (2002, 2003) and Nicolaou et al. (1995) in their studies:

$\Delta\rho_{b}(z)=-\frac{\rho_{0}\varepsilon}{2}\tanh{\frac{2z}{\delta}}.$
(4)

The central nonlinear region of the 614 profile is obtained from the I1I
profile. The linearly stratified regions of profile 614 are added by matching
the slopes at $z=1.25\delta$ and $z=-\delta$ so that the density profile is
continuously differentiable in $z$. This results in Froude numbers of $6.13$
and $3.74$ above and below the pycnocline layer respectively, for profile 614.

Figure 2: Variation of background density for (a) 111, (b) 614, (c) I1I, (d)
$614-$, (e) $614+$, along with (f) respective buoyancy frequency profiles.
Solid black line at the center of (a-e) represents the disk (not to scale).
Profiles are summarized in table 1 and also discussed in text. For $614-$, the
upper half of the disk is in a constant-$N$ region and the lower is in the
pycnocline. For $614+$, the disk is in the bottom constant-$N$ region with its
upper edge $0.5D$ below the pycnocline.

The profiles $614-$ and $614+$ are shifted with respect to 614 and their
central nonlinear region is as follows:

$\Delta\rho_{b}(z)=-\frac{\rho_{0}\varepsilon}{2}\tanh{\frac{2(z-s)}{\delta}},$
(5)

where $s$ denotes the shift and takes values $-2.5D$ and $3D$ for $614-$ and
$614+$, respectively.

The local Froude number corresponding to this pycnocline region can be
calculated using the local background buoyancy frequency:

$Fr(z)=\frac{U_{\infty}}{\sqrt{\varepsilon
gD}}\sqrt{\frac{\delta}{D}}\;\Big{|}\cosh{\frac{2(z-s)}{\delta}}\Big{|}.$ (6)

The minimum value of $Fr(z)$ occurs at $z=s$. Thus, the minimum Froude number
for profiles 614 and I1I is

$Fr_{min}=Fr(z=s=0)=\frac{U_{\infty}}{\sqrt{\varepsilon
gD}}\sqrt{\frac{\delta}{D}}.$ (7)

Eq. 7 highlights two important non-dimensional parameters for a body moving
through a pycnocline: $U_{\infty}/\sqrt{\varepsilon gD}$, which is the
conventional Froude number defined using reduced gravity $\varepsilon g$, and
$\delta/D$, which is the nondimensional thickness of the pycnocline. In the
simulations, $U_{\infty}/\sqrt{\varepsilon gD}=1/\sqrt{2}$ and $\delta/D=2$
for profiles 614 and I1I to obtain $Fr_{min}=1$. Note that for 111, which has
constant linear stratification throughout, $Fr_{min}=1$ still holds.
Parameters for each case are given in Table 1.

###  2.3 Domain and Grid

The domain extends from $x/D=-L_{x}^{-}=-30$ to $x/D=L_{x}^{+}=102$ in the
streamwise direction and from $r/D=0$ to $r/D=L_{r}=60$ in the radial
direction. The number of grid points employed for discretizing the domain is
$N_{x}=2176$, $N_{\theta}=128$, and $N_{r}=479$ in the streamwise, azimuthal
and radial directions respectively, resulting in approximately $130$ million
grid points. The LES grid is non uniform in the streamwise and radial
direction and is designed to have high resolution. In terms of the Kolmogorov
length ($\eta=(\nu^{3}/\epsilon_{k}^{T})^{1/4}$), the maximum value of $\Delta
x/\eta$ is $4.39$ and that of $\Delta r/\eta$ is $5.03$. Also, the turbulent
dissipation rate used for the calculation of $\eta$ includes the resolved
scale dissipation ($\epsilon=2\nu\langle
s_{ij}^{\prime}s_{ij}^{\prime}\rangle$) as well as the subgrid scale
dissipation
($\epsilon_{sgs}=-\langle\tau_{sgs_{ij}}^{\prime}s_{ij}^{\prime}\rangle$ with
$\tau_{sgs_{ij}}=-2\nu_{sgs}S_{ij}$), where $s_{ij}^{\prime}$ and $S_{ij}$ are
the fluctuating and mean strain rates respectively. Most of the turbulent
dissipation resides in the resolved scales.

Case | $\frac{\delta}{D}$ | $\mbox{{Fr}}_{min}$ | $\mbox{{Fr}}_{z\rightarrow+\infty}$ | $\mbox{{Fr}}_{z\rightarrow-\infty}$ | $\frac{s}{D}$ | Comment
---|---|---|---|---|---|---
$111$ | N/A | $1$ | $1$ | $1$ | $0$ |  $N(z)=$ constant (Fig. 1a)
I$1$I | $2$ | $1$ | $\infty$ | $\infty$ | $0$ |  $N(z\rightarrow\pm\infty)=0$ (Fig. 1c)
$614$ | $2$ | $1$ | $6.13$ | $3.74$ | $0$ |  $N(z\rightarrow+\infty)=0.163$, $N(z\rightarrow-\infty)=0.267$ (Fig. 1b)
$614-$ | $2$ | $1$ | $6.13$ | $3.74$ | $-2.5$ |  $N(z)$ of $614$ shifted down (Fig. 1d)
$614+$ | $2$ | $1$ | $6.13$ | $3.74$ | $3$ |  $N(z)$ of $614$ shifted up (Fig. 1e)
Table 1: Physical parameters of the simulated cases. For each case,
$\Rey=5000$ and $Pr=1$. The domain with $L_{r}=60$, $L_{\theta}=2\pi$,
$L_{x}^{-}=30$, and $L_{x+}=102$ is discretised using $N_{r}=479$,
$N_{\theta}=128$, and $N_{x}=2176$ points.

##  3 Linear vs. nonlinear stratification

###  3.1 Steady lee waves

Disturbances in stratified environments generate internal waves. These can be
steady lee waves generated by the body or unsteady internal waves generated by
turbulence in the wake of body. In this section, the characteristics of the
body generated lee waves on the vertical center plane for 111 and 614 will be
analyzed using linear asymptotic theory. Note that since I1I is essentially
unstratified away from the pycnocline layer, it does not show any lee waves.

Instantaneous contours of vertical velocity ($w$) for 111 (top half) and 614
(top and bottom halves) on a vertical plane passing through the centerline are
plotted in Fig. 3(a-c). The amplitude of steady lee waves decays moving away
from the body. The lee waves observed in 111 are symmetric with respect to
$z=0$ while 614 has two different sets of lee waves (top and bottom) owing to
the two different local values of $Fr(z)$ ($6.13$ above and $3.74$ below) in
the linearly stratified regions above and below the pycnocline layer. The
waves in the 614 case have a much smaller amplitude than in the 111 case and a
larger wavelength.

Linear theory, given by Voisin (1991, 2003, 2007) is used to analytically
predict the wavelength and amplitude of the lee waves. The theory computes a
wave function, $\chi$ from a linearized set of inviscid equations involving a
source term $q$ on the right hand side of the continuity equation to model the
moving body using a Green’s function approach for large times ($Nt>>1$), which
in our case translates to $x/D>>Fr$. Once $\chi$ is known, $w$ is calculated.
The expression for $w$ for the case of a Rankine ovoid as derived by Ortiz-
Tarin et al. (2019) and employed for a 4:1 spheroid wake can be reduced to
give the wave pattern on the vertical centerline plane as:

$w(x,y=0,z)\sim-\frac{mN}{\pi
U_{\infty}r_{xz}}\cos{\psi}\sin{\Big{(}\frac{Na}{U_{\infty}}\cos{\psi}\Big{)}}\sin{\Big{(}\frac{N}{U_{\infty}}r_{xz}\Big{)}},$
(8)

where $r_{xz}=\sqrt{x^{2}+z^{2}}$ and $\psi=\arctan(z/x)$. Also, $m$ and $a$
are calculated from the potential flow stagnation point solution for a Rankine
ovoid of length $L_{ro}$ and cross-sectional diameter $D_{ro}$:

$(L_{ro}^{2}-4a^{2})^{2}=\Big{(}\frac{8ma}{\pi U_{\infty}}\Big{)}L_{ro},\qquad
D_{ro}^{2}=\Big{(}\frac{8ma}{\pi
U_{\infty}}\Big{)}\frac{1}{\sqrt{4a^{2}+D_{ro}^{2}}}$ (9)

Note that Ortiz-Tarin et al. (2019) calculated $m$ and $a$ using the potential
flow solution of a Rankine oval, which is a 2-dimensional solution. Eq. (9) is
based on the solution for a Rankine ovoid (which is 3-dimensional), and serves
as a correction to Ortiz-Tarin et al. (2019). Note that, as per Eq. 8, the
wave amplitude decreases with decreasing $N$ or increasing Fr.

The justification for using the expression for a Rankine ovoid instead of a
disk is as follows. To apply linear theory to a bluff body such as a disk, the
separation bubble created by the disk should also be included as a part of the
extended body involved in wave generation. The complete extended body for 111,
I1I and 614 cases resembles an ovoid with $L_{ro}\approx 2.5D$ and
$D_{ro}\approx 1.5D$. It should be noted that the exact body shape is
insignificant as far as the wavelength of the waves is considered, but the
wave amplitude depends on the body shape.

Figure 3: Vertical velocity contours showing lee waves and their amplitudes
along the dashed arrow labeled by $d/D$. (a,d) - 111 ($90^{\circ}$ plane);
(b,e) - 614 ($90^{\circ}$ plane); (c,f) - 614 ($270^{\circ}$ plane).

Fig. 3(d) contrasts the wave amplitude between simulation and the theoretical
estimate of Eq. 8. The comparison is along the dashed line with an arrow in
Fig. 3(a), which is at an angle of $45^{\circ}$ to the $x$-axis. The theory,
when used with the Rankine ovoid approximation for the disk and the separated
flow behind it, is able to accurately capture the amplitude variation
specially moving away from the body for the 111 case which has linear
stratification with constant $N$.

A similar analysis is performed for the lee waves of 614 after taking the
Froude number in the analysis to be that of the linearly stratified region,
i.e. $Fr=6.13$ for Fig. 3(e) and $Fr=3.74$ for Fig. 3(f). This procedure
results in the correct prediction of the wavelength of the propagating lee
waves. Theory is able to capture the order of magnitude of the significantly
reduced (relative to 111) wave amplitude but the theoretical estimate of the
wave amplitude is an under-prediction. The disk moves through a local
stratification corresponding to $\mbox{{Fr}}=1$ but the lee waves form and
propagate in the linearly stratified regions where the Froude number is larger
($\mbox{{Fr}}=6.13$ and $\mbox{{Fr}}=3.74$). Therefore, the true wave
amplitude corresponds to Fr somewhat lower than that used in the theory and,
thus, larger than the theoretical estimate. Since the wavelength is comparable
to the variability scale of $N(z)$, WKB theory cannot be used and we do not
proceed further with linear analysis.

###  3.2 Mean defect velocity

All cases show a strong effect of buoyancy on the wake since $\mbox{{Fr}}=1$
at the disk is common to them but there are differences among the cases as
elaborated below. Fig. 4(a) shows the mean centerline defect velocity for the
three cases plotted against the streamwise distance. The 111 case shows strong
oscillatory modulation of the wake and its defect velocity, similar to the
sphere wake (Pal et al., 2017). In contrast to 111, the modulation in 614 and
I1I, although present, is small. The wavelength of modulation for 111 is $2\pi
DFr$ as in the linear-theory result (Eq. 3.1), and as noted in previous work.

All cases show a similar decay law of $U_{0}\propto x^{-0.18}$ in the NEQ
regime, which was also observed by Chongsiripinyo & Sarkar (2020) for a disk
at a higher Reynolds number of $5\times 10^{4}$. The decay law transitions to
$x^{-0.45}$ around $x/D=50(Nt=50)$ signaling the beginning of the Q2D regime.
A notable difference in the defect velocity profiles occurs at $x/D\sim 4$,
where the initial dip in the profile for 111 is larger in magnitude and occurs
earlier than that of 614 and I1I, owing to the stronger lee wave in 111 as
compared to 614 and I1I. The mean velocity contours on a vertical streamwise
cut (Fig. 4b), show that the lee wave modulation in the wake for 111 is
stronger and its separation bubble is smaller, which is consistent with the
earlier dip in the profile of defect velocity (Fig. 4a).

Figure 4: (a) Mean defect velocity at centerline. (b) Mean streamwise velocity
contours on vertical plane passing through centerline. Dashed yellow line
represents the separation bubble ($U=0$).

###  3.3 Turbulent kinetic energy

The nonlinearly stratified cases exhibit stronger levels of turbulent kinetic
energy (TKE). This difference is illustrated by the contours of TKE at
$x/D=20$ in Fig. 5(a), plotted for the three stratification profiles. The
pycnocline cases have a higher TKE (by 50 - 100 %) relative to the 111 case
over the entire wake length as shown by the area-averaged values of TKE
plotted in Fig. 5(b). The area averaging at any streamwise location is done in
the region covered by a circle of radius $3D$, therefore containing most of
the turbulent zone. The TKE is substantially larger in I1I and 614 relative to
111, especially in the NEQ regime. For example, I1I has almost twice the TKE
of 111 at $x/D=45$. It is worth noting the slight increase in TKE for the case
111 at $x/D\approx 50$ which is where the decay law for the mean defect
velocity changes (Section 3.2), again pointing towards the transition to Q2D
regime.

To explain the relatively high turbulence, we inspect the terms in the TKE
budget, especially the lateral production (the production by vertical Reynolds
shear stress is small) and the wave flux. Although the hyperbolic tangent
profile is close to a linear profile near the centerline, the weaker
stratification regions above and below are not as effective in suppressing
near-wake turbulence stresses, which results in higher TKE production further
downstream. Also, the weaker far-field stratification in 614 and I1I relative
to 111 does not allow the full frequency range of waves generated in the
$Fr=1$ central region to escape into the far field. The reduction of wave
energy flux traps fluctuations in the wake and results in a TKE increase. The
two quantities, lateral TKE production and wave flux, are diagnosed in the
following subsections.

=

Figure 5: (a) Turbulent kinetic energy contours at $x/D=20$. (b) Area averaged
values of TKE as a function of streamwise distance.

####  3.3.1 TKE production

Contours of lateral production $P_{xy}$, which is the dominant component of
TKE production for all three cases, are plotted in Fig. 6(a) at $x/D=20$. The
increased $P_{xy}$ in 614 and I1I is partly a result of the higher lateral
Reynolds shear stress $\langle-u^{\prime}v^{\prime}\rangle$.

Figure 6: (a) Lateral production contours at $\frac{x}{D}=20$. (b) Area-
averaged values of lateral production as a function of streamwise distance.

Area averaged-values of the lateral production over a circle of radius $3D$
are shown in Fig. 6(b). The lateral production is consistently larger in the
614 and I1I cases over $10<x/D<50$ and is as much as twice at some locations.

####  3.3.2 Energy flux of wake generated unsteady waves

The dispersion relation for internal gravity waves in a stratified non-
rotating medium,

$\Omega=N\cos{\Theta}\,,$ (10)

where $\Theta$ is the angle of the phase line with the vertical, limits the
maximum allowable frequency to $N$.

Figure 7: (a) Radial waveflux contours at $\frac{x}{D}=20$. (b) Total waveflux
integrated over rectangular perimeter as a function of streamwise distance.

Since weaker stratification supports a narrower frequency band of wave
propagation, 614 and I1I can be expected to allow a smaller wave energy flux
to radiate through the pycnocline layer as compared to 111. This is indeed the
case as can be seen by the contours of the radial wave flux $\langle
p^{\prime}u_{r}^{\prime}\rangle$ in Fig. 7(a), where it is evident that the
wave flux is restricted within the pycnocline layer ($-2<z/D<2$) instead of
being radiating away as in 111. Note that the mean of the radial velocity and
pressure are subtracted when computing the flux so as to discard the component
from the steady lee waves.

A quantitative comparison among the cases is performed by computing the line
integral of the internal wave flux over a rectangular perimeter. Specifically,
the integral

$\oint_{C}\frac{1}{\rho_{0}}\langle{\bf u}^{\prime}p^{\prime}\rangle.{\bf
n}\,dl$ (11)

which is a function of $x/D$, is calculated using a rectangle $C$ with
vertical sides at $y/D=\pm 4$ and top/bottom sides at $z/D=\pm 2.5$ (thereby
containing the pycnocline layer). Fig.7(b) shows that the flux of wave energy
radiated away from the wake of 111 is much larger than 614 and I1I in the
early part of the NEQ regime. However, at $x/D=40$ in the late NEQ stage, the
fluxes for 614 and I1I overtake the flux of 111. Both of these trends can be
explained by the schematic shown in Fig. 8(a). The wake generated internal
waves that are trapped inside undergo complex wave-wave interaction after
getting vertically restricted, resulting in sideward escape of the wave flux.
The sideward escape can be verified by using an appropriate scaling argument
for the vertical as well as lateral wave flux at the top and side boundaries
of the rectangle, respectively.

$\langle p^{\prime}w^{\prime}\rangle\sim p_{rms}w_{rms},\langle
p^{\prime}v^{\prime}\rangle\sim p_{rms}v_{rms}$ (12)

The values of $w_{rms}$ at the point $P_{top}$, $v_{rms}$ at the point
$P_{side}$ and $p_{rms}$ at $P_{top}$ and $P_{side}$ are shown in Fig. 8(b-e),
respectively. Comparison of Fig. 8(d) and 8(e) shows that, although pressure
fluctuations at the top boundary are highest initially for 111 among all cases
and boundaries, the pressure fluctuations for 614 and I1I increase with time
at the side boundaries because of the wave trapping and the resulting sideward
escape. The velocity fluctuations at the boundaries also show similar trends.
The sideward escape of the wave energy takes some time to manifest which, in
the simulation frame, translates to streamwise distance, hence causing the
slight delay in the increase of integrated wave flux for 614 and I1I. Fig.
8(b-e) are consistent with the following result: most of the wave flux for 614
and I1I that escapes the wake core does so in the form of lateral wave flux
$\langle p^{\prime}v^{\prime}\rangle$ unlike 111 where the vertical wave flux
$\langle p^{\prime}w^{\prime}\rangle$ is the major contributor.

Figure 8: (a) Schematic showing trapping of internal waves for 614 and I1I
compared to 111 (b) $w_{rms}$ at $P_{top}$ (c) $v_{rms}$ at $P_{side}$ (d)
$p_{rms}$ at $P_{top}$ (e) $p_{rms}$ at $P_{side}$.

##  4 Centered vs. shifted nonlinear stratifications

###  4.1 Taylor-Goldstein Equation and Kelvin wake waves

In cases involving pycnoclines, there is a family of steady waves which
resemble a Kelvin ship wave pattern, on horizontal planes inside the
pycnocline layer. Intuitively, a stratification profile with large central
value of buoyancy frequency ($\mbox{{Fr}}\leq O(1)$) that changes rapidly in
the vertical (over $\delta/D\leq O(1)$) to a small value at the edges can be
thought of as a smoothed density jump and thus one can expect a pattern
resembling that generated on the air-water interface by a moving ship. Both
centered and shifted stratification profiles show Kelvin wake wave patterns.
However, as demonstrated here, the wake wave structure differs qualitatively
from the centered 614 wake when the stratification profile is shifted with
respect to the disk center as in $614-$ and $614+$ .

The dispersion relation of air-water surface gravity waves leads to the
waveforms of air-water Kelvin ship waves. To obtain an appropriate dispersion
relation for the nonlinear continuously stratified profiles, the Taylor-
Goldstein Equation is numerically solved:

$\frac{d^{2}\phi(z)}{dz^{2}}+k^{2}\Big{(}\frac{N^{2}(z)}{\Omega^{2}}-1\Big{)}\phi(z)=0.$
(13)

Here, $\phi(z)$ is the vertical displacement eigenfunction, and $k$ and
$\Omega(k)$ are the wavenumber and angular frequency corresponding to the
dispersion relation. The TGE is numerically solved as an eigenvalue problem
for $\phi(z)$ by treating $\Omega$ as an eigenvalue and $k$ as a parameter,
eg. by Robey (1997). The boundary conditions used are

$\phi(z=-L_{r})=\phi(z=L_{r})=0.$ (14)

Also, based on Sturm-Liouville theory, the eigenfunctions are normalized
using:

$\int_{z=-L_{r}}^{z=L_{r}}\phi_{m}N^{2}(z)\phi_{n}\,dz=\delta_{mn}.$ (15)

Fig. 9(a) and 9(b) respectively show the mode 1 and mode 2 eigenfunctions for
different values of $k$. The numerically calculated dispersion relation is
also compared in Fig. 9(c) to the approximation given by Barber (1993):

$\Omega(k)=\frac{c_{p_{0}}k}{1+\frac{c_{p_{0}}k}{N_{max}}}.$ (16)

Here, $c_{p_{0}}$ is the limiting long-wave phase speed of a given mode at
$k=0$, which is obtained by extrapolation, i.e., $c_{p_{0}}$ is the slope at
the origin of the mode-specific curve in Fig. 9(c). For the chosen profiles,
$c_{p_{0}}/U_{\infty}$ for mode 1 and mode 2 are $2.23$ and $0.68$,
respectively. Since Eq. 16 provides a very good approximation to the
dispersion relation, it will be used for the remaining analysis of this
section.

Figure 9: (a) Mode-1 and (b) Mode-2 eigenfunctions for different values of
wavenumber in Eq. 13. (c) Numerically calculated dispersion relation compared
with approximation given by Barber (1993).(d) Mode-1 and (e) Mode-2 waveforms
as calculated using Eq. 17.

Using the dispersion relation, phase and group velocities associated with each
mode can be calculated using $c_{p}=\Omega/k$ and $c_{g}=d\Omega/dk$, and then
substituted in the expressions given by Keller & Munk (1970) to obtain the
modal wave patterns corresponding to each eigenfunction:

$x=\frac{\varphi
U_{\infty}\Big{(}1-\frac{c_{p}c_{g}}{U_{\infty}^{2}}\Big{)}}{k(c_{p}-c_{g})},y=\frac{\varphi
c_{g}\Big{(}1-\frac{c_{p}^{2}}{U_{\infty}^{2}}\Big{)}^{\frac{1}{2}}}{k(c_{p}-c_{g})},$
(17)

where $(x,y)$ corresponds to a locus of points on the horizontal plane
parametrically given in terms of $k$ for a constant value of phase, $\varphi$.
Fig. 9(d) and 9(e) respectively show mode 1 and mode 2 wave patterns,
respectively, plotted using Eq. 17.

Fig. 10 shows the instantaneous radial velocity contours for 111, 614, $614-$
and $614+$ plotted on a half-horizontal plane passing through the respective
pycnocline center. (For 111, the plot is on $\theta=0^{\circ}$ plane). The
waves in 614 have constant-phase lines which resemble mode 2 (shown earlier in
Fig. 9e) but not mode 1. In contrast, the phase lines for $614-$ and $614+$
are more complex. In the region $0<z<10$ they resemble the mode 1 pattern of
Fig. 9(d) but, for $z>10$, they resemble the mode 2 pattern. The manifestation
of any mode by a disturbance in the pycnocline layer depends on the location
of the disturbance with respect to the pycnocline layer. In 614, the disk
center travels along the center of the pycnocline layer to displace fluid
symmetrically in the upward and downward direction, corresponding to the
antisymmetric mode 2 eigenfunction for displacement and therefore generates
mode 2 waves. Any vertical offset from the pycnocline center modifies the
antisymmetric pattern of the displacement so as to also involve some
contribution from the mode 1 waveform. The presence of a mode 2 wave pattern
in the central region of the lateral plane and a mode 1 pattern otherwise is
in agreement with the visualizations in the experiment of Robey (1997). For
profile $614+$, the contours look similar to that of $614-$ but with slightly
lower intensity because of the disk being further away from the center of the
pycnocline. (Note that the wavepattern for I1I is the same as that of 614
because they have the same hyperbolic tangent function modeling the
nonlinearity.)

Figure 10: Instantaneous contours of radial velocity on a half-pycnocline-
center plane: (a) 111, (b) $614-$, (c) 614, and (d) $614+$. Dashed lines
correspond to mode 2 waves which were shown in Fig. 9(e).

###  4.2 Distinction between lee waves and Kelvin wake waves

Figure 11: Instantaneous contours of radial velocity for the centered cases on
vertically offset horizontal planes for 111 (left column), 614 (middle) and
I1I (right). The vertical offset of each plane increases from the top to
bottom row.

Kelvin wake waves are steady in the frame of reference of the moving body like
the lee waves but constitute a distinct family of waves and appear only when
the stratification is nonuniform. Fig. 11 shows instantaneous contours of
radial velocity on horizontal planes for the centered profiles 111, 614, and
I1I at different vertical locations. At $z/D=-2.5$ (Fig. 11 (a-c)), the
contours have superficial resemblance but they represent fundamentally
different waves. For 111, since there is no nonlinear gradient, Kelvin wake
waves are not formed on the horizontal plane passing through the domain
centerline at $z/D=0$, (e.g. the previously shown Fig. 10a), and so Fig. 11(a)
has the horizontal imprint of only lee waves. However, for 614 and I1I, the
contours at $z/D=-2.5$ in Fig. 11(b,c) are mostly the imprints of the Kelvin
wake waves in Fig. 10(c). For 614, Kelvin wake waves disappear farther away
from the region of nonlinear stratification and the contours transition to
those of lee waves corresponding to $\mbox{{Fr}}=3.74$ (Fig. 11(e,h,k)). For
I1I, once the Kelvin wake waves disappear, lee waves are not observed since
the profile is unstratified outside the pycnocline region (Fig. 11(f,i,l)).
Turning back to case 111, we see the lee waves weakening as $z/D$ is increased
(Fig. 11(d,g,j)) but their pattern is still clear at $z/D=-30$.

###  4.3 Wake generated internal gravity waves

Figure 12: $\partial w^{\prime}/\partial z$ contours showing wake generated
internal gravity waves on the $\theta=90^{\circ}$ plane for 111. Black lines
are inclined at an angle of $39^{\circ}$ with the vertical.

The internal gravity waves generated by the turbulent wake are unsteady. These
waves are visualized in Fig. 12 for case 111 by contours of $\partial
w/\partial z$ on the $\theta=90^{\circ}$ vertical plane. The phase lines of
the radiated waves in the far field are seen to cluster around a a
characteristic inclination angle which suggests a narrow frequency band in the
far field according to the dispersion relationship, Eq. 10, for the intrinsic
(in a frame where the background has zero velocity) wave frequency. Fig. 12
also shows solid black lines plotted at an angle of $39^{\circ}$ from the
vertical, which also seems to be the angle for the waves ($39^{\circ}\pm
2^{\circ}$). This gives $\Omega/N=0.78\pm 0.03$.

Figure 13: $\partial w^{\prime}/\partial z$ contours showing wake generated
internal gravity waves on the $\theta=90^{\circ}$ and $180^{\circ}$plane for
$614+$.

For profile $614+$, the $\partial w^{\prime}/\partial z$ contours on the
entire vertical plane ($\theta=90^{\circ}$ and $270^{\circ}$) are plotted in
Fig. 13. Phase lines with a dominant inclination angle are seen in the
negative $z/D$ region where the waves with downward group velocity travel in
an entirely linear stratification ($\mbox{{Fr}}=3.74$). The inclination angles
in this region of $z/D<0$ cluster around $\Theta=44^{\circ}$ which gives
$\Omega/N\approx 0.72$ (Note that $N$ here corresponds to the local
$\mbox{{Fr}}=3.74$). The upward moving waves generated by the disk enter a
region of nonlinear stratification where there is significant small-scale
variability. As the waves exit the pycnocline layer from the top, they have
near-zero inclination with respect to the vertical, i.e. $\Theta\approx
0^{\circ}$, which implies that only the high-frequency near-$N$ ($N$ here
again corresponds to the local $\mbox{{Fr}}=6.13$) part of the wake generated
waves escapes into the upper region of weak stratification. Thus, a
significant portion of the wake generated waves is trapped within the
pycnocline leading to the small-scale variability in that region.

###  4.4 Wake characteristics

Figure 14: Evolution of centerline mean defect velocity for the five simulated
cases.

Wake characteristics, namely, the mean defect velocity and turbulent kinetic
energy, show qualitatively different behavior for the shifted pycnocline cases
$614-$ and $614+$ relative to the centered 614 profile. Fig. 14 compares the
streamwise evolution of centerline mean defect velocity ($U_{0}/U_{\infty}$)
among all five simulated cases. Recall that for $614+$, the disk is entirely
in the $\mbox{{Fr}}=3.74$ region with its upper edge in the pycnocline and,
for $614-$, half of the disk is in the $\mbox{{Fr}}=6.13$ region and half in
the pycnocline. Since, for both $614-$ and $614+$, the disk wake evolves in
relatively weaker stratification relative to the $Fr=1$ stratification seen by
the disk in 614, there is a faster rate of decay of $U_{0}$ ($\propto
x^{-0.35}$) in these cases (Fig. 14). $614-$ has a weaker effective
stratification relative to $614+$ and, therefore, has a somewhat lower
$U_{0}$.

Figure 15: Turbulent kinetic energy contours at $x/D=20$ for all five cases.

The TKE contours plotted in Fig. 15 also show that the buoyancy effect is
weaker for the shifted profiles. $614-$ and $614+$ show far higher TKE than
the other three cases that have $\mbox{{Fr}}=1$ at the disk center. Since the
bottom half of the disk of $614-$ is inside the pycnocline layer, the wake
core is not symmetric about the horizontal axis and appears to be somewhat
compressed from the bottom where the effective stratification is higher. The
TKE profile at $x/D=20$ is more symmetric for $614+$.

##  5 Summary and Conclusions

We present results from a body-inclusive LES of turbulent wakes in a
background stratification which is nonlinear instead of the linear
stratification that is typically considered. Four density profiles with a
hyperbolic tangent pycnocline are selected along with a standard constant
linear stratification profile. The nonlinear profiles are 614 (weak upper
stratification with $Fr\approx 6$ and weak lower stratification with
$Fr\approx 4$ bounding the pycnocline), I1I (no stratification surrounding
pycnocline), $614-$ (614 vertically shifted by $-2.5D$), and $614+$ (614
vertically shifted by $3D$). A disk of diameter $D$ is chosen as the canonical
bluff body, the relative velocity is $U_{\infty}$ and the Reynolds number is
$5000$. The Froude number $Fr(z)=U_{\infty}/N(z)D$, which is based on local
value of buoyancy frequency $N(z)$, varies and takes a minimum value of $1$
for the cases with nonuniform $N(z)$ and is set to 1 in the benchmark
constant-$N$ case. The inverse of $Fr(z)$ can be viewed as the local buoyancy
frequency, nondimensionalized with the flow scales, so that nondimensional
$N_{max}$ is the same among all cases. The pycnocline thickness normalized by
$D$ is also held constant.

Our main conclusion is that the non-constant $N$ profile studied here
substantially alters both the turbulent wake and the internal wave field. In
the first part of the study, results with the nonlinear profiles 614 and I1I
are compared with the linear profile 111. The mean defect velocity
($U_{0}/U_{\infty}$) in the near wake ($x/D<8$) is quite different for 111
which shows an initial oscillation with stronger amplitude than 614 and I1I.
This initial oscillation of the defect velocity is the oscillatory modulation
of the wake, which has been shown to be an imprint of the steady lee wave
field in linear stratification (Pal et al., 2017). The I1I case, which is
unstratified outside the pycnocline, does not support a steady lee wave in the
far field. Nevertheless, its near wake exhibits oscillatory modulation
(wavelength is slightly larger than that for 111) because of an evanescent
steady wave. In the NEQ regime, which follows the initial oscillation,
$U_{0}/U_{\infty}$ in all three cases was found to decay as $x^{-0.18}$.
Chongsiripinyo & Sarkar (2020) found the same power-law exponent of $-0.18$ in
their disk wake, which had an order of magnitude larger $\Rey=5\times 10^{4}$.
The decay rate transitions to $x^{-0.45}$ around $x/D=45$, suggesting the end
of the NEQ regime and beginning of the regime with Q2D power law. The Q2D
regime was not found by Chongsiripinyo & Sarkar (2020) at their lowest Froude
number, $\mbox{{Fr}}=2$.

Wakes formed within a pycnocline layer are found here to be more turbulent
than in linear stratification given the same value of $N(z)$ at the disk
center. The TKE for cases 614 and I1I exceeds that for 111 even though $N(z)$
is very similar in the wake core among the three cases. The difference is
specially large in the NEQ regime where I1I has about twice the TKE of 111.
The reasons are as follows. First, the most dominant production term,
$P_{xy}$, is higher for 614 and I1I. The turbulent momentum flux is sensitive
to the weak stratification, even though the weakening occurs away from and not
in the wake core and, therefore, the suppression of the turbulent flux by
stabilizing buoyancy is weaker in 614 and I1I. Second, the internal wave
energy flux out of the wake, $\langle p^{\prime}u_{r}^{\prime}\rangle$ is
smaller for 614 and I1I because there is a range of internal waves that is
generated but cannot propagate out vertically since their frequency exceeds
the smaller $N$ outside the pycnocline. This wave energy is trapped inside the
pycnocline layer. Both of these effects dominate at the beginning of the NEQ
regime, resulting in a more turbulent wake. It is worth noting that, a side
lobe of wave radiation forms in 614 and I1I (center and left panels of Fig.
7a), but the overall integrated wave flux in the NEQ regime is still smaller
than 111.

Based on the theoretical asymptotic analysis for steady lee waves by Voisin
(2007), the potential flow solution for an oval was used by Ortiz-Tarin et al.
(2019) to develop an expression for the lee wave field of a spheroid in linear
stratification. We find that the disk, along with its separation bubble, can
be approximated as a Rankine ovoid and have used the ovoid shape for the
potential flow solution instead of an oval. There is excellent agreement
between the 111 simulation and the analytical expression with regards to both
wave length and wave amplitude. The 614 simulation has a lee wave field whose
wavelength is well predicted by constant-$N$ linear theory upon using
$\mbox{{Fr}}=6.13$ above and $\mbox{{Fr}}=3.74$ below the pycnocline to
approximate the steady lee wave in those regions. However, the simulation
amplitude for 614 is higher than that given by this approximation since the
wave-generation region, which includes the pycnocline, has an effective $N$
that is larger then the constant far-field $N$ and, therefore, a stronger wave
field.

In the second part of the study, the influence of shifting the pycnocline
layer relative to the disk center is found to be strong. In 614-, half of the
disk lies in the upper region of weak ($\mbox{{Fr}}=6.13$) linear
stratification and, in 614+, the disk lies entirely in the lower linearly
stratified region ($\mbox{{Fr}}=3.74$) but still very close to the
nonuniform-$N$ region. Overall, the relative shift between pycnocline and disk
weakens the buoyancy effect on the mean wake, since the wake feels a weaker
effective stratification. Thus, different from the 614 centered case, the
defect velocity in the shifted cases do not show an oscillation in $x/D<8$ and
also decays more quickly than 614 in the far wake. The wakes are also more
turbulent in $614-$ and $614+$ as shown by the TKE. The asymmetric placement
(with respect to the profile) of the disk in $614-$ is also clearly manifested
in its TKE contour where the lower half, which is inside the pycnocline layer
is vertically compressed relative to the upper half.

With regards to the wave field, in addition to the usual steady far-field lee
waves and unsteady wake generated waves, the pycnocline supports a third
family of waves, the steady Kelvin wake waves. The Kelvin wake waves can be
analytically described by solving the Taylor-Goldstein equation as an
eigenvalue problem and they take a fundamentally different form than the lee
waves as visualized by Fig. 11. Furthermore, shifting the disk alters the
modal wave form. The dominant waveform in $614$ corresponds to the mode-2
eigenfunction while $614-$ and $614+$ exhibit a mix of two waveforms,
corresponding to mode-1 and mode-2 eigenfunction. In 614, the disk moves right
at the center of the pycnocline layer leading to the symmetry property of the
mode-2 wave. Any vertical offset from the center leads to the appearance of
mode-1 waves. These modal waveforms are well corroborated by the experiments
of Robey (1997) as well as Nicolaou et al. (1995).

Phase lines of unsteady internal gravity waves in the linearly stratified far
field are found to cluster around a characteristic inclination angle which
takes the value of $\theta\approx 39^{\circ}$ in 111. This result of phase-
angle clustering is consistent with several previous studies of wave radiation
from a turbulent flow (Sutherland & Linden, 1998; Dohan & Sutherland, 2003;
Taylor & Sarkar, 2007; Pham et al., 2009) but it is also worth noting that as
$\Rey$ increases, Abdilghanie & Diamessis (2013) find a broader range of wave
phase angles. In case $614+$, the wake generated waves that propagate downward
do so in a linear stratification and cluster around $\theta\approx
44^{\circ}$. The upward propagating waves have to go through the pycnocline
layer and only the near-$N$ waves are able to propagate in the weak linear
stratification above the pycnocline. Some of the wave energy that reflects off
the pycnocline boundaries is trapped and leads to small-scale variability
while a fraction escapes through the sides of the wake.

We have explored a limited portion of the parameter space applicable to
turbulent wakes in nonlinear background stratification. Depending on the
application and the environmental setting, there can be wide variability in
the nondimensional pycnocline thickness ($\delta/D$), the shape of the
nonuniform stratification, and the variability of $N(z)$. Body characteristics
such as its nominal Froude number, Reynolds number and its shape can vary.
Future examination of this wide parameter space in laboratory experiments as
well as simulations (body-inclusive hybrid spatial or temporal as well as the
cheaper body-exclusive temporal with a good choice for initial conditions)
will be useful.

[Funding]The authors gratefully acknowledge the support of Office of Naval
Research Grant N00014-20-1-2253.

[Declaration of interests] The authors report no conflict of interest.

## References

  * Abdilghanie & Diamessis (2013) Abdilghanie, Ammar M & Diamessis, Peter J 2013 The internal gravity wave field emitted by a stably stratified turbulent wake. Journal of Fluid Mechanics 720, 104–139.
  * Balaras (2004) Balaras, E. 2004 Modeling complex boundaries using an external force field on fixed Cartesian grids in large-eddy simulations. Comput. Fluids 33 (3), 375–404.
  * Barber (1993) Barber, BC 1993 On the dispersion relation for trapped internal waves. Journal of Fluid Mechanics 252, 31–49.
  * Bonneton et al. (1993) Bonneton, P, Chomaz, JM & Hopfinger, EJ 1993 Internal waves produced by the turbulent wake of a sphere moving horizontally in a stratified fluid. Journal of Fluid Mechanics 254, 23–40.
  * Brandt & Rottier (2015) Brandt, A & Rottier, JR 2015 The internal wavefield generated by a towed sphere at low froude number. Journal of Fluid Mechanics 769, 103–129.
  * Brucker & Sarkar (2010) Brucker, K. A. & Sarkar, S. 2010 A comparative study of self-propelled and towed wakes in a stratified fluid. J. Fluid Mech. 652, 373–404.
  * Chomaz et al. (1993) Chomaz, J. M., Bonneton, P. & Hopfinger, E. J. 1993 The structure of the near wake of a sphere moving horizontally in a stratified fluid. J. Fluid Mech. 254, 1–21.
  * Chongsiripinyo et al. (2017) Chongsiripinyo, Karu, Pal, Anikesh & Sarkar, Sutanu 2017 On the vortex dynamics of flow past a sphere at re= 3700 in a uniformly stratified fluid. Physics of Fluids 29 (2), 020704\.
  * Chongsiripinyo & Sarkar (2020) Chongsiripinyo, K. & Sarkar, S. 2020 Decay of turbulent wakes behind a disk in homogeneous and stratified fluids. J. Fluid Mech. 885, A31.
  * Diamessis et al. (2005) Diamessis, Peter J, Domaradzki, JA & Hesthaven, Jan S 2005 A spectral multidomain penalty method model for the simulation of high reynolds number localized incompressible stratified turbulence. Journal of Computational Physics 202 (1), 298–322.
  * Dohan & Sutherland (2003) Dohan, K & Sutherland, BR 2003 Internal waves generated from a turbulent mixed region. Physics of Fluids 15 (2), 488–498.
  * Dommermuth et al. (2002) Dommermuth, Douglas G, Rottman, James W, Innis, George E & Novikov, Evgeny A 2002 Numerical simulation of the wake of a towed sphere in a weakly stratified fluid. Journal of Fluid Mechanics 473, 83–101.
  * Ermanyuk & Gavrilov (2002) Ermanyuk, Eugeny V & Gavrilov, Nikolai V 2002 Force on a body in a continuously stratified fluid. part 1. circular cylinder. Journal of Fluid Mechanics 451, 421.
  * Ermanyuk & Gavrilov (2003) Ermanyuk, Eugeny V & Gavrilov, Nikolai V 2003 Force on a body in a continuously stratified fluid. part 2. sphere. Journal of Fluid Mechanics 494, 33.
  * Germano et al. (1991) Germano, M., Piomelli, U., Moin, P. & Cabot, W. H. 1991 A dynamic subgrid‐scale eddy viscosity model. Phys. Fluids 3 (7), 1760–1765.
  * Gilreath & Brandt (1985) Gilreath, HE & Brandt, A 1985 Experiments on the generation of internal waves in a stratified fluid. AIAA journal 23 (5), 693–700.
  * Gourlay et al. (2001) Gourlay, M. J., Arendt, S. C., Fritts, D. C. & Werne, J. 2001 Numerical modeling of initially turbulent wakes with net momentum. Phys. Fluids 13 (12), 3783–3802.
  * Grisouard et al. (2011) Grisouard, Nicolas, Staquet, Chantal & Gerkema, Theo 2011 Generation of internal solitary waves in a pycnocline by an internal wave beam: a numerical study. Journal of Fluid Mechanics 676, 491–513.
  * Hopfinger et al. (1991) Hopfinger, Emil J, Flor, J-B, Chomaz, J-M & Bonneton, P 1991 Internal waves generated by a moving sphere and its wake in a stratified fluid. Experiments in Fluids 11 (4), 255–261.
  * Keller & Munk (1970) Keller, Joseph B & Munk, Walter H 1970 Internal wave wakes of a body moving in a stratified fluid. The Physics of Fluids 13 (6), 1425–1431.
  * Lin & Pao (1979) Lin, Jung-Tai & Pao, Yih-Ho 1979 Wakes in stratified fluids. Annual Review of Fluid Mechanics 11 (1), 317–338.
  * Meunier et al. (2018) Meunier, P, Le Dizès, Stéphane, Redekopp, L & Spedding, GR 2018 Internal waves generated by a stratified wake: experiment and theory. Journal of Fluid Mechanics 846, 752–788.
  * Nicolaou et al. (1995) Nicolaou, D, Garman, JFR & Stevenson, TN 1995 Internal waves from a body accelerating in a thermocline. Applied scientific research 55 (2), 171–186.
  * Nidhan et al. (2022) Nidhan, Sheel, Schmidt, Oliver T & Sarkar, Sutanu 2022 Analysis of coherence in turbulent stratified wakes using spectral proper orthogonal decomposition. Journal of Fluid Mechanics 934.
  * Orlanski (1976) Orlanski, I. 1976 A simple boundary condition for unbounded hyperbolic flows. J. Comput. Phys. 21 (3), 251–269.
  * Orr et al. (2015) Orr, Trevor S, Domaradzki, J Andrzej, Spedding, Geoffrey R & Constantinescu, George S 2015 Numerical simulations of the near wake of a sphere moving in a steady, horizontal motion through a linearly stratified fluid at re= 1000. Physics of Fluids 27 (3), 035113.
  * Ortiz-Tarin et al. (2019) Ortiz-Tarin, Jose L, Chongsiripinyo, KC & Sarkar, S 2019 Stratified flow past a prolate spheroid. Physical Review Fluids 4 (9), 094803.
  * Pal et al. (2016) Pal, Anikesh, Sarkar, Sutanu, Posa, Antonio & Balaras, Elias 2016 Regeneration of turbulent fluctuations in low-froude-number flow over a sphere at a reynolds number of 3700. Journal of Fluid Mechanics 804.
  * Pal et al. (2017) Pal, Anikesh, Sarkar, Sutanu, Posa, Antonio & Balaras, Elias 2017 Direct numerical simulation of stratified flow past a sphere at a subcritical reynolds number of 3700 and moderate froude number. Journal of Fluid Mechanics 826, 5.
  * Pham & Sarkar (2010) Pham, Hieu T & Sarkar, Sutanu 2010 Internal waves and turbulence in a stable stratified jet. Journal of fluid mechanics 648, 297–324.
  * Pham et al. (2009) Pham, Hieu T, Sarkar, Sutanu & Brucker, Kyle A 2009 Dynamics of a stratified shear layer above a region of uniform stratification. Journal of fluid mechanics 630, 191–223.
  * Redford et al. (2015) Redford, JA, Lund, TS & Coleman, GN 2015 A numerical study of a weakly stratified turbulent wake. Journal of Fluid Mechanics 776, 568.
  * Redford et al. (2012) Redford, J. A., Castro, I. P. & Coleman, G. N. 2012 On the universality of turbulent axisymmetric wakes. J. Fluid Mech. 710, 419–452.
  * Robey (1997) Robey, Harry F 1997 The generation of internal waves by a towed sphere and its wake in a thermocline. Physics of Fluids 9 (11), 3353–3367.
  * Rossi & Toivanen (1999) Rossi, T. & Toivanen, J. 1999 A Parallel Fast Direct Solver for Block Tridiagonal Systems with Separable Matrices of Arbitrary Dimension. SIAM J. Sci. Comput. 20 (5), 1778–1793.
  * Rowe et al. (2020) Rowe, KL, Diamessis, PJ & Zhou, Q 2020 Internal gravity wave radiation from a stratified turbulent wake. Journal of Fluid Mechanics 888.
  * Spedding (1997) Spedding, GR 1997 The evolution of initially turbulent bluff-body wakes at high internal froude number. Journal of Fluid Mechanics 337, 283–301.
  * de Stadler & Sarkar (2012) de Stadler, Matthew B & Sarkar, Sutanu 2012 Simulation of a propelled wake with moderate excess momentum in a stratified fluid. Journal of fluid mechanics 692, 28–52.
  * de Stadler et al. (2010) de Stadler, Matthew B, Sarkar, Sutanu & Brucker, Kyle A 2010 Effect of the prandtl number on a stratified turbulent wake. Physics of Fluids 22 (9), 095102.
  * Stevenson et al. (1986) Stevenson, TN, Kanellopulos, D & Constantinides, M 1986 The phase configuration of trapped internal waves from a body moving in a thermocline. Applied scientific research 43 (2), 91–105.
  * Sturova (1974) Sturova, IV 1974 Wave motions produced in a stratified liquid from flow past a submerged body. Journal of Applied Mechanics and Technical Physics 15 (6), 796–805.
  * Sutherland (2016) Sutherland, Bruce R 2016 Internal wave transmission through a thermohaline staircase. Physical Review Fluids 1 (1), 013701.
  * Sutherland & Linden (1998) Sutherland, Bruce R & Linden, Paul F 1998 Internal wave excitation from stratified flow over a thin barrier. Journal of Fluid Mechanics 377, 223–252.
  * Sutherland & Yewchuk (2004) Sutherland, Bruce R & Yewchuk, Kerianne 2004 Internal wave tunnelling. Journal of Fluid Mechanics 511, 125–134.
  * Taylor & Sarkar (2007) Taylor, John R & Sarkar, Sutanu 2007 Internal gravity waves generated by a turbulent bottom ekman layer. Journal of Fluid Mechanics 590, 331–354.
  * Voisin (1991) Voisin, Bruno 1991 Internal wave generation in uniformly stratified fluids. part 1. green’s function and point sources. Journal of Fluid Mechanics 231, 439–480.
  * Voisin (1994) Voisin, Bruno 1994 Internal wave generation in uniformly stratified fluids. part 2. moving point sources. Journal of Fluid Mechanics 261, 333–374.
  * Voisin (2003) Voisin, Bruno 2003 Limit states of internal wave beams. Journal of Fluid Mechanics 496, 243–293.
  * Voisin (2007) Voisin, Bruno 2007 Lee waves from a sphere in a stratified flow. Journal of Fluid Mechanics 574, 273–315.
  * Voropayeva & Chernykh (1997) Voropayeva, OF & Chernykh, GG 1997 Numerical model of the dynamics of a momentumless turbulent wake in a pycnocline. Journal of applied mechanics and technical physics 38 (3), 391–406.
  * Yang & Balaras (2006) Yang, J. & Balaras, E. 2006 An embedded-boundary formulation for large-eddy simulation of turbulent flows interacting with moving boundaries. J. Comput. Phys. 215 (1), 12–40.
  * Zhou & Diamessis (2019) Zhou, Qi & Diamessis, Peter J 2019 Large-scale characteristics of stratified wake turbulence at varying reynolds number. Physical Review Fluids 4 (8), 084802.