Decay of turbulent wakes behind a disk in homogeneous and stratified fluids

Body-inclusive large-eddy simulations of disk wakes are performed for a homogeneous fluid and for different levels of stratification. The Reynolds number is 5 × 104 and the Froude number ($Fr$) takes the values of $\infty$, 50, 10 and 2. In the axisymmetric wake of a disk with diameter $L_{b}$ in a homogeneous fluid, it is found that the mean streamwise velocity deficit ($U_{0}$) decays in two stages: $U_{0}\propto x^{-0.9}$ during $10<x/L_{b}<65$ and, subsequently, $U_{0}\propto x^{-2/3}$. Consequently, none of the simulated stratified wakes is able to exhibit the classical 2/3 decay exponent of $U_{0}$ in the interval before buoyancy effects set in. Stratification affects the wake within approximately one buoyancy time scale, after which, we find three regimes: weakly stratified turbulence (WST), intermediately stratified turbulence (IST) and strongly stratified turbulence (SST). WST begins when the turbulent Froude number ($Fr_{h}$) decreases to $O(1)$, spans $1\lesssim Nt_{b}\lesssim 5$ and, while the mean flow is strongly affected by buoyancy in WST, turbulence is not. During IST, which commences at $Nt_{b}\approx 5$ when $Fr_{h}=O(0.1)$, the mean flow has arrived into the non-equilibrium (NEQ) regime with $U_{0}\propto x^{-0.18}$, but the turbulence state is still in transition, as indicated by progressively increasing turbulence anisotropy. When $Fr_{h}\sim O(0.01)$ at $Nt_{b}\approx 20$, the wake transitions into SST, where the turbulent vertical Froude number ($Fr_{v}$) asymptotes to a $O(1)$ constant. There is strong anisotropy ($u_{z}^{\prime }\ll u_{h}^{\prime }$), and both $u_{h}^{\prime }$ and $U_{0}$ satisfy $x^{-0.18}$ decay, signifying the arrival of the NEQ regime for both turbulence and mean flow. Turbulence is patchy and temporal spectra are broadband in the SST wake. The wake height decreases as $L_{V}\sim O(U_{0}/N)$ in IST/SST. Energy budgets reveal that stratification prolongs wake life during WST/early-IST by both an energy transfer from mean potential energy to mean kinetic energy and reduction of turbulent production. In the late-IST/early-SST stages, production is enhanced and, additionally, there is injection from turbulent potential energy slowing down turbulent kinetic energy (TKE) decay. Only in the SST stage, when NEQ is realized for both the mean and turbulence, does the turbulent buoyancy flux become negative again, acting as a sink of TKE.

at low Fr (Pal et al. 2016;Ortiz-Tarin, Chongsiripinyo & Sarkar 2019) and the associated change in vortex dynamics . These features of the flow could not have been captured with a temporal model. The statistics to be presented in this work demonstrate the impact of these features on wake evolution.
Stratified wakes are inhomogeneous turbulent flows with mean shears that are subject to stratification. There has been much work -a non-exhaustive list is Lilly (1983), Billant & Chomaz (2001), Riley & de Bruyn Kops (2003), Brethouwer et al. (2007), Kimura & Herring (2012), Maffioli & Davidson (2016) and de Bruyn Kops & Riley (2019) -in the related problem of fluctuations that evolve without mean shear in a stratified fluid. The interest is in the behaviour of 'stratified turbulence' defined as fluctuations that have low values of Froude number (Fr h defined by (3.4) is <1) but the Reynolds number is not low (Re b defined by (3.3) is >1).
A limitation of the previous body-inclusive simulations is that Re = O(10 3 ) was not sufficiently large to have a developed region of stratified turbulence with fluctuations at low Froude number and high Reynolds number. Furthermore, wakes with high Fr (>O(10)) were not simulated. We are therefore motivated to simulate wakes at a higher Re = 5 × 10 4 and over a wide range of stratifications, {Fr = 2, 10, 50, ∞}. We will explore links between the findings in our simulations and those in the general topic of stratified turbulence. Furthermore, we consider a disk rather than a sphere to broaden the stratified-wake literature. The unstratified Fr = ∞ case allows examination of power laws at a higher Re than in prior simulations and assesses the possibility of non-canonical power laws for U 0 and L.
The primary objective of the present study is to improve our understanding of the behaviour of relatively high-Re stratified wakes. In particular, we consider the following questions. What are the power laws that are satisfied by characteristic velocities/length scales as the wake progresses? Are there differences between the progression of mean and r.m.s. turbulence as they react to stratification? How does the development of wake turbulence relate to the broader area of stratified turbulence decay? Lastly, what are the reasons underlying the slower decay of the wake in the NEQ regime?
The numerical set-up and the parameters of the simulated cases are given in § 2. Definitions used in the analysis and interpretation of results are given in § 3. Presentation of the results begins with § 4, which introduces the effect of stratification by visual contrasts among the different cases. Section 5 reports on the evolution of centreline values of the mean and the three r.m.s. components for each of the cases, separately. Section 6 presents a consolidated picture of the three stratified cases in phase space. We discuss links with previous work on regimes (characterization and transitions) of stratified wakes and stratified turbulence in both § § 5 and 6. Section 7 concerns the behaviour of two distinct wake sizes; one is derived from profiles of U 0 and the other is based on TKE. Section 8 discusses the progression of area-integrated kinetic and potential energy in stratified wakes and contrasts with the unstratified wake. Section 9 delves deeper into the evolution of mean kinetic energy (MKE) and TKE through their budgets. Section 10 discusses possible scaling laws for the dissipation (ε) of TKE. We end with a summary and conclusions in § 11.

Numerical simulations
We address the questions introduced in § 1 with large-eddy simulations (LES) of flow past a circular solid disk placed perpendicular to the uniform free stream. The background is taken to have a uniform stratification. These body-inclusive simulations of flow into the far wake (up to x/D = 125) are at a relatively high Re of 50 000 and computationally intensive because of the necessity of resolving the separated flow at the body and the turbulent recirculation region. The following non-dimensional filtered Navier-Stokes equations under the Boussinesq assumption for density effects along with the filtered advection-diffusion equation for density are numerically solved: continuity: ∂ i u i = 0, (2.1) momentum: ∂ t u i + u j ∂ j u i = −∂ i p + Re −1 ∂ j [(1 + ν s /ν)∂ j u i ] − Fr −2 ρ d δ i3 , (2.2) density: ∂ t ρ + u j ∂ j ρ = (RePr) −1 ∂ j [(1 + κ s /κ)∂ j ρ]. (2. 3) The symbol ∂ i denotes a spatial derivative with respect to x i , where i = 1, 2, 3 respectively refer to streamwise (x), lateral (y) and vertical (z) directions with u x , u y and u z the corresponding velocities. The governing equations are non-dimensionalized with the following parameters: the background free-stream velocity (U b ), the disk diameter referred as the body length scale (L b ), advection time (L b /U b ), dynamic pressure (ρ 0 U 2 b ) and characteristic change in background density deviation across the body (−L b ∂ 3 ρ b ). Under the Boussinesq approximation, the density deviation from its equilibrium alters the momentum equation only through the last term on the right-hand side of (2.2). Here δ i3 is the Dirac delta function. The density is decomposed into a constant reference density (ρ 0 ), the linearly varying deviation of the background ρ b (x 3 ) and the flow-induced deviation, ρ d (x i , t), as follows: (2.4) is not necessarily zero but takes a value of ∂ n ρ d (x i , t) = −∂ n ρ b (x 3 ) at a solid surface; n is a surface normal direction. Background density and static linear variation are absorbed into the modified pressure. The body Reynolds number is Re = U b L b /ν, where U b is the free-stream velocity, L b is the characteristic length taken to be the disk diameter, and ν is the fluid kinematic viscosity. The body Froude number Fr = U b /NL b is the ratio of buoyancy time scale to the characteristic advection time scale, L b /U b ; here, N is the constant buoyancy frequency defined by N 2 = −(g/ρ o )∂ 3 ρ b . The background is stably stratified if ∂ 3 ρ b is negative, as is the case here. The disk thickness of 0.01L b is small. The Prandtl number, Pr = ν/κ, that is, the ratio of velocity and density (temperature) diffusivities, is assumed to be unity. Additional variables from small unresolved scales are subgrid kinematic viscosity, ν s , and density diffusivity, κ s . The Prandtl number based on these subgrid variables is also assumed to be unity. The study of the Prandtlnumber effects in a stratified wake can be found in de Stadler, Sarkar & Brucker (2010).
Since the turbulent wake under weak-to-intermediate stratification is quasiaxisymmetric in the 3-D and in the early-NEQ regimes, a cylindrical coordinate system is employed. The choice of cylindrical coordinates allows efficient distribution of grid points especially from the core to the wake periphery. Spatial numerical derivatives are obtained with second-order-accurate central differences, while temporal marching is done with a fractional-step method that combines a third-order Runge-Kutta explicit scheme with the second-order Crank-Nicolson implicit scheme. To alleviate stiffness of the discretized system, especially near the coordinate streamwise axis, implicit marching is performed for viscous and advection terms (both velocities and density) that involve spatial discretization in the azimuthal direction. In the fractional-step method, the Poisson equation is formed by taking the divergence Case no. Re Note that all lengths are normalized by the diameter, L b , of the disk.
of the momentum equation for predicted velocity. The periodic boundary condition in the azimuthal direction transforms the discretized Poisson equation into inversion of a pentadiagonal matrix. The pentadiagonal matrix system is then inverted using a direct solver (Rossi & Toivanen 1999). The disk is represented by the immersed boundary method of Balaras (2004) and Yang & Balaras (2006). Kinematic subgrid viscosity, ν s , is obtained using the eddy viscosity model of Germano et al. (1991), a variant of the dynamic Smagorinsky model (Smagorinsky 1963). The coefficient C, as in ν t = C ∆ 2 | S|, where ∆ 3 = V is a measure of the local cell volume and | S| is the instantaneous strain-rate magnitude of filtered velocity, is dynamically computed using a method of Lilly (1992) that takes the cumulative average of C along flow trajectories with an exponential weighting function chosen to give more weight to recent times in flow history. Boundary conditions at the inlet and outlet of the domain are Dirichlet inflow and Orlanski-type convective outflow, ∂ t φ + C∂ 1 φ = 0. The magnitude of convective velocity, C, is set at U b . Note that a localized mean value of C as an alternative to U b is tested with negligible consequence. The Neumann condition is used at the outer radial boundary. Unlike  and Pal et al. (2017), who use a 'sponge' region for all of their cases to absorb internal gravity waves, the domain's radial extent in this present study is enlarged to 80L b , allowing free propagation of internal gravity waves while minimizing reflection from the boundary without artificially adjusting the far field to a background state. In the present study, a sponge region of a Rayleigh-damping functional form, which gradually relaxes the velocities and density to their respective values at the boundaries, is applied only in the F02 case.
Parameters of the simulations are given in table 1. The distribution of the computational grid is optimized on the basis of the unstratified case (UNS) where the turbulent dissipation rate is known to decay as ε ∝ x −m , m 5/2, in the self-similar solution (George 1989), thus leading to an estimate of the Kolmogorov length. The streamwise-dependent radial distribution of ε is estimated from the direct numerical simulation (DNS) of Dairay et al. (2015) and the LES of . The resulting number of grid points for the weakly stratified cases is 624 million cells. In the unstratified wake, x/η is at worst ≈17 near the centreline at x/L b ≈ 2.5; here η is the Kolmogorov scale calculated directly from the turbulent dissipation rate by η = (Re 3 ε) −1/4 . Beyond approximately x/L b = 10, the centreline x/η is smaller than 10 and gradually decreases to approximately 6 at x/L b = 125. The distributions of x/η and r/η are given in figure 1. To ensure adequacy of the high spatial resolution utilized here, a grid sensitivity study is conducted. In the grid sensitivity study, the grid spacings in the azimuthal and Decay of stratified turbulent wakes behind a disk 885 A31-7 10 1 10 0 10 1 10 0 10 2 10 1 10 0 10 -1 FIGURE 1. Grid quality. (a) Streamwise variation of non-dimensional streamwise grid spacing. (b) Radial profiles of non-dimensional radial grid spacing. The Kolmogorov length scale, η = (ν 3 /ε k ) 1/4 , is computed using the rate of TKE dissipation (ε k ) of the unstratified wake.
streamwise directions are both increased by a factor of 2 (N θ = 128 and N x = 2304). None of the qualitative or quantitative results on wake statistics are altered. The approximate time step for all the cases lies in the range (2-9) × 10 −4 L b /U b . Statistics are collected over an averaging time (T) during a statistically steady state determined by convergence of moving-window averages of centreline streamwise and radial velocities at x/L b = 100. The steady state is reached after an initial transient that takes approximately two flow-through times. For the unstratified case, we use an averaging time of T = 500L b /U b after the transient has subsided. The effect of averaging time on the statistics is found to be small to negligible as the averaging time is changed from one to two flow-through times. The CPU usage, between 750 000 and 10 6 core hours distributed over 512 cores for each case, is large because of the long integration time required to reach statistical steady state and then obtain converged statistics. The averaged drag coefficient in the unstratified case is found to be C d = 1.145, comparable to the value of C d = 1.12 in Fail, Eyre & Lawford (1959).

Analysis and interpretation
3.1. Statistics Once statistical steady state is reached, the dependent variables are Reynoldsdecomposed, φ = φ + φ , where φ is an instantaneous realization. Representing an ensemble average, the angle bracket · denotes an appropriate averaging operator that involves averaging in all applicable homogeneous directions. Simulated in a spatially evolving domain, the turbulent wake in a homogeneous fluid is statistically homogeneous in time and in the azimuthal direction. The stratified wake is statistically homogeneous in time and is statistically symmetric across vertical and horizontal centre planes, φ (x, y, z) = φ (x, ±y, ±z). The other type of statistic is obtained from cross-wake area integration indicated by braces {·}, where the wake width (L k ) is described in § 3.2 and dA = θ r dr with θ = 2π/N θ . In order to include only wake turbulence and exclude internal wave contributions, a small integral number (2L K ) of wake widths is chosen for the radial limits of the integral. Lee waves, especially at intermediate stratification, generated by a wake generator can propagate considerably into the far field and we prefer to exclude the far field in the area integral.
3.2. Mean and turbulent length scales In a homogeneous fluid, an axisymmetric half-length L measures an azimuthally averaged distance from the centreline to a position where the streamwise mean velocity deficit is reduced to half of its centreline value, U 0 | r=L = U 0 | r=0 /2. The other half-length, L k , is based on profiles of the TKE, K = u i u i /2, as opposed to U 0 , taking K| r=L k = K| r=0 /2. In a stratified fluid, horizontal and vertical length scales are defined separately but again based on U 0 and K profiles, and denoted as L H , L V , L Hk , L Vk . (3.2) Here, the horizontal width (L H ) and vertical height (L V ) represent half-lengths in the lateral (y) and vertical (z) directions based on U 0 , while L Hk and L Vk are derived from profiles of the TKE. Since U 0 and K are obtained from the · operator, a half-length is calculated based on an appropriately averaged ensemble. Note that, at steady state of all cases, although instantaneous wakes can meander around a mean position, the centreline peak values of U 0 and K remain close to the axis of the cylindrical grid. The turbulent integral length scale is not easy to calculate in a spatially evolving stratified wake where spatial homogeneity is absent. We use the wake width (L Hk ) as a surrogate for the horizontal integral length scale. The turbulent vertical length scale, l v , is calculated along the centreline and its method of calculation is adopted from Riley & de Bruyn Kops (2003) The quantity, l v , has dynamic significance to shear instability. As shown by the scaling analysis of Riley & de Bruyn Kops (2003), Ri v is the Richardson number of the fluctuating shear in the flow.

Reynolds numbers and Froude numbers
The buoyancy Reynolds number (Re b ), the horizontal-motion Reynolds number (Re h ) and the microscale Reynolds number (Re λ ) are defined as where u h = (u 2 x + u 2 y ) 1/2 is the r.m.s. horizontal velocity fluctuation and λ is the Taylor microscale defined by λ 2 = 15νu 2 x /ε k , with ε k = ν ∂ j u i ∂ j u i + ν ∂ j u i ∂ i u j − τ s ij ∂ j u i the rate of TKE dissipation. The mean vertical Froude number (Fr V ), the mean horizontal Froude number (Fr H ), the turbulent vertical Froude number (Fr v ) and the turbulent horizontal Froude number (Fr h ) are defined as follows:

Visualization
Before considering the detailed analysis of wake statistics, we begin with a visualization that contrasts flow structures of the stratified wakes with those under unstratified conditions. Snapshots of the streamwise velocity (u x ) for the Fr = ∞, 10 and 2 wakes are shown in figures 2, 3 and 4, respectively. The body moving through the background leaves its imprint, which grows and lasts for a long distance, with its visual presentation being strongly dependent on the strength of stratification.   When the stratification is weak (Fr = 10), the near wake remains relatively unchanged relative to Fr = ∞ but the anisotropy between vertical and horizontal cuts is readily apparent in the intermediate and far wakes. When the stratification is strong (Fr = 2), the wake quickly responds to the background stratification; even the near wake is vertically contracted and the intermediate/far wake is more coherent than at Fr = ∞.
As buoyancy forces become relatively stronger, characterized by the decreased local Froude number, the far wake of the Fr = 10 case behaves similarly to the near wake of the Fr = 2 case by becoming horizontally wide but vertically thin. In fact, the Fr = 2 and Fr = 10 wakes are qualitatively equivalent if we compare the two at a normalized distance given by (x/L b )Fr −1 = Nt b . This aspect will be quantitatively elaborated in § § 5 and 6. As vertical motions become progressively suppressed, the wakes become 20 40 60 80 100 120 FIGURE 5. Unstratified wake. Evolution of centreline values of mean streamwise velocity deficit (U 0 ), r.m.s. velocity fluctuation (streamwise u x , spanwise u y , vertical u z ) and turbulent velocity (K 1/2 ). The inset shows the evolution of turbulence velocity quantities normalized with the local U 0 (x).
increasingly two-dimensional, with the appearance of horizontal waviness that drives the formation of large-scale horizontal but vertically thin 'pancake' eddies that will emerge later in the Q2D regime.

Evolution of mean deficit velocity and turbulence intensities
Centreline values of mean streamwise velocity deficit (U 0 ) and r.m.s. velocity fluctuations in the streamwise (u x ), spanwise (u y ) and vertical (u z ) directions are shown in figure 5 for the unstratified wake. The initial increase of U 0 in the recirculation zone to exceed the free-stream velocity (U b ) is accompanied by an increase in all r.m.s. components signifying the establishment of turbulence. The r.m.s. fluctuations peak at x/L b ≈ 2.5, which lies in the recirculation region. The U 0 /U b value decays over 2 x/L b 10 with a rate that gradually decreases with increasing x. The subsequent evolution of U 0 in figure 5 reveals a break in slopes at x/L b ≈ 65 that separates two stages with different power-law exponents. The first stage of 10 < x/L b < 65 exhibits approximately U 0 ∝ x −0.9 while the second stage of 65 < x/L b < 125 satisfies U 0 ∝ x −0.6 , close to the classical U 0 ∝ x −2/3 behaviour. A similar transition between two different power laws was found by Dairay et al. (2015) in the axisymmetric wake of a fractal plate. They identified U 0 ∝ x −1 (an exponent of −0.94 in their Re = 5000 DNS and −1.03 in their Re = 50 000 experiment) for 10 < x/L b < 50 that transitions to a different power law reported as U 0 ∝ x −0.81 in the Re = 5000 DNS. In our previous simulations of sphere wakes, namely DNS at Re = 3700 by Pal et al. (2017) and LES at Re = 10 000 by Chongsiripinyo et al. (2019), U 0 showed behaviour close to x −1 scaling. Although it is tempting to attribute the x −1 power law of U 0 to low-Re viscous decay (George 1989), that is not the case. As we will show later during the analysis of MKE, the turbulent production that acts as a sink of MKE is much larger in magnitude than the viscous dissipation term in the MKE equation. Furthermore, the microscale Reynolds number, which varies between Re λ = 200 at x/D = 10 and Re λ = 120 at x/D = 100, is not small.  In contrast to U 0 , the behaviour of the turbulent velocity scale (K 1/2 ) does not break into two separate power laws. The decay of K 1/2 (green circles in figure 5) is found to be K 1/2 ∝ x −0.71 , close to x −2/3 from x/L b = 10 to the end of the computational domain. Thus, the mean velocity scale in the intermediate wake (10 < x/L b < 65) decays with a different power-law exponent (−0.9 as shown in figure 5) than the exponent satisfied by the turbulence velocity scale. Similarity theory for the turbulent wake (Tennekes & Lumley 1972) assumes the same velocity scale for both mean and turbulence. In the present simulation, it is only beyond x/L b = 65 that the power-law exponents for U 0 and K 1/2 become similar and close to the classical value of −2/3. Near-wake turbulence is found to be anisotropic, with u x substantially smaller than u y and u z near the body; however, the r.m.s. velocity components (inset of figure 5) become more isotropic for x/D > 40.
The behaviour of U 0 in case F50 (figure 6a) shows little deviation from the UNS case until x/L b 50 or Nt b 1, consistent with theoretical arguments such as given by Riley & Lelong (2000) that it takes a time interval of Nt ∼ 1 for buoyancy forces to become comparable to inertial forces in a flow that originates with weak buoyancy effects. Figure 6 emphasizes the importance of Nt b = 1 by showing that U 0 /U 0∞ (figure 6b) abruptly increases and u /U 0 (figure 6c) reaches the peak at Nt b ≈ 1. Here, U 0∞ denotes the centreline deficit in the UNS case. Figure 6(b) also shows that stratification increases the rate of decay of U 0 as early as x/L b = 10 or Nt b = 0.2 (observe that U 0 /U 0∞ starts decreasing around x/L b = 10). Thus, there is a mild effect of buoyancy in the weakly stratified wake before Nt b = 1. Buoyancy increases the decay rate of K 1/2 (green circles) relative to the unstratified case (dashed green line). Figure 6(b) shows that all three turbulence components in F50 are reduced with respect to UNS, but the difference relative to UNS is smaller than that in U 0 . The delay in the turbulence response to stratification in comparison to the mean velocity 885 A31-12 K. Chongsiripinyo and S. Sarkar  deficit will be more apparent in the F10 case and the reason for this delay will be made clear in the discussion of the F02 case.
In the F10 case (figure 7), buoyancy forces become comparable with inertial forces at a location closer to the body relative to F50 but at a similar value of Nt b . The U 0 value deviates from the unstratified case at x/L b ≈ 10 (equivalently, Nt b ≈ 1) and, as emphasized by figure 7(b,c), U 0 /U 0∞ (b) increases sharply and u /U 0 (c) reaches the peak at x/L b ≈ 10. There is a strong effect of buoyancy on r.m.s. velocity fluctuations but it is delayed until x/L b ≈ 50 when, as shown in figure 7(b), the r.m.s. values of the horizontal components increase while the vertical r.m.s. value decreases relative to UNS. It is clear that the stratified wake exhibits a regime between Nt b ≈ 1 and Nt b ≈ 5 wherein the effect of buoyancy on the mean velocity is much stronger than that on turbulence. Beyond Nt b ≈ 5, the Reynolds stresses deviate strongly from isotropy towards a 'pancake' with u x ≈ u y u z .
The behaviour of Fr = O(1) wakes is quite different in the near and intermediate wakes from the Fr O(10) cases that have been discussed so far. Consider U 0 at Fr = 2 (case F02) in figure 8(a) together with the evolution of Froude numbers in figure 9. The mean recirculation region (which has negative velocity and deficit velocity U 0 > 1) is shorter because the incoming fluid that is vertically displaced by the body has sufficient buoyancy so that it plunges back in the separated flow towards the centreline. In contrast to F10 and F50, there is a region after the recirculation region where U 0 increases instead of decreasing. There is a local minimum of U 0 at Nt b π a half lee-wave period behind the disk; subsequently, U 0 increases and reaches a local maximum at Nt b 2π. This behaviour is the manifestation of the lee-wave-induced 'oscillatory modulation' of U 0 reported in the DNS of the Re = 3700 sphere wake by Pal et al. (2017). In the NEQ regime, U 0 is found to decay with a rate that is close to x −0.18 . The rate of decay is smaller than the U 0 ∝ x −0.25 behaviour in the sphere-wake experiments of Spedding (1997).
We now turn to the Froude numbers. Consider Fr V = U 0 /NL V based on mean wake deficit velocity. The value of Fr V in the case F02 (red curve in figure 9c) reaches 1 at x/D ≈ 5, a location at which U 0 commences a rapid deviation from the unstratified result. The Fr V value remains close to unity further downstream. Interestingly, in Decay of stratified turbulent wakes behind a disk 885 A31-13 10 0 10 -1 10 -2 10 1 10 1 10 0 10 2  The role of Froude number, Fr v (based on fluctuation rather than mean-flow statistics), in the evolution of r.m.s. turbulence is analogous to that of Fr V in the evolution of the mean flow. When Fr v decreases to O(1) (x/L b ≈ 20), the evolution of r.m.s. turbulence deviates strongly from unstratified behaviour since buoyancy becomes comparable to inertial forces in the inertial-range turbulent motions.
There is a distinct region in all the stratified cases where Fr V = O(1) so that the mean is strongly affected by buoyancy but Fr v > O(1) so that r.m.s. turbulence is not, for example, between x/L b = 5 and x/L b = 20 for case F02 (figure 8a). Since the Reynolds number (both Re λ and Re b ) is sufficiently high in this distinct region between x/L b = 5 and 20 that has Fr v > O(1), the decay of turbulence is close to the classical Kolmogorov decay law for unstratified turbulence, u 2 ∝ t −10/7 , as can be seen by comparing the evolution of K 1/2 with the x −5/7 line in figure 8(a). Despite the weak effect of buoyancy on turbulence, U 0 in this same distinct region exhibits strong buoyancy effects since Fr V has decreased to O(1). In the vicinity of Fr v 1 and beyond, the decay rates of u h and K 1/2 are reduced with respect to their unstratified counterparts. Once their decay has 'settled' in the NEQ regime, it is found that K 1/2 ∝ x −0.18 . In other words, turbulence and mean evolve with the same power law. This can be regarded as the official arrival of the NEQ regime. The value of u z continues to monotonically decrease with an approximate power-law decay of u z ∝ x −1 , agreeing with the t −1 scaling observed in the experiment of Lin & Pao (1979) and close to the t −1.08 scaling observed in the numerical simulation of . Notice that u y > u x in the late wake, a result that is consistent with coherent vortex patches that are located laterally off-centreline.

Turbulence in phase space
The Fr h -Re h Fr 2 h phase presents a consolidated view of the progression of wake turbulence under the influence of stratification. The phase-space map (figure 10) of wake evolution shows the progression of each of the simulated cases through different stages, which, as elaborated below, exhibit qualitatively different effects of buoyancy. Decreasing Fr h implies an increasing effect of buoyancy. Large scales of motion are preferentially affected by buoyancy and, as Fr h decreases, the largest scale affected by buoyancy also decreases. The quantity Re h Fr 2 h has been introduced by Billant & Chomaz (2001) as a parameter that must be >1 to prevent viscous effects from directly affecting turbulent motions, and Riley & de Bruyn Kops (2003) relate Re h Fr 2 h to an inverse Richardson number of fluctuating motions so that Re h Fr 2 h > 1 implies the possibility of local shear instability. The quantity Re h Fr 2 h can also be related to the buoyancy Reynolds number, denoted by Re b , via the inviscid estimated dissipation scaling, ε k ∼ u 3 h /L Hk . The Re b can be expressed as a ratio of inertial and viscous forces, i.e. the Reynolds number of Ozmidov-scale eddies. Eddies with size less than the Ozmidov length, l o , are not restrained from overturning by buoyancy. Since Re b = (l o /η) 4/3 , where η = (ν 3 /ε k ) 1/4 is the Kolmogorov scale, it follows that when Re b 1 there is a large scale separation between the Ozmidov and the Kolmogorov scales, which is necessary for the presence of an inertial subrange that is unaffected by either buoyancy or viscous forces. The turbulent horizontal Froude number, Fr h = u h /NL Hk , is the ratio of the buoyancy time scale (N −1 ) to the turbulent horizontal time scale (L Hk /u h ), and is a measure of the strength of buoyancy effects on the energy-containing scales; a decrease in Fr h implies an increase in the relative strength of buoyancy. The Froude number Fr h can also be interpreted based on length scales using ε k ∼ u 3 h /L Hk , which leads to Fr h = (l o /L Hk ) 2/3 . The Fr h -Re h Fr 2 h phase portrait serves to distinguish among the different manifestations of the influence of buoyancy as the wake evolves through progressively decreasing values of both Fr h and Re h Fr 2 h . As Fr h decreases, the lower bound (l o ) on the length scale of buoyancy-affected motions decreases. As Re h Fr 2 h decreases, the scale separation between l o and η decreases so that, for a given Fr h , the parameter Re h Fr 2 h is a measure of the range of scales that can support turbulent motions. The wake is found to pass from a state of weak buoyancy to stratified turbulence at Fr h = 1 when the large-eddy length scale, L Hk , becomes equal to the Ozmidov length scale, l o = (ε k /N 3 ) 1/2 . Figure 11(a) shows that Fr h = 1 is achieved at approximately the first buoyancy adjustment period (Nt b = 1) for all the simulated wakes. As was discussed in § 5, U 0 /U 0∞ is found to deviate sharply from unity at Nt b ≈ 1 marking a point when buoyancy has begun to significantly affect the largest scales of the flow.
Stratified wake turbulence can be further subcategorized into three regimes: weakly stratified turbulence (WST), intermediately stratified turbulence (IST) and strongly stratified turbulence (SST) as marked in figure 10. In the WST stage, the effect of buoyancy on the mean flow is significant but its effect on turbulence is not. In particular, turbulence anisotropy is hardly affected in the WST regime. The value of Fr h has to decrease from unity by almost an order of magnitude before there is a trend of increasing turbulence anisotropy associated with the r.m.s. in the horizontal progressively becoming larger than in the vertical. As discussed in § 5, turbulence anisotropy increases in both Fr = 2 and 10 wakes at Fr h ≈ 0.1. This suggests that WST transitions at Fr h ∼ O(0.1) to a regime of IST that is distinguished by progressively increasing turbulence anisotropy. The final stage of SST is based on Fr h becoming Fr h ∼ O(0.01). In particular, we consider SST to commence at Fr h = 0.03, based on the value of Fr v = u h /Nl v approaching an O(1) constant. It is worth noting that Fr h = 0.03 is close to the prediction of Lindborg (2006) that the critical horizontal Froude number Fr h,crit = 0.02. Figure 10 shows that the entry of stratified wakes into the three stratified turbulence regimes depends on the Fr of the wake. The F50 trajectory enters WST relatively late and is not able to enter IST. Notice that Fr h is, however, close to 0.1 (near the entry of IST) by the end of the computational domain and turbulence starts becoming anisotropic as seen in figure 6(a). Both F10 and F02 trajectories access the IST regime with a sufficiently large Re h Fr 2 h ≈ 100. Turbulence anisotropy progressively increases as the wake traverses the IST regime, as was shown by the divergence of u z from K 1/2 in figures 7(a) and 8(a). In the IST regime, the inertial force on the mean flow dynamically adjusts to balance the buoyancy force such that Fr V becomes close to an O(1) constant (see figures 9b and 9c). However, buoyancy becomes progressively stronger in the case of turbulent eddies, as indicated by the continued decrease of Fr v as the flow evolves. The Fr = 2 wake is able to cross the Fr h = 0.03 boundary to access the SST regime and remains in the regime until the end of the simulation domain. Figure 9(c) shows that Fr v = u h /Nl v asymptotes to an O(1) constant in the SST regime. Thus, stratification imposes a characteristic buoyancy length scale, l v = u h /N, on small-scale turbulence in the SST regime. Since l v = u h /∂ z u h , it follows that vertical shear of the fluctuations asymptotes to O(N). In case F02 (figure 12) l v becomes approximately constant when Nt b > 20. The value of Nt b = 20 also marks the entry of case F02 into the SST regime as indicated by the crossing of Fr h 0.03 at Nt b ≈ 20 by the F02 trajectory in figure 11.
Since Re h Fr 2 h ∼ O(10) as the F02 wake enters into the SST regime, figure 13(a) reveals small-scale patchiness of turbulence displayed by a contour of instantaneous turbulent dissipation. While F02 turbulence has decayed in amplitude during the SST regime, figure 13(b) shows that the power spectra of centreline streamwise velocity remains broadband during the SST regime. Figure 14 displays the three regimes of stratified turbulence. At F02, the SST spans a substantial streamwise distance (almost 85 body diameters from x/L b ≈ 40 to the end of the computational domain) indicating that the dynamics of strongly stratified flow is important to the evolution of high-Re wakes.

Wake length scales
The evolution of the geometrical dimensions of the wake is closely related to the decay of U 0 owing to conservation of momentum deficit. The horizontal (L H ) and vertical (L V ) wake extents derived from profiles of mean streamwise velocity deficit  are shown in figure 15. The width and height of the unstratified case are identical by construction, as they are obtained from additional azimuthal averaging. It is found that, during the interval 10 < x/L b < 65, the wake satisfies L ∝ x 0.45 in accord with momentum-deficit conservation, U 0 L 2 = const., and the U 0 ∝ x −0.9 law found in that interval. During the interval 65 < x/L b < 125, the wake length grows with an exponent L ∝ x 0.29 , close to L ∝ x 1/3 , the classical similarity law for L in the turbulent far wake.
Consider the wake length L k , derived from profiles of TKE in the UNS case, and plotted in figure 16 (solid black line). Unlike L, which exhibits a breakpoint at x/L b ≈ 65, there is no corresponding breakpoint in the evolution of L k . Rather, a least-squares power-law fit to L k yields L k ∝ x 0.31 , close to x 1/3 , during the entire domain after x/L b > 10. Therefore, it is a combination of quantities derived from turbulence (K 1/2 and L k ) and not those derived from the mean flow (U 0 and L) that appears to satisfy the classical high-Re similarity result. We return to figure 15 and discuss how buoyancy sets in to create anisotropy between L H and L V in stratified wakes. For the weakly stratified Fr = 50 and 10 cases, both L H and L V behave in the near wake as expected with little to no variation compared to the UNS case. Further downstream, the growth of L V reduces at Nt b ≈ 1 (recall that Nt b = Fr −1 x/L b ). As far as we know, this is the first body-inclusive stratified wake simulation that finds a continuous decrease of vertical thickness. This continuous decrease in vertical length scale is consistent with the strongly stratified regime of turbulence decay, which was theoretically derived by Davidson (2010). Also, the finding is consistent with Billant & Chomaz (2001), who state that the vertical length scale is O(u h /N) in a strongly stratified flow; here, u h denotes the r.m.s. of horizontal velocity fluctuation. The horizontal width L H , for both Fr = 10 and 50, deviates away from the UNS case before L V does. Beyond x/L b = 10 in the Fr = 10 wake, the growth of L H is close to L H ∝ x 1/3 over a long interval until the end of the computational domain.
The wake dimensions in the Fr = 2 case are significantly different from the UNS case. Both L H and L V start to deviate from the UNS case close behind the body, indicating an absence of the '3-D, unstratified' regime. The deviation is consistent with conservation of momentum deficit; L H is enhanced to counteract the contraction of L V as buoyancy returns vertically displaced fluid towards its neutral equilibrium. After x/L b = 2, which marks the end of the recirculation zone, L H commences to increase with a rate slightly larger than that of the UNS case until x/L b = 5 or Nt b = 2.5, where the growth rate decreases. At x/L b ≈ 2π, where the mean wake enters into the next phase of the oscillatory modulation depicted by the rapid increment of L V , the velocity deficit (U 0 ) increases and, as a consequence of conservation of momentum deficit, L H decreases sharply to satisfy U 0 L H L V ≈ const. This sharp lateral contraction of L H leads to the F02 trajectory crossing its unstratified counterpart at x/L b = 10. Beyond x/L b = 20, the growth rate of L H progressively approaches L H ∝ x 1/3 . After the rapid growth of L V over 2π < x/L b < 4π (π < Nt b < 2π), the oscillatory modulation remains visible, albeit with diminished amplitude, throughout the wake evolution until the end of the computational domain.
The Fr = 2 power-law exponents of L H and L V can, in fact, be inferred based on Fr V → const. and the conservation of momentum. The proportionality Fr V ∝ x 0 implies that L V ∝ U 0 . Invoking the conservation of streamwise momentum deficit, where U 0 L V L H ≈ const. and assuming that the shape of deficit profiles in the y-z plane is characterized by the two length scales L H and L V , yields L H L 2 V ≈ const. or L H ∝ L −2 V . Since U 0 ∝ x −0.18 in the Fr = 2 wake, it directly follows that L V ∝ x −0.18 and L H ∝ x 0.36 , which are close to the present result.
The horizontal extent (figure 16a) of the TKE profile becomes comparable to that of the UNS case, even in the NEQ regime, with a grow rate of L Hk ∝ x 1/3 . In the vertical direction, the half-height (L Vk ) based on TKE decreases similar to L V . For example, L Vk at F10 begins to deviate away from the UNS case around x/L b = 10 (Nt b = 1), practically identical to where L V deviates from its UNS counterpart. In case F02, once wake turbulence enters into the NEQ stage, it is found that L Vk ∝ K 1/2 (∝ x −0.18 ) in agreement with the predicted vertical length scaling of strongly stratified flow.

Histories of kinetic and potential energies
In § 5, we discussed buoyancy effects on centreline values of mean and r.m.s. velocity. The evolution of area-integrated values of the kinetic and potential energy, each split into mean and turbulent constituents, is also of interest. Figure 17 shows the histories of area integrals, denoted by {·}, of MKE, E M k = u d 2 /2, TKE, E T k = u i u i /2, mean available potential energy (APE), E M ρ = γ ρ d 2 /2, and turbulent APE, E T ρ = γ ρ 2 /2, where γ = g 2 ρ −2 0 N −2 . The subscript 'd' denotes deviation of velocity from background velocity (U b , 0, 0) or density from the linearly varying background density.
Consider first the decay of E M k (figure 17a). The decay of E M k in the stratified wakes slows down once the downstream location in each wake has reached the corresponding Nt b = 1. The inset shows that E M k , when normalized by the corresponding UNS value, increases as approximately x 0.5 , which can be explained as follows. If U 0 , L H and L V are the characteristic velocity and lateral wake length scales for the mean deficit velocity profile, then (8.1) The second proportionality comes from invoking conservation of deficit momentum and, again, assuming that the shape of the deficit profile is characterized by L H and 18 /x −2/3 ∝ x 0.49 , which is close to the x 0.5 scaling in the inset. Similarly, scaling E T k shows that K. Chongsiripinyo and S. Sarkar In the WST regime, the evolution of K is unchanged while U 0 increases relative to the unstratified case, so that (8.2) implies that {E T k } decreases relative to UNS once the stratified wakes enter the WST regime. This can be seen in figure 17(b), where the Fr = 10 and 50 wakes access the WST regime at their corresponding Nt b ≈ 1, at which point they decrease faster than the UNS case. The inset shows that the reduction of the {E T k }/{E M k } ratio relative to the unstratified case largely takes place in the WST and the early-IST regimes and reaches a plateau in the SST regime. In the unstratified case, the mean component becomes much smaller than the turbulent component. However, as inferred from the inset of figure 17(b), the mean wake remains at long distance from the body in stratified cases even when the turbulence has decayed.
Consider the evolution of potential energies in figure 17(c,d) figure 17(d) shows that the TKE and turbulent potential energy are more evenly partitioned with {E T ρ }/{E T k } ≈ O(1) over a substantial distance in the Fr = 2 and Fr = 10 wakes. This is a consequence of hydrostatic balance of vertical motion in buoyancydominated high-Re flow, where the vertical pressure gradient is balanced by buoyancy, giving ρ ∼ u 2 h ρ 0 /gl v and, since l v ∼ u h /N, it follows that ρ 2 ∝ (u h N) 2 or E T ρ ∼ E T K .

Budgets of mean and turbulent kinetic energies
The velocities, length scales and energies discussed in the previous sections are statistical representations of mean flow and turbulence whose evolution is modified in stratified wakes. We further elaborate on the dynamical processes responsible for these modifications in stratified wakes by means of the evolution equations for MKE and TKE. The budget of MKE is written as where the advection, production, dissipation, transport and buoyancy fluxes are defined as follows: Similarly, the budget of TKE is written as where the advection, dissipation, transport and buoyancy fluxes are defined as follows: The only additional term arising from the subgrid viscosity is u i ∂ j τ s ij , where τ s ij = 2ν s S ij with S ij = (∂ j u i + ∂ i u j )/2. The contribution of this term is found to be negligible.
We start by contrasting area-integrated kinetic energy budgets between the Fr = 2 and Fr = ∞ wakes. Terms in the MKE budget, (9.1), are plotted in figure 18(a,c). The figure shows that the magnitude of {P k } {D M k } (recall that {·} denotes an integral over cross-sectional area in the y-z plane), indicating that the production of TKE is the primary sink of MKE. This is another manifestation of high-Re behaviour of this flow; the close to x −1 scaling of U 0 in the UNS case cannot be attributed to a low-Re similarity result. The large pressure transport, {T M k,p } in the insets, is a consequence of the relaxation of the flow from the low-pressure recirculation zone with increasing streamwise distance. While the Fr = ∞ pressure transport decreases significantly and becomes virtually negligible after the relaxation, the Fr = 2 pressure transport manifests wave-like behaviour with a wavelength of gravity lee waves, λ ≈ 2πFr, and progressively decreasing amplitude. This oscillatory modulation is observed in the behaviour of every term in the MKE budget except the dissipation. The oscillations of the advection, {A M k }, and the pressure transport have opposing phase, and the wavelength of the mean buoyancy flux, {B M k }, is λ ≈ πFr. Terms in the TKE budget, equation (9.2), are plotted in figure 18(b,d). Turbulence is established early on, as shown by the large positive production, {P k }, that exceeds the dissipation, {D T k }, in the region before x/L b ≈ 5. While the decay of TKE in the 885 A31-22 K. Chongsiripinyo and S. Sarkar FIGURE 18. Terms in the MKE budget (a,c) and in the TKE budget (b,d) for the Fr = ∞ wake (a,b) and the Fr = 2 wake (c,d). The budget terms are normalized with U 3 b L b . The insets show budgets in the near field (0.5 x/L b 10) while each main panel depicts the budgets in 10 x/L b 100. The insets also separate the transport into pressure and advection transports. The remaining transport terms are negligible. Fr = ∞ wake via {D T k } is counteracted solely by {P k }, the turbulent buoyancy flux, {B T k }, plays a role in stratified wakes. In particular, {B T k } is found to be significant in the Fr = 2 wake after x/L b = 5. The term {B T k } has an order of magnitude that is comparable to the production and the dissipation, and it is in phase with the advection term.
To help understand how MKE and TKE evolve, the production, buoyancy flux and dissipation are normalized by the Lagrangian rate of change of their corresponding energy, taken to be δ t E k ∼ O(E k U b /x), and plotted in figure 19. It is clear that all the simulated wakes are, again, of the high-Re type until the end of the domain since the compensated production, {P k } M = {P k }/{δ t E M k }, is a much larger sink for MKE than its dissipation counterpart, {D M k } M . The inset in figure 19(a) shows that the production is suppressed in the presence of stratification after Nt b ≈ 1 but not necessarily for the entire wake. Thus, {P k }/{P k } ∞ is <1 for the Fr = 50 and 10 wakes but becomes >1 in the Fr = 2 wake beyond approximately x/L b = 60. However, regardless of the fact that the Fr = 2 production is promoted later in the Fr = 2 wake ({P k } F02 /{P k } ∞ > 1), the compensated production for all the stratified wakes remains relatively small at all times after Nt b ≈ 1 ({P k } M < {P k } M∞ ). Figure 19(c) shows that the mean buoyancy flux takes positive values, indicating that energy is transferred from mean APE to x/L b 10 1 10 0 10 2 FIGURE 19. Ratio of cross-wake integrated budget terms for the evolution equation for E M K (a,c,e) and E T K (b,d,f ). The budget terms are normalized by the Lagrangian rate of change of their corresponding energy, MKE. Observe the Fr = 2 wake; the oscillatory modulation induces a huge jump in the mean buoyancy flux at x/L b ≈ 2π, the same location where the production is most suppressed. This relatively large B M k where {B M k } M | peak ∼ O(1) implies that B M k is an important contributor to E M k being larger in stratified wakes. Therefore, reduction in the decay rate of E M k in the stratified Fr = 2 wake has two reasons: (i) the early suppression of production, and (ii) the significant positive mean buoyancy flux.
TKE is reduced in stratified wakes, as was seen in figure 17(b), and the reasons can be deduced from figure 19(b,d,f ). Interestingly, it is primarily the turbulent buoyancy flux (B T k ) and, additionally, the dissipation of TKE (D T k ) that cause the initial reduction of E T k , and it is only later that the direct effect of reduced production sets in to further reduce the TKE. This is seen from the following points: In other words, turbulence in the Fr = O(1) wakes is energized by the potential energy that was injected during the intermediate-wake period and not only the enhanced production. Later (after x/L b ≈ 70), although {B T k } T has become negative, the significant increased {P k } T is sufficient to counteract {D T k } T so that the reduced decay rate of {E T k } remains intact. It is worth noting that (consider the Fr = 2 and Fr = 10 cases) if the NEQ regime is characterized solely by the reduced decay rate of {E M k } in figure 17(a), the present result of partially positive B T k will be in contrast with the finding of Redford et al. (2015) for weakly stratified wakes that B T k is a sink of E T k throughout the NEQ regime. However, as described in § 5, the official arrival of the NEQ regime should be based on the reduced decay rate of {E T k } in figure 17(b) at Nt b ≈ 20 where the wake enters into the SST regime. If we accept Nt b ≈ 20 as the official arrival of the wake (both mean and turbulence) into NEQ, then B T k behaves as a sink in the NEQ regime in agreement with Redford et al. (2015).

Dissipation
The evolution of the rate of dissipation (ε), characterized by centreline values, is discussed and possible scaling laws are explored in this section. The rates of dissipation (ε k ) of TKE in lines and TPE (ε ρ ) in symbols are shown in figure 20. Given that the unstratified wake has high Re, classical theory suggests the scaling ε ∝ x −7/3 , which is a direct result from substituting the self-similarity solution u ∝ x −2/3 and l ∝ x 1/3 into the inviscid estimate ε ∼ u 3 /l. To test the preceding inviscid scaling of ε, we take u = K 1/2 and l = L k to be the half-width of the TKE profile. The power-law exponent, m, as in ε ∝ x m , is obtained by least-squares fits to the simulation data for ε. For the unstratified wake, we find m ≈ −2.25 followed by m ≈ −2 (after a slope change at x/L b ≈ 50). Thus, the ε power law in the unstratified case is not consistent with the u 3 /l scaling prediction of x −2.44 that follows from the power-law exponents (u = K 1/2 ∝ x −0.71 and l = L k ∝ x 0.31 ) found in the simulations as discussed in § 5. Figure 21(a) shows the evolution of ε using various choices for the normalization. Note that, interestingly, ε(U 3 0 /L) −1 ≈ const. beyond x/L b ≈ 65 in figure 21(a). Figure 21(a) also emphasizes the point that ε evolves differently from u 3 /l by showing that the scaled dissipation, C ε = ε(u 3 /l) −1 (in green triangles), is not constant but increases with streamwise distance in the wake. This apparent contradiction between ε and u 3 /l was previously found in an experiment on turbulence generated by a fractal grid (Seoud & Vassilicos 2007), where C ε = const. despite an observed −5/3 inertial range in the energy spectrum. Vassilicos (2015) points to accumulating evidence in some turbulent shear flows that C ε = Re m Global /Re n Local ( =const.) in a so-called 'non-equilibrium region' wherein the rate at which the TKE cascades from large-to small-scale motion, u 3 /l (formed by a dimensional argument), is not strictly equal to the dissipation, ε, felt by the TKE. Our result shown in figure 21(b) presents additional support that the centreline dissipation can reasonably be scaled by the form (Re n Global /Re n Local )K 3/2 /L, where Re Local = U 0 L/ν at least beyond x/L b = 40. In stratified flows, the dissipation of TKE and turbulent potential energy have often been suggested to scale using the inviscid estimate, with u and l replaced by u h and l h , i.e. ε k,ρ ∼ u 3 h /l h . This comes from the notion that the nonlinear forward energy cascade is driven largely from the fluctuating energy present in horizontal motions. However, if the non-equilibrium region exists in the stratified wakes, we should expect that ε k,ρ = C ε k ,ε ρ u 3 h /l h , where C ε k ,ε ρ = const. While C ε is found to be close to a constant in the homogeneous stratified turbulence of Maffioli & Davidson (2016), assessing this hypothesis can be difficult for localized turbulence such as in a core of the stratified wake where the lack of a homogeneous direction presents a challenge for the estimation of l h . If ε k ∼ u 3 h /l h , l h ∼ L Hk and centreline estimates are employed, the decay rate of ε k should reduce once the stratified wake has entered into the IST regime where the decay rate of u h is reduced. This is seen in figure 20 for the Fr = 2 wake but not for the Fr = 10 case. In the Fr = 2 and 10 wakes, ε ρ appears to decay like its ε k counterpart but only after Nt b 1, i.e. after the initial adjustment between the APE of the displaced fluid and its KE is completed. Scaling of the Fr = 2 dissipation of TKE is presented in figure 22(a). Similar to the Fr = ∞ wake, ε k ∼ (U 3 0 /L) after approximately x/L b ≈ 65. The alternative scaling, ε k ∼ u 2 h N, appears valid as the wake approaches the SST regime in which there is turbulence (Re b O(1)) and it is largely controlled by buoyancy. As the Fr = 2 wake approaches towards the viscous regime, our result (not shown) reveals that ε k becomes increasingly dominated by the components associated with the vertical shear of horizontal motions. This implies that the scaling ε k ∼ νu 2 h /l 2 v , deduced from a dimensional argument, should become applicable. Indeed, ε k (νu 2 h /l 2 v ) −1 flattens in the downstream direction. The non-equilibrium dissipation scaling is shown in figure 22(b). Unlike the unstratified wake, there is no evidence that the non-equilibrium dissipation scaling works for the simulated Fr = 2 wake.

Summary and conclusions
The present work is devoted to an examination of relatively high-Reynolds-number turbulent wakes that evolve in stratified fluids. Body-inclusive LES of unstratified and stratified flows past a circular thin disk placed perpendicular to a background stream are conducted. The wakes are simulated at a Reynolds number of Re = U b L b /ν = 50 000 and over a wide range of stratifications: Fr = U b /NL b = ∞, 50, 10 and 2.
In the axisymmetric wake at Fr = ∞, it is found that the centreline deficit velocity (U 0 ) does not necessarily scale with the turbulence velocity scale (K 1/2 ), contrary to classical theory. The mean flow decays in two stages with a break in slopes at x/L b ≈ 65, but no such transition is found for turbulence. The first stage of 10 < x/L b < 65 exhibits U 0 ∝ x −0.9 and L ∝ x 0.45 . It is only after x/L b ≈ 65 that U 0 exhibits a power law that is close to the high-Re scaling U 0 ∼ K 1/2 ∝ x −2/3 and L ∼ L k ∝ x 1/3 . As a result, none of the stratified cases considered here display the U 0 ∝ x −2/3 scaling in the 3-D regime before buoyancy effects set in. It is worth noting that the x −0.9 behaviour of U 0 in the intermediate stage of the disk wake is similar to the approximately x −1 behaviour found in the intermediate stage of sphere wakes at Re = 3700 by Pal et al. (2017), Re = 10 000 by Chongsiripinyo et al. (2019), and fractal-plate wakes at Re = 5000 (DNS) and Re = 50 000 (laboratory experiment) by Dairay et al. (2015).
Stratified wakes evolve with scant effect of buoyancy until approximately one buoyancy time scale (Nt b ≈ 1). At this point, the vertical mean Froude number decreases to and remains at Fr V ∼ O(1) and the stratified wake then enters into a stratified turbulence regime, which can be further subcategorized into three regimes: WST, IST and SST. We use a turbulence-based classification, unlike the mean-based categorization from Spedding (1997), to delineate the regime boundaries by the turbulent horizontal Froude number Fr h = u h /NL Hk ; here, u h and L Hk are r.m.s. horizontal velocity fluctuations and TKE-based horizontal wake width. At the entry of the WST stage, which is marked by Fr h = 1, the large-eddy length scale (L Hk ) becomes equal to the Ozmidov length scale (l o = (ε k /N 3 ) 1/2 ). As the wake progresses in the WST stage, Fr h reduces from O(1) to O(0.1) and, while the effect of buoyancy on the mean flow is significant, its effect on turbulence intensities is insignificant. Thus, during WST that spans 1 Nt b 5, the mean flow progressively transitions into the NEQ regime but turbulence remains approximately isotropic (u h ≈ u z ). A transition towards the next stage of IST takes place at Fr h ∼ O(0.1) in the proximity of Nt b ≈ 5. The IST stage is distinguished by progressively increasing turbulence anisotropy. In IST, while the mean flow has arrived at the NEQ-regime behaviour of U 0 ∝ x −0.18 , turbulence is still in transition towards a different regime, namely, SST.
The stratified wake enters into the third stage of SST, demarcated by the value of Fr v asymptoting to an O(1) value. We find that, at Fr h ≈ 0.03, which occurs at Nt b ≈ 20, the wake enters the SST regime with Fr v → 0.65. Indeed, a key feature of the SST regime of stratified flows in general is that buoyancy imposes a vertical length scale (l v ∼ u h /N), which follows from the result Fr v → O(1) const. (Billant & Chomaz 2001), as has been seen in previous simulations of decaying homogeneous turbulence (Maffioli & Davidson 2016;de Bruyn Kops & Riley 2019). It is expected that SST lasts as long as the buoyancy-modified Reynolds number Re h Fr 2 h > 1 ( Riley & de Bruyn Kops 2003;Brethouwer et al. 2007), which remains true in the F02 wake from x/D ≈ 40 to the end of the domain at x/D = 125. In the SST regime of the wake, turbulence is strongly anisotropic (u z u h ) and arrives into the NEQ regime with u h ∼ U 0 ∝ x −0.18 while u z ∝ x −1 . Furthermore, the SST wake is patchy, as inferred from visualizations of the fluctuating dissipation rate, and the temporal spectra are broadband.
Throughout the IST-SST period, the vertical mean length scale (wake height) is naturally selected by the flow so that mean inertial and buoyancy forces remain in balance. Therefore, in the IST-SST period, L V ∼ U 0 /N, which yields Fr V ≈ const. The constraint of Fr V ≈ const. leads to the continuous decrease in L V to accommodate the ongoing U 0 decay. The need to obey conservation of momentum (U 0 L H L V ≈ const.) indicates that the wake width (L H ) is thus indirectly selected by buoyancy as L H ∝ U −1 0 L −1 V ∝ U −2 0 . It is worth noting that L V exhibits a decreasing trend throughout the computational domain, indicating that, unlike prior low-Re studies of stratified wakes, the effect of viscous diffusion to increase L V is not dominant in the present simulations at comparable Nt.
It is worth noting that, in the present work, it is L k and not L that possesses the power law close to the x 1/3 growth in the unstratified wake and the 3-D region of stratified wakes. It is intriguing that the results for power laws satisfied by the mean velocity deficit and horizontal wake width from temporal-model simulations (e.g. Dommermuth et al. 2002;Diamessis et al. 2011) are similar to power laws in the present work that are based on fluctuation ( √ K) and not mean profiles. A possible reason for the difference in the power laws for U 0 between temporal and body-inclusive simulations is the consequence of a large-scale structure introduced by the body into the wake that is not present in the initial conditions of the temporal simulations.
The present results suggest re-examination of the NEQ regime of stratified wakes. At Nt b ≈ 1, which has been taken to be the commencement of the NEQ wake, it is only the mean flow that is sufficiently influenced by buoyancy. In particular, at Nt b ≈ 1, the decay of U 0 is slowed down and the mean-based Froude number (Fr V ) becomes O(1). In other words, L V becomes O(U 0 /N). It is much further downstream at Nt b ≈ 20, when the turbulence-based Froude number (Fr v ) becomes O(1). At this point, turbulent fluctuations experience strong buoyancy effects: (1) turbulence becomes strongly anisotropic with u z u h , and (2) the decay of turbulence is also reduced so that u h ∼ U 0 .
The U 0 ∝ x −0.18 found for the NEQ disk wake is different from the previously found sphere-wake behaviour of U 0 ∝ x −0.25 in the NEQ regime. Further study will be useful to ascertain how the geometry of the wake generator affects the exponent of the U 0 power law in the NEQ region of stratified wakes.
An overall picture of stratified wakes emerges when the results are plotted in phase space. We choose centreline values of (Re h Fr 2 h , Fr h ) as the phase-space coordinates, which is analogous to the study of stratified turbulence in a periodic box by de Bruyn Kops & Riley (2019), who use (Re b , Fr h ), where Fr h is defined using the horizontal r.m.s. velocity and a horizontal integral length scale. The Fr = 10 and 2 wakes at Re = 50 000 come close in phase space once they enter the regime of stratified turbulence at Fr h = 1. Therefore, the Fr = 10 wake can be expected to access the SST regime at Fr h = O(0.01) similar to what is found for the Fr = 2 wake. The SST regime of the Fr = 2 wake commences at x/L b ≈ 40 and continues until the end of the domain, which is approximately 85 body diameters away. For a given body Fr, the span of the SST regime is expected to increase with increasing body Re of the wake.
Stratification prolongs the lifetime of the wake. The mechanisms leading to the preservation of wake MKE originate in the WST/early-IST stages, where the NEQ regime of the mean flow is being established. During the early-NEQ stage, the mean buoyancy flux is positive, so that energy is transferred from APE to MKE. This feature was not observed in our prior temporal-model simulations since the temporal model does not capture the correct potential energy change as the fluid goes past the body. Additionally, production of TKE that acts as a sink of MKE is reduced. In the late-IST/SST stages, the decay rate of MKE remains relatively smaller as the production in the late wake remains small in normalized form. The decay rate of TKE is initially increased in the WST/early-IST stages primarily because of negative buoyancy flux and increased dissipation. It is only later that the direct effect of reduced production sets in to further reduce the TKE. In the late-IST/early-SST stages, not only is production enhanced but also turbulence is energized by the injection from turbulent potential energy. Consequently, the TKE decay slows down similar to the buoyancy-induced slowdown of U 0 that had occurred closer to the body. As mean and turbulent buoyancy fluxes are positive during the establishment of the mean and turbulence NEQ regimes, respectively, it is clear that stratification has a direct and positive effect on the prolongation of wake life. Only in the SST stage does the turbulent buoyancy flux become negative, acting as a sink of TKE and to mix the buoyancy field.
The lee-wave-induced 'oscillatory modulation' of U 0 reported in the DNS of the Re = 3700 sphere wake by Pal et al. (2017) is a manifestation of the transfer between APE and MKE in the NEQ regime. Oscillatory modulation is also found in the present simulations, showing its importance in the case of a disk as the wake generator and at an order of magnitude higher Re = 50 000. It is worth noting that the lee-wave-induced oscillatory modulation is absent in temporal models.
In the unstratified wake, it is found that the decay of ε is not consistent with the classical K 3/2 /L k scaling but is described better with the NEQ scaling (Vassilicos 2015) of the form ε = C ε (K 3/2 /L k ) with C ε = Re n Global /Re n Local ( = const.). In the Fr = 2 wake, the scaling ε k ∼ u 2 h N appears valid as the wake approaches the SST regime in which there is turbulence (Re h Fr 2 h O(1)) but it is largely controlled by buoyancy. Further research is required to explore the SST regime of the stratified turbulent wake. Although the SST regime is recognized in the Fr = 2 wake simulated here, it will be beneficial to further increase the buoyancy Reynolds number. A follow-up study of disk wakes is planned at a larger Re and a smaller Fr, a change that would increase Nt for the same x/L b and Re h Fr 2 h for the same value of Nt.