The high-Reynolds-number stratified wake of a slender body and its comparison with a bluff-body wake

The high-Reynolds number stratified wake of a slender body is studied using a high-resolution hybrid simulation. The wake generator is a 6:1 prolate spheroid with a tripped boundary layer, the diameter-based body Reynolds number is $Re= U_\infty D/\nu = 10^5$, and the body Froude numbers are $Fr=U_\infty/ND=\{2,10,\infty\}$. The wake defect velocity ($U_d$) decays following three stages with different wake decay rates \citep{Spedding1997} as for a bluff body. However, the transition points among stages do not follow the expected $Nt = Nx/U_\infty$ values. Comparison with the wake of a circular disk in similar conditions \citep{Chongsiripinyo2020} quantifies the influence of the wake generator - bluff versus slender - in stratified flow. The strongly stratified $Fr=2$ wake is in a resonant state. The steady lee waves strongly modulate the mean flow and, relative to the disk, the $Fr=2$ spheroid wake shows an earlier transition from the non-equilibrium (NEQ) stage to the quasi-two-dimensional (Q2D) stage. The NEQ-Q2D transition is followed by a sharp increase in the turbulent kinetic energy and horizontal wake meanders. At $Fr=10$, the start of the NEQ stage is delayed. Transfers between kinetic energy and potential energy reservoirs (both mean and turbulence) are analyzed and the flows are compared in phase space (local Froude and Reynolds number as coordinates). Overall, the results of this study point to the difficulty of finding a universal framework for stratified wake evolution, independent of the features of the body, and provide insights into how buoyancy effects depend on the wake generator.


Introduction
Due to their low drag coefficients, slender bodies are extensively used in aerospace and naval applications. Multiple studies have described the flow around these bodies focusing on the drag force, the boundary layer, and the flow separation (Wang 1970;Costis et al. 1989;Wang et al. 1990;Chesnakas & Simpson 1994;Fu et al. 1994;Constantinescu et al. 2002;Wikström et al. 2004). However, despite their presence in many underwater applications, only a few works have looked into the wake of a slender body (Chevray 1968;Jiménez et al. 2010;Kumar & Mahesh 2018) and, only recently, the far wake of a slender body has been studied (Ortiz-Tarin et al. 2021).
The near wake of a slender body with a turbulent boundary layer (TBL) is characterized 2 Ortiz, Nidhan and Sarkar by having a small recirculation region. The recirculation region is surrounded by a ring of small scale turbulence that emerges from the boundary layer and does not show strong vortex shedding (Jiménez et al. 2010;Posa & Balaras 2016;Kumar & Mahesh 2018;Ortiz-Tarin et al. 2021). As a result, the wake is thin and develops slowly compared to the wake of bluff bodies. These particular features of the slender body high-Re near wake lead to interesting effects further downstream: (i) despite having a smaller drag coefficient than bluff bodies, the defect velocity ( = ∞ − ) of the slender body wake can be larger than that of a bluff body for a long downstream distance, (ii) the turbulent kinetic energy of the wake shows an off-center radial peak at the location where the turbulent boundary layer separates -instead of a Gaussian profile with a central peak, and (iii) helical instabilities come into play only in the intermediate and far field of the wake. These particularities affect the scaling laws of the wake. In a domain spanning 80 the defect velocity, the kinetic energy, and the dissipation do not follow the classic high-Re scaling and they decay differently than bluff body wakes (Ortiz-Tarin et al. 2021) exhibiting a non-equilibrium scaling of dissipation (Vassilicos 2015;Dairay et al. 2015).
The few studies that look into slender body wakes assume that the body moves in an unstratified environment, where the density of the surrounding fluid is constant. However, in a realistic underwater marine environment the effect of density stratification due to salinity and temperature can become relevant. Density stratification suppresses vertical motions, triggers the formation and sustenance of coherent structures, and leads to the radiation of internal gravity waves. More importantly, in a stratified environment, the wake of a submersible lives longer than in an unstratified environment, i.e., it takes more time for the flow disturbance to die out (Spedding 2014). The study of stratified wakes has been nearly exclusively focused on the flow past bluff bodies (Lin & Pao 1979;Hanazaki 1988;Lin et al. 1992;Chomaz et al. 1992;Orr et al. 2015;Pal et al. 2017) and underwater topography (Drazin 1961;Castro et al. 1983;Baines 1998). Here, we study the influence of stratification on the high-Re wake of a prolate spheroid with a turbulent boundary layer.
The strength of ambient stratification is measured by the body-based Froude number Fr = ∞ / . This is the ratio between the convective frequency of the flow, ∞ /where ∞ is the freestream velocity and the diameter of the body -and the buoyancy frequency . In the wake of ocean submersibles, Fr ∼ (1 − 10 2 ). However, since the velocity deficit in the wake ( ) decays with the streamwise distance and the wake width ( ) increases, the Froude number defined with local variables Fr = / decreases as the flow evolves. Thus, even in a weakly stratified environment, eventually all wakes are affected by stratification.
Since the relative strength of stratification increases locally as the flow develops, the evolution of the stratified wake is multistage in nature. Based on the measurements of and in high-Fr (i.e. initially weak stratification) bluff body wakes, Spedding (1997) identified three regimes in stratified wake evolution based on the power-law decay rates of . These regimes are generally identified by empirically fitting decay rates to in different temporal (or spatial) regions and are as follows: (i) Three-dimensional (3D) regime: Close to the generator, wake decay is similar to the unstratified wake of the corresponding body shape. This is the so-called 3D regime and lasts until the buoyancy time defined by = / = / · 1/Fr approaches (1), equivalently until / ∼ Fr.
(ii) Non-equilibrium (NEQ) regime: As the wake evolves, buoyancy effects become progressively stronger. The decay of slows relative to the 3D regime and, furthermore, anisotropy between the vertical and horizontal velocity components increases. Spedding (1997)  ∼ −0.25±0.04 during the NEQ regime according to Spedding (1997). However, there has been some variability in the observed NEQ decay rate in later studies. Bonnier & Eiff (2002) reported a NEQ regime with ∼ −0.38 for 1.5 < Fr < 5. Brucker & Sarkar (2010) and Diamessis et al. (2011) reported ∼ −1/4 during the NEQ regime in their temporal simulations. Chongsiripinyo & Sarkar (2020) found that ∼ −0.18 during the NEQ regime of their Fr = 2 and 10 disk wakes. For wakes with Fr ∼ (1), the NEQ decay rate is preceded by a pronounced oscillatory modulation in which is linked to lee waves (Pal et al. 2017;Chongsiripinyo & Sarkar 2020;Ortiz-Tarin et al. 2019).
(iii) Quasi two-dimensional (Q2D) regime: After the NEQ regime, the stratified wake enters into the quasi two-dimensional regime (Q2D regime). The Q2D regime is characterized by transition in the power law to a significantly increased decay rate, e.g. Spedding (1997) reports a transition to ∼ −3/4 . In the Q2D regime, the wake progressively organizes into vortices that meander primarily in the horizontal plane (Gourlay et al. 2001;Dommermuth et al. 2002;Brucker & Sarkar 2010) and take the form of 'pancakes'. Although the wake motion in this regime is primarily in the horizontal plane, there is a variability in the vertical direction in the form of layers (Spedding 2002), hence the prefix 'quasi'.
In recent literature, stratified wakes have been characterized using turbulence features (Zhou & Diamessis 2019; Chongsiripinyo & Sarkar 2020) instead of the -based criteria of Spedding (1997). These studies are motivated by an attempt to connect buoyancy-related wake transitions to the broader stratified turbulence field.
Notice that the arrival of the wake into each of the three stages in its evolution depends on the value of , which is equivalent to a downstream distance of /Fr from the wake generator. At high Froude number, the downstream distance required to reach the NEQ and Q2D regions can become very large. Consequently, the size of the computational domain required to access these regimes rapidly becomes computationally unfeasible. To circumvent these limitation temporal simulations were used in the study of Gourlay et al. (2001). Temporal simulations use a reference frame moving with the wake where time correlates with streamwise distance in a fixed reference frame. By assuming that the streamwise development of the flow is slow, periodic boundary conditions can be used and the equations are advanced in time without the need of introducing the wake generator. This reduces the computational cost significantly. Most of the studies that have contributed to our current understanding of stratified wakes use temporal simulations (Gourlay et al. 2001;Dommermuth et al. 2002;Brucker & Sarkar 2010;Diamessis et al. 2011;De Stadler & Sarkar 2012;Abdilghanie & Diamessis 2013;Redford et al. 2015;Rowe et al. 2020).
The main drawback of the temporal model is the influence of its initialization. Since the flow at the wake generator is not solved, the starting profiles of the mean and turbulence have to be assumed. These simulations lack some specific features that are generated due to the body, e.g., steady lee waves, near-wake buoyancy effects, and the vortical structures shed from the boundary layer. Even when it is tempting to assume that body specific features are lost far from the body, the universality of the wake decay has remained elusive to experiments (Bevilaqua & Lykoudis 1978;Wygnanski et al. 1986;Redford et al. 2012), even in unstratified wakes. An alternative to temporal simulations are body inclusive simulations that retain the wake generator dependent features at the expense of a higher computational cost and a limited domain size (Orr et al. 2015;Chongsiripinyo et al. 2017;Pal et al. 2017;Nidhan et al. 2019;More et al. 2021).
To the best of authors' knowledge, Ortiz-Tarin et al. (2019) performed the first study of a stratified flow past a slender body that investigates the near and intermediate wake dynamics. Their analyses reveal that at Fr ∼ (1) the type of separation and the subsequent wake establishment is strongly dependent on the characteristic frequency of the lee waves and the aspect ratio of the body. When half the wavelength of the steady lee waves ( = 2 Fr) matches the length of the slender body ( ), the separation of the boundary layer is inhibited by buoyancy effects. Based on this condition, a critical Froude number can be defined Fr = / . When Fr > Fr stratification suppresses the generation of turbulence in the near wake, when Fr ≈ Fr buoyancy strongly limits the flow separation and can lead to a relaminarization of the wake at low Reynolds numbers. Finally, when Fr < Fr , the lee waves enlarge the separation region and there might be an increase in the turbulence intensities in the wake. When Fr ≈ Fr the wake is in a resonant state with both the separation and the wake dimensions being strongly modulated by the steady lee waves (Hunt & Snyder 1980;Chomaz et al. 1993;Ortiz-Tarin et al. 2019).
As mentioned before, the use of body inclusive simulations has one major limitation, i.e., the high computational cost. Due to the high resolution required to solve the boundary layer of the wake generator, the downstream domain is limited and thus the possibility of looking into the far wake gets significantly restricted, particularly at high Re. VanDine et al. (2018) presented a hybrid spatially-evolving model, which builds on the hybrid temporally-evolving model of Pasquetti (2011), and addresses most of the aforementioned problems. The hybrid method uses inflow conditions generated from a well-resolved body-inclusive simulation to perform a separate temporal simulation in the case of Pasquetti (2011) or spatially-evolving simulation in the work of VanDine et al. (2018) without including the body. By doing so, the amount of required points is substantially reduced since the flow near the body does not have to be resolved. This important reduction of the computational cost allows one to extend the domain farther downstream to gain insight in the far wake.
Here, we use a hybrid method that combines a body-inclusive simulation and a spatially evolving body-exclusive simulation to study the stratified high-Re far wake of a slender body for the first time. The Reynolds number is set to Re = ∞ / = 10 5 and two levels of stratification are used, Fr = ∞ / = 2 and 10. The simulation at Fr = 10 allows us to study the evolution of a weakly stratified wake in a domain that spans 80 . Additionally, Fr = 2 is chosen because it is close to the critical Froude number for a 6:1 prolate spheroid Fr = ( / )/ = 6/ . At the critical Froude number, the size of the separation region is strongly reduced by the lee waves (Ortiz-Tarin et al. 2019). These choices also allow us to compare our results with the findings of Chongsiripinyo & Sarkar (2020) (hereafter referred as CS20) regarding the stratified wake of a disk.
In CS20, the stratified wake of a disk at Re = 5 × 10 4 is studied at Fr = 2, 10, 50, ∞. Apart from a detailed analysis of the decay rates of the mean and turbulent quantities, CS20 links the general evolution of stratified homogeneous turbulence (Brethouwer et al. 2007;de Bruyn Kops & Riley 2019) with the evolution of the wake turbulence. As the disk wake evolves, the influence of buoyancy is 'felt' by the turbulent motions at progressively smaller scales. First the mean flow and the large scales and later the r.m.s. velocities are affected by stratification. Simultaneously, the horizontal eddies start gaining energy. Based on the strength of these effects, three distinct stages can be identified: weakly, intermediate and strongly stratified turbulence. In CS20, the transition between these regimes is examined and parameterized using local Froude and Reynolds numbers. Zhou & Diamessis (2019) also examined these transitions and their link with the evolution of stratified homogeneous turbulence using temporal simulations.
The present work is the continuation of Ortiz-Tarin et al. (2021) -referred to as ONS21 -where the unstratified wake of a 6:1 prolate spheroid with a turbulent boundary layer was studied and compared with a large number of simulations and experiments. In our previous study we found that the particularities of the slender body wake, e.g., small recirculation region, low entrainment, large defect velocity, bimodal distribution of the turbulent kinetic energy, among others affect the wake decay significantly. In this study we are analyzing how Focus on Fluids articles must not exceed this page length Stratified slender body wakes 5 these features affect the evolution of the stratified wake. We also analyze the simulations of CS20 to closely compare our results with the stratified bluff body wake.
Some of the questions we want to answer are: do the stratified decay laws and their transition points depend on the shape of wake generator? how does the turbulence evolve in stratified slender body wakes and are there difference with bluff body wakes? how does the phase-space evolution of the stratified turbulence compares between bluff and slender body wakes? In broader terms, we attempt to find whether a turbulent stratified wake retains some imprint of the wake generator in the mean and turbulence evolution.
A description of the solver and the methodology is given in §2. The wakes are visualized in §3. The decay of the mean wake properties is analyzed in §4. Finally, the evolution of the turbulence and the phase-space analysis of the wake are presented in sections §5 and §6, respectively. The study is concluded in §7.

Methodology
To study the far wake of a slender body at a high Reynolds number we use a hybrid simulation. The hybrid model combines two simulations: body-inclusive (BI) that solves the flow past the wake generator and body-exclusive (BE) that resolves the intermediate and far wakes. Here, we use a spatially-evolving simulation following the procedure validated by VanDine et al. (2018). In the implementation, data from a selected cross-plane in the BI simulation is interpolated on to a new grid and used as an inlet boundary condition for the BE stage. This procedure allows us to alleviate the natural stiffness of the wake problem. Whereas the BI simulation is designed to capture the turbulent boundary layer and the flow separation, the BE simulation resolves the turbulence in the wake. Both the grid size and the time step required to solve the turbulent boundary layer are much smaller than those needed in the intermediate and far wakes. This method leads to significant savings in computational cost without compromising accuracy.
The setup and the solver here are the ones used in ONS2021 with the addition of stratification. Both simulations solve the three-dimensional Navier-Stokes with the Boussinesq approximation in cylindrical coordinates. The solver uses a third-order Runge-Kutta method combined with second-order Crank-Nicolson to advance the equations in time. Secondorder-accurate central differences are used for the spatial derivatives in a staggered grid. A wall-adapting local eddy viscosity (WALE) is used to properly capture the turbulent boundary layer dynamics (Nicoud & Ducros 1999). Both the BI and BE simulations use Dirichlet boundary conditions at the inflow, convective outflow and Neumann at the radial boundary. Similarly to Ortiz-Tarin et al. (2019) a sponge layer is added to the boundaries to avoid the spurious reflection of gravity waves.
An immersed boundary method (Balaras 2004;Yang & Balaras 2006) is used to resolve the flow past a 6:1 prolate spheroid at zero angle of attack. A numerical bump is introduced on the surface of the body to accelerate the transition of the boundary layer to turbulence. The annular bump is located where the surface favorable pressure gradient is nearly zero. This location is found at approximately 0.5 from the nose. The radial extent of the bump is 0.002 (∼ 15 wall units) and the streamwise extent is 0.1 .
The stratification is set by a linear background density profile characterized by the Froude number, Fr = ∞ / , where is the buoyancy frequency. Three levels of stratification are simulated: Fr = 2, 10 and ∞. Fr = 2 is close to the critical Froude number Fr = 6/ for the 6:1 spheroid at which the suppression of turbulence in the wake by stratification is optimal (Ortiz-Tarin et al. 2019). Fr = 10 is a moderate level of stratification closer to oceanic values. Finally Fr = ∞ is the unstratified case which will be used as a reference (ONS21). The cylindrical coordinate system is ( , , ) with the origin at the body center. For  Table 1: Parameters of the body-inclusive simulation of prolate 6:1 spheroid. − and + are the upstream and downstream distances from the wake generator.
The BI grid is designed to resolve the turbulent boundary layer and the small-scale wake turbulence. The turbulent boundary layer is resolved with Δ + = 40, Δ + = 1, and Δ + = 32. There are 10 points in the viscous sublayer and 130 across the buffer and log layers. The mean velocities and turbulence intensities within the boundary layer were validated against existing studies (Posa & Balaras 2016;Kumar & Mahesh 2018) and the law of the wall. Additionally a grid refinement study was performed to guarantee the independence of the statistics to the grid choice.
In the wake, the peak ratio between the grid size and the Kolmogorov length = ( 3 / ) 1/4 , in both BI and BE domains is max(Δ / ) = 7.5, max(Δ / ) = 6, and max( Δ / ) = 5. A figure showing the ratio between the Kolmogorov scale and the grid resolution can be found in ONS21 (figure 2). Additionally, the unstratified wake decay coincides with all the previous existing numerical and experimental works on slender body wakes (see figure 1 of ONS21).
The domain size in the stratified cases is large so that internal gravity waves are weak before reaching the sponge region near the walls. The total number of grid points across BI and BE domains is approximately 1.5 billion in the unstratified case and 2 billion in the stratified simulations. Tables 1 and 2 include the most relevant parameters of BI and BE simulations, respectively. Further details on the grid design can be found in section 2 of ONS21. Once the flow has reached statistically steady state, the statistics are obtained by temporal averaging, denoted by · . Instantaneous quantities are written with lower case, mean quantities with upper case, and fluctuations with prime. In the stratified cases the average is performed over 270 / ∞ , approximately three flow-throughs. In the unstratified simulation flow statistics are obtained through the temporal (over 100 / ∞ ) as well as azimuthal averaging. Apart from temporal averaging, some statistics are obtained from cross-wake area integration denoted by {·}. Unless otherwise indicated, the integral is performed over a cross-section of radius 4 . All the flow statistics presented here die out well before they reach the limit of the integrated region.

Stratified slender body wakes
Reported velocities and lengths are normalized with the free-stream velocity ∞ and the body minor axis , respectively. The normalized streamwise distance from the center of the body is also measured as a function of the buoyancy frequency and the time. The time in the axis refers to time measured by an observer attached to the mean flow that sees the body move at a speed − ∞ . A Galilean transformation yields / = .
To compare the stratified wake of the 6:1 spheroid with that of a bluff body we use the body-inclusive disk wake simulations of CS20. The solver used in CS20 is the same as the one used here although, instead of using the WALE closure model, CS20 uses a variant of dynamic Smagorinsky. The eddy-viscosity model was changed in the spheroid simulations since WALE was demonstrated to capture the behavior of the turbulent boundary layer with the resolution used in the present wall-resolved LES. Both sets of simulations are very well resolved and have a small subgrid contribution -see ONS21 and CS20 -hence the validity of the comparison. Further details of the simulations can be found in ONS21 and CS20. The main parameters of disk simulations are shown in table 3. Figure 1 shows instantaneous snapshots of the near wake of a spheroid and a disk at Fr = 2 and Fr = 10. At both Fr, the near wake structure of the two bodies is very different. Compared to the spheroid with turbulent boundary layer (TBL), the disk wake has a large recirculation region (∼ 2 ), as shown by the red isolines in figure 1. This large recirculation region oscillates (Rigas et al. 2014) and generates a vortex shedding structure that is advected downstream (Nidhan et al. 2020). In a spheroid with TBL, the recirculation region is very small (∼ 0.1 ) and is surrounded by the small scale turbulence of the boundary layer. As a result, the near wake is highly organized and large scale oscillations are not observed in the near wake (Jiménez et al. 2010;Kumar & Mahesh 2018;Ortiz-Tarin et al. 2021). Only further downstream, does the wake begin to show a helical structure. This change in the structure of the slender body wake has been found to lead to a change in the decay rate and dissipation scaling in the unstratified wake (ONS21). In the following sections, we will analyze how the differences between the near wake of a disk and that of a spheroid lead to distinct trends of mean and turbulence evolution in a stratified environment. But first, let us describe different snapshots of the spheroid intermediate and far wakes. Snapshots of the disk intermediate and far wakes can be found in CS20. Figure 2 shows an instantaneous visualization of the spheroid Fr = 2 wake in the centervertical and center-horizontal planes. One of the distinctive features of the spheroid wake is that at Fr ∼ (1) the separation of the boundary layer can be strongly modulated by the steady lee waves (Ortiz-Tarin et al. 2019). This interaction between the lee waves and the wake is particularly strong when the Froude number is close to a critical Froude number Fr = / , where is the body aspect ratio. When Fr ≈ Fr , half the wavelength of the lee wave ( / = 2 Fr) coincides with the length of the body and the size of the separation region is reduced. The flow is then in what is called resonant or saturated lee wave regime, (Hanazaki 1988;Chomaz et al. 1992). At low Reynolds numbers, this effect can lead to the relaminarization of the turbulent wake (Pal et al. 2016;Ortiz-Tarin et al. 2019). In the present case, figure 2(a) reveals that, even at Re = 10 5 , the wake height is strongly modulated by the waves, although the wake is not relaminarized due to the high Re of the flow. For example, the wake height exhibits oscillations with a wavelength of / = 2 Fr = 4 . The modulation of the wake by the waves leads to an unusual configuration in the intermediate wake ( = 20 − 40) where the wake width is smaller than the wake height , figure 2(a) and (c). In these figures, sinuous oscillations are observed only in the horizontal plane (figure 2d) due to strong stratification. These horizontal sinuous instabilities contrast with the lee wave induced varicose modulation in the vertical plane. As the wake evolves, the < configuration transitions to the expected > . In this late region, the smallscale turbulence of the boundary layer has been dissipated and a layered-layer structure is observed in the vertical plane, figure 2(b). The qualitative trends of and discussed here are quantified in §4.
The main features of the Fr = 10 wake can be observed in the instantaneous snapshots of figure 3. The near wake, figure 3(a,c), is thin and carries the small scale turbulence generated in the boundary layer. Similar to the unstratified wake of ONS21, in the < 20 region, it has a quasi-cylindrical structure. Only after ≈ 20, a helical structure develops. In the unstratified wake, the oscillation found at ≈ 20 is present until the end of the domain. Here, the Fr = 10 wake does not show major oscillations after ≈ 30. Stratification restrains the vertical motions in the wake and enhances the horizontal spread as can be seen in the visualization of the late wake, figure 3(b,d). Unlike the Fr = 2 wake, the horizontal and vertical Fr = 10 wake extent grow monotonically with increasing downstream distance.

Evolution of the mean flow in spheroid and disk wakes
4.1. Evolution of the mean defect velocity ( ) The decay rate of the mean defect velocity = ∞ − shows the different stages in the evolution of a wake. In a stratified environment, wakes transverse the 3D, NEQ and Q2D regimes (Spedding 1997). Figure 4 compares the decay of among the unstratified, Fr = 10, and Fr = 2 spheroid and disk wakes. To facilitate a one-to-one comparison, we present the disk data in the domain 3 80, coinciding with the domain of the spheroid wake. Since is measured from the center of the body, the location of = 3 is in the near wake for the disk and is at the terminus of the body for the spheroid. The unstratified spheroid wake (figure 4a) shows a transition between the classical high-Re decay ∼ −2/3 to ∼ −6/5 at ≈ 20 coinciding with the development of a helical structure (ONS21). The Fr = ∞ disk wake decays as ∼ −0.9 from 10 65 and transitions to the classical high-Re decay of −2/3 afterward (CS20), as shown in figure 4(b). Compared to the disk, the Fr = 10 and Fr = ∞ spheroid wakes have a higher value of , owing to weaker near-wake entrainment and the slower development of slender body wakes. In the weakly stratified Fr = 10 regime, the defect velocity of the spheroid wake (figure 4a) evolves similarly to the unstratified wake until ≈ 3.5, when the decay rate slows down. However, at the same value of Fr = 10 but for the disk wake (figure 4b),

Ortiz, Nidhan and Sarkar
deviates from the unstratified case at ≈ 1. Based on , the end of the 3D region and the beginning of the NEQ region of the spheroid wake occurs at ≈ 30 whereas in the disk it occurs at ≈ 10. We discuss the reason behind this delayed deviation of the spheroid-wake from its unstratified counterpart in §5.
At Fr = 2, there are significant differences in evolution between the disk and spheroid wakes. In the Fr = 2 spheroid wake, shows an increased decay rate from the beginning. Although not shown here, the wake establishment is affected similarly to the Fr = 1 wake of the 4:1 spheroid of Ortiz-Tarin et al. (2019), where there was no 3D regime. The boundary layer evolution on the body and the separation are affected by stratification. At ≈ , there is a sudden change in the decay rate due to the lee-wave-induced oscillatory modulation (Pal et al. 2017) observed in the 3 10 region. This oscillatory modulation gets weaker downstream as the lee-wave amplitude decreases with the distance from the source.

At
≈ , the wake transitions to the NEQ stage where exhibits a slower decay compared to both the preceding stage and the following stage which commences at ≈ 15. Fitting a power law to the NEQ stage between = 6 − 25 results in a decay with −0.266 . This decay is close to the -1/4 decay in the NEQ regime found in the experiments of Spedding (1997) and later in numerical simulations (Diamessis et al. 2011;Brucker & Sarkar 2010;Redford et al. 2015;Pal et al. 2017). More details about the fitting strategy can be found in ONS21. At ≈ 15, the spheroid wake transitions to the Q2D regime with a sharper decay and a power-law fit between = 30 − 80 results in ∼ −0.72 , which is close to the −0.75 behavior established by Spedding (1997) for the Q2D regime. The Fr = 2 disk wake shows a very different behavior. Until at least = 125 ( = 62.5) -the full extent of the computational domain -the disk wake exhibits no transition to the Q2D regime. Instead, after transitioning to the NEQ regime with a power law of ∼ −0.18 , the disk wake stays in that regime.
Thus, the NEQ regime in the spheroid wake at Fr = 2 is significantly shortened compared to the disk wake, with it starting at ≈ and ending early at ≈ 15 when Q2D commences. In the experiments of Spedding (1997), the NEQ regime is reported to last until ≈ 40. Other temporal studies have found that the span of the NEQ regime depends on the Reynolds number. For example, in temporal simulations, Diamessis et al. (2011) found an increase of the NEQ duration to ≈ 50 when the Reynolds number increased to Re = 10 5 . Brucker & Sarkar (2010) found a transition to a Q2D-type power law at ≈ 100. Only Redford et al. (2015) observed an earlier transition around ≈ 25. The reasons behind the early arrival of the Q2D regime in the spheroid Fr = 2 wake will be discussed in §5.

Evolution of the mean horizontal ( ) and vertical ( ) lengthscales
The evolution of the mean wake dimensions in the spheroid and disk wakes is shown in figure 5(a,c) and (b,d), respectively. is defined such that ∞ − ( ) = 1 2 . The subscripts { , } indicate that these measures have been taken in the vertical and horizontal planes so that they represent the half height and the half width.
The wake of a slender body is generally thinner than its bluff body counterpart. Compared to the disk wake of CS20, the present unstratified wake is smaller by a factor of about 3 -contrast black lines in figure 5(a) to (b). The difference in wake size stems from the different near-wake features. Here, the initial non-dimensional wake width is around 0.2 . We find that disk ≈ 1.11 and spheroid ≈ 0.13 resulting in Besides the initial dimensions, the near-wake growth rates of the spheroid and disk are also very different. Whereas in the = 3 − 20 region the spheroid unstratified wake grows as ∼ 0.2 , the disk wake grows as ∼ 0.45 . Later, the growth rate of both wakes becomes comparable but the difference in size is already established and dictated by the near wake. The evolution of the Fr = 10 spheroid wake height ( ) is similar to its unstratified counterpart until ≈ 3.5 where the growth of slows down. While remains almost constant beyond ≈ 4, keeps increasing with a growth rate of ∼ 1/3 . In the Fr = 10 disk wake, the deviation from the Fr = ∞ case happens at ≈ 20 ( ≈ 2). Interestingly, after ≈ 2, the disk wake shows a continuous decrease in wake height. of both spheroid and disk wakes at Fr = 10 closely follow the trend of the corresponding unstratified wake. See figures 5(c,d).
The wake dimensions at Fr = 2 for both disk and spheroid show oscillations with a wavelength of / = 2 Fr. This reveals the influence of the steady lee waves especially on the wake height, see figure 5(a,b). From = 1 − 15 the oscillations of and are consistent with the conservation of momentum deficit, i.e., to counteract the contraction of caused by buoyancy, is enhanced. Note that these initial oscillations are of similar amplitude in both disk and spheroid. However, the relative change over the initial wake dimensions is much more pronounced in the spheroid wake (∼ 10 times) owing to its initial thinness. The influence of the lee waves on the spheroid Fr = 2 wake dimensions is illustrated by radial velocity contours in figure 6. The wake width contracts significantly in the region where the vertical velocity of the lee wave induces a rapid increase of wake height. Starting at = 20, the width of the spheroid wake shows a rapid growth ∼ 1.3 corresponding with: (i) the development of the horizontal wavy motions observed in figure 2(d) and (ii) the arrival of the Q2D stage with ∼ −3/4 . The growth of remains constant at a rate of 0.25 . Both, the very rapid growth of and the sustained increase in of the spheroid wake, are very different from the trends in the Fr = 2 disk wake. In the disk wake, we find that, instead of an increase, the vertical height exhibits a decrease ( ∼ −0.18 ) at 20. Furthermore, grows at 1/3 , a moderate rate relative to its rapid growth rate in the disk wake.

Comparison of flow topology between stratified spheroid and disk wakes
The difference between the spheroid and disk with regards to the evolution of mean length scales ( and ), particularly at Fr = 2, points toward qualitative differences in the flow topology. To further characterize these differences in the Fr = 2 wake, figure 7 shows contours of mean streamwise velocity ( ) at different streamwise locations for the spheroid (top two rows) and the disk (bottom row). In each panel of figure 7, the right half shows and the left half shows turbulent kinetic energy (TKE), = ( 2 + 2 + 2 )/2. For the sake of brevity, we have not included contours of the Fr = 10 disk and spheroid wakes since their topology is similar -the mean can be well approximated by a verticallysqueezed two-dimensional Gaussian while the TKE evolves from a bimodal (off-center peaks) distribution in the radial direction to a Gaussian at intermediate to late streamwise distances. In the case of the disk, the TKE evolves as a two-dimensional Gaussian right beyond the recirculation region.
We first discuss the disk wake (bottom row). The mean shows a monotonic spread in the horizontal direction and resembles the shape of an ellipse or a two-dimensional Gaussian squeezed in the vertical direction. This shape does not change until the end of the computational domain at ≈ 125. The TKE for the disk wake also has a similar vertically squeezed appearance.
Turning to the spheroid wake, we find that its turbulence topology is different from that of the mean. In the region 5 15, TKE shows two off-center peaks reminiscent of the TBL shedding from a slender body ( 2019), where the Fr = 4/ ≈ 1 wake of a 4:1 spheroid was studied. Note that in this stage, the wake is thinner than taller, i.e., < . At ≈ 20, TKE starts transitioning from a bimodal distribution to a single peak near the center-horizontal plane. In the next section, we will analyze how the horizontal contraction Contours of mean streamwise velocity in the right half and TKE in the left half of each contour. Contour limits are between the minimum (red) and maximum (white) values of the respective quantity at a given with ten levels in between. Radial extent span till = 1 and = 4 for the spheroid and disk contours, respectively. The disk wake is larger than the spheroid wake as can be inferred from the = 1 circle in the bottom row.
of the mean wake between = 20 and = 40 results in an increased horizontal mean shear, resulting in the maximum TKE being produced close to the center-horizontal plane. This leads to a transition in the TKE topology from bimodal distribution to a squeezed-Gaussian distribution at 40. By ≈ 60, has been organized into two distinct layers while TKE is sustained between these two vertically off-center layers. Note that the multi-layered mean flow structure at late in the Fr = 2 spheroid wake is reminiscent of the layered structure of the Q2D regime Spedding (1997).
Previously, temporal simulations (Gourlay et al. 2001;Brucker & Sarkar 2010;Redford et al. 2015) have shown that the mean and the turbulence can evolve differently. Indeed, the effect of buoyancy is 'felt' very differently by the large and the small scales in the flow. The general trend is that, in the late wake, the turbulence occupies a smaller and smaller vertical fraction of the mean defect as time passes (Redford et al. 2014). Instead, the finding here for the spheroid wake is the combined effect of having a wake in saturated-lee-wave state and initial off-center peaks of TKE, established by the TBL separation. These characteristics of the flow have not been not captured in temporal simulations since they have not accounted for the wake generator and the steady lee waves. x 10 −4

Evolution of the turbulent flow in spheroid and disk wakes
The energy of the flow is contrasted between spheroid and disk wake in this section. The turbulent kinetic energy (TKE also denoted as ), turbulent potential energy (TPE, ), mean kinetic energy (MKE, ) and mean potential energy (MPE, ) are defined as: where = 2 / 2 2 . In what follows, trends of area-integrated values, denoted by {·}, of these energy measures are reported. Area-integrated quantities are preferred because the peaks of mean and turbulence in stratified slender wakes are often times off-centered as can be seen in figure 7. The integration allows for a uniform comparison across cases. Also note that temporally averaged quantities are denoted by the angled brackets · .

Evolution of TKE, spectra and PE-to-KE ratios
The evolution of the area-integrated TKE of the disk and spheroid wakes is shown in figure  8. The most noticeable aspect is that, for all Fr numbers, { } in spheroid wakes is an order of magnitude smaller than in corresponding disk wakes. Although it is not shown here, we found a similar result for the area-averaged TKE (over an area of = 4 cross-section), instead of area-integrated TKE.
The decay of the unstratified wake is studied in more detail in ONS21 for the spheroid and CS20 for the disk. In unstratified flow, { } decays following { } ∼ −2/5 for the spheroid and the disk wake follows { } ∼ −2/3 . Note that these fits are empirical. Both decay rates are consistent with the decay of the peak TKE ( ) and the growth of wake width ( ) under the self-similarity framework, i.e., { } ∼ 2 . Specifically, in the case of the spheroid ∼ 3/5 and ∼ −8/5 (Ortiz-Tarin et al. 2021) and, in the case of the disk, ∼ 1/3 and ∼ −4/3 (Chongsiripinyo & Sarkar 2020). At Fr = 10 and for both the spheroid and the disk, { } deviates from the unstratified case at between ≈ 1 and ≈ 3, see figure 8(a,b).Interestingly, while the TKE is affected by stratification in the similar range for both wake generators, the mean flow showed a different behavior. In the disk Fr = 10 wake, deviated from its unstratified counterpart at ≈ 1 while, in the spheroid, this change occurred later at ≈ 3. The onset of deviations from the unstratified case will be explained in more detail in the following subsections.
At Fr = 2, there is a striking influence of the wake generator on the evolution of TKE. Whereas in the disk wake, the TKE decays monotonically, the far wake of the spheroid displays a rapid increase.
The disk wake shows a monotonic decay in { } and its individual components throughout 3 < < 80. Compared to the horizontal components, { } shows a sharper decay after ≈ 10 and turbulence anisotropy progressively increases. In the spheroid Fr = 2 wake, { } decays rapidly until ≈ 10. However, after ≈ 10, the decay slows down and is followed by a period of sustained growth starting at ≈ 20 and lasting until the end of the domain. The region of TKE growth coincides with the development of the large scale horizontal motions observed in figure 2(d) and the rapid growth of shown in figure 5(c). It also coincides with the accelerated decay rate of starting at ≈ 20. Notice that right before the start of the rapid increase of TKE, the Fr = 2 wake has a configuration where < (figure 7d). The horizontal response of the flow to the strong lee waves is what sustains this configuration. It is only after their strength subsides that the control on the wake is released to allow the horizontal wavy motion to develop. The rapid development of this motion coincides with the rapid increase in horizontal TKE, namely { } and { } as seen in figure 8(c,e). To the best of the authors' knowledge this is the first wake study, resolving the flow at the body, that observes an increase of fluctuation energy with downstream distance instead of its usual decrease. To further quantify the horizontal wavy motions observed in figure 2(d), the energy spectra of the spanwise velocity are computed at the centerline. These spectra are compared between locations before (figure 9a) and after (figure 9b) the start of the the TKE increase associated with the Q2D regime. The spectra before < 40 do not show preferential energization of the low frequencies. This finding is consistent with the visualizations in figure 2  In the Fr = 2 disk wake, the vortex shedding mode at St = / ∞ ≈ 0.13 − 0.14 is dominant throughout the whole domain. The horizontal meanders, which are prevalent in the spheroid wake, are absent in the disk wake at least until the end of the domain at / = 125. Since the vortex shedding mode, its long-lasting effect on the wake, and its internal wave field are described in detail by Nidhan et al. (2022), we do not discuss these aspects further.
The spheroid wake has significantly lower TKE content relative to the disk wake and also a different distribution of the mean momentum. The content of potential energy relative to that of kinetic energy is also of interest. Figure 10 shows the ratio of area-integrated potential energy (PE) to kinetic energy (KE). Both fluctuating (10a) and mean (figure 10b) components are shown at Fr = 2 and 10 for both wake generators. At Fr = 10, the turbulent PE-to-KE ratio increases steadily in both disk and spheroid cases indicating an increasing influence of buoyancy on turbulence (figure 10a). In the Fr = 2 wakes, the turbulent PE-to-KE ratio peak around ≈ 30 in both cases and decay afterward. By the end of the measurement region at = 80, the turbulent PE-to-KE ratios are similar across Fr = 2 and 10 wakes. Thus, stratification and body shape do not qualitative affect the ability of turbulence to stir the density field in the intermediate and far wake. Quantitatively, the TPE-to-TKE ratios are somewhat higher for the disk relative to the spheroid for both wakes.
The mean PE-to-KE ratios in the Fr = 10 wakes (figure 10b) are minuscule compared to their turbulent counterparts. At Fr = 2, the mean-based ratio is much larger, particularly close to the wake generators, pointing towards a strong influence of the steady lee waves. Both disk and spheroid mean-based ratios oscillate with a characteristic lengthscale corresponding to the wavelength of steady lee waves at Fr = 2. It is particularly revealing as to how much larger is the magnitude of { }/{ } in the spheroid wake compared to the disk, as it explains why the spheroid flow is so strongly modulated by the lee waves. It is worth noting that comparison of the absolute value of MPE between the disk and the spheroid reveals that it is the disk that has the larger MPE, an order of magnitude larger. The amplitude of the lee wave generated by the disk is larger than that of the spheroid by about a factor of 2.

Analyses of the spheroid TKE budget terms
To quantitatively understand the origin of TKE increase in the Fr = 2 spheroid wake, we look into the different terms of the TKE transport equation: where is the turbulent production, is the turbulent dissipation and denotes the turbulent buoyancy flux. These quantities are defined by: is the strain-rate tensor and = −2 is the subgrid stress tensor. The contribution of the subgrid term to the TKE transport equation is found to be small.
The turbulent production is a source of TKE and a sink in the MKE equation. The turbulent buoyancy flux transfers energy between TKE and TPE and the turbulent dissipation is a sink of TKE. Along with the sinks and sources of energy, there is a term responsible for the spatial redistribution of TKE, the turbulent transport, The turbulent transport redistributes energy, primarily in the -plane, and its contribution to the area-integrated budget is negligible. The production term in equation (5.3) can be further expanded as follows: For the stratified wakes at hand, , << . Hence equation (5.6) can be further simplified as below: In a wake developing in the direction, is primarily dominated by the transverse shear terms since | / |, | / | >> | / |. Hence in the discussions on { } to follow, we focus on the transverse terms, namely = and = . Figure 11 shows the evolution of the area-integrated terms in the TKE transport equation of the spheroid wake. We do not present the TKE transport terms of the disk wakes here as they are presented and discussed in great detail in CS20. We first discuss the trends of the area-integrated production, { }, and its components in Fr = 2 and 10 spheroid wakes. This is followed by a discussion of the area-integrated dissipation, { }. Thereafter, the variation of { }/{ } sheds light on why an increase in TKE is observed in the Fr = 2 spheroid wake for > 30. { } in the Fr = 10 spheroid wake decays at a rate comparable to its unstratified counterpart until ≈ 3 when it starts to deviate (figure 11a). in the Fr = 10 spheroid wake also starts deviating around ≈ 3 (figure 4a). At ≈ 3, the contribution of = starts declining due to a reduction in (not shown here for brevity). At ≈ 4.5, the contribution of the lateral production exceeds that of the vertical production . The reduction of { } in the Fr = 10 wake coincides with the transition of to a slower decay rate. In figure 11(b), we note that { } starts decaying faster than its unstratified counterpart between 3 < < 4 (or 30 < < 40). In figure 4(a), we see that this the location where of the Fr = 10 spheroid wake starts deviating from the Fr = ∞ spheroid wake. At the same location that { } starts decaying rapidly, there is a maximum in the buoyancy flux { } (figure 11e) and a slowdown in the growth of (figure 5a). Buoyancy is starting to affect the wake decay.
The initial value of { } in the Fr = 2 wake is higher than in the unstratified case. The distinct separation pattern and the vertical contraction of the mean flow (figure 5a) lead to an increase in the vertical shear and hence { }, see red dashed lines in figure 11(b). This initially high value of { } rapidly decays as increases, leading to a reduction in the vertical shear. Note that the beginning of the -based NEQ regime at ≈ also coincides with this reduction of { }. However, the mechanism at Fr = 2 is different from that of the Fr = 10 wake. Whereas at Fr = 10, the reduction of causes the decay of { }, at Fr = 2 the decay of is accompanied by a sudden reduction of (not shown here). The reduction in vertical shear is driven by the expansion of due to lee wave modulation in figure 5(a).
As the modulation of the wake by the lee waves continues, the horizontal production is enhanced by the strong reduction in (figure 5c) in the NEQ region during 7 < < 20. The value of { } overtakes { } at ≈ 5. The magnitude of { } remains nearly constant until the end of the domain. Figure 11(f) shows the maximum of the horizontal mean shear ( / ), i.e., max ( , =0) / , in the central streamwise-horizontal plane for both wakes. The mean shear in the Fr = 2 wake increases in the region 10 < < 30, exactly in the region where the TKE decay starts plateauing (see figure 8a). The enhanced mean horizontal shear at Fr = 2 prevents the horizontal production from decaying monotonically unlike Fr = 10. Figure 11(c) shows the evolution of { } as a function of for Fr = ∞, 10 and 2 spheroid wakes. In the Fr = 10 case, the decay of is similar to that of the unstratified wake until ≈ 3 after which the decay rate increases slightly. The Fr = 2 wake dissipation shows a sharper decay until ≈ 20. After ≈ 20, the decay rate appears to be closer to the Fr = ∞ decay rate of −7/5 (Ortiz-Tarin et al. 2021).
Since the turbulent dissipation keeps decreasing monotonically in the Fr = 2 spheroid wake (figure 11c), while the production tends to asymptote for > 30 (figure 11a), the value of { }/{ } becomes greater than 1 in the Fr = 2 spheroid wake, see figure 11(d), explaining the rapid increase of TKE beyond ≈ 30. Note that the { }/{ } ratio oscillates with the characteristic wavelength of the lee waves revealing the strong influence of wave-related buoyancy effects on the energetics of the Fr = 2 wake.
One of the features of the arrival of the NEQ regime reported in previous studies is the radiation of internal gravity waves (Abdilghanie & Diamessis 2013;De Stadler & Sarkar 2012;Rowe et al. 2020). In these spheroid wakes, we find that the integrated wave flux remains negligible compared to the other terms in the TKE budget and hence is not shown here. The small magnitude of the turbulent wave flux is consistent with the findings of Meunier et al. (2018) who found that the magnitude of the wake-generated waves depends on the body drag coefficient -and the 6:1 spheroid has a very low drag compared to bluff bodies.

Early arrival of the Q2D regime in the Fr = 2 spheroid wake
Remember that the Q2D regime in stratified wakes is characterized by a faster decay of the mean defect velocity ( ∼ −3/4 ) and the organization of wake into distinct layers in the vertical direction. Due to the strong effects of buoyancy, the fluid motions are primarily confined to the horizontal plane. However, the vertical variations in the velocity profiles still exist during the Q2D regime, primarily in the form of layered structures.
A key difference between the Fr = 2 disk and spheroid wakes is the early arrival of the Q2D regime in the spheroid wake. Whereas in the disk wake, CS20 did not observe a transition to the Q2D in a domain that extended up to = 125 ( = 62.5), here we observe a transition at ≈ 40 ( ≈ 20). The early transition in the spheroid wake is a consequence of the strong modulation of the intermediate wake by buoyancy, an effect that occurs for bodies with large aspect ratio ( / ) and, specifically, when Fr is near its critical value Fr = / . Figure 5(c) shows that in the Fr = 2 spheroid wake contracts in the region 5 < < 30 as a response to the expansion of (figure 5a) by steady lee waves. This phenomenon leads to the 'butterfly' shaped structure of mean ( figure 7(d,e)) which was also observed at a lower Reynolds number, Re = 10 4 in Ortiz-Tarin et al. (2019). In that study, which was performed at Re = 10 4 , the spheroid wake at critical Fr relaminarized. Here, at a higher Re = 10 5 , the flow response at the resonant state is quite different. The constriction of leads to an enhancement in the mean horizontal shear shown in figure 11(f). This enhancement significantly slows down the decay of horizontal production (figure 11b). While continues to decay, { }/{ } becomes > 1 leading to an increase of TKE for > 40 (figure 8a). In figure 7(d-f), we also see the maximum TKE location moving to the wake axis from its off-center location in the near wake. It is worth noting that the enhanced { } that acts as a source of TKE is a sink for the mean energy. The sharp increase in TKE leads to a faster decay of , and the Fr = 2 wake transitions to the Q2D regime early on, at ≈ 40.
In the disk wake, is initially 3 − 4 times larger than in the spheroid wake. While the lee wave modulation of is present in the disk wake as well (figure 5d), its amplitude relative the original value of is quite small. Hence the horizontal mean shear (not shown here) and the horizontal production in the Fr = 2 disk wake continue to decay unlike in the spheroid wake.
The arrival of the Q2D regime in the spheroid wake is also accompanied by distinctive features of the Q2D regime reported in previous literature. Figure 2  Ortiz, Nidhan and Sarkar 5.4. Late transition to the NEQ regime in the Fr = 10 spheroid wake At Fr = 10, the main difference between the disk and spheroid wakes is the location at which deviates from the unstratified counterpart, i.e., the transition point to the NEQ regime. In the spheroid wake this transition occurs around ≈ 3, whereas in the disk it occurs at ≈ 1, compare the Fr = 10 curves between figure 4 (a) and (b). The beginning of the NEQ regime is marked by a slowdown in the decay rate of in the stratified wake as compared to its unstratified counterpart. The NEQ regime is also characterized by the progressively increasing anisotropy between the horizontal and vertical components of TKE.
To understand better how this transition occurs we look into the mean kinetic energy (MKE) transport equation. The TKE transport equation was given in equation (5.3). In a similar fashion, the MKE transport equation is given by: where superscript denotes the mean counterparts of terms in equation (5.3). Broadly, the mean and turbulence quantities in stratified wakes deviate from their unstratified counterparts as a result of buoyancy. In the MKE transport equation, buoyancy can affect (a) the turbulence production as an indirect effect and (b) the MKE ↔ MPE transfer, through the mean buoyancy flux . In both the disk and spheroid Fr = 10 wakes, we find that the contribution of the mean buoyancy flux to the MKE transport is significantly smaller than the contribution of the production, particularly for 40. Now turning to the area-integrated production in figure 12(a), we find that the production at Fr = 10 displays a strong reduction from the Fr = ∞ wake later in the spheroid at ≈ 30 ( ≈ 3) compared to the disk where a similar strong reduction occurs at ≈ 10 ( ≈ 1). This explains the late transition of the spheroid wake to the NEQ regime, i.e. deviates from unstratified behavior later ( ≈ 3) for the spheroid than the ≈ 1 transition point for the disk. The decreased production in the Fr = 10 wakes of both bodies is a consequence of the decreased correlation (not shown here for brevity) in stratified turbulent shear flows (Jacobitz et al. 1997). In stratified wakes, this buoyancy effect has been observed experimentally by Spedding (2002) and numerically by Brucker & Sarkar (2010). The deviation between production of the Fr = 10 and ∞ wakes (figure 11a) is indeed caused by a reduction in its { } component, see blue dashed lines in figure 11(b).

Evolution of the local flow state and its trajectory in the phase space
In previous sections, we showed how the spheroid and the disk wake do not transition between the 3D, NEQ and Q2D regimes at the same . In this section, we examine the evolution of key local non-dimensional numbers describing the mean and fluctuating state to explore their roles. These non-dimensional numbers are local, streamwise-varying measures of stratification (Froude number) and the dynamical range of inertially dominated scales (Reynolds number). We also plot the trajectory of each wake in the Froude-Reynolds phase space. The mean vertical Froude number (Fr ) and the turbulent vertical and horizontal Froude numbers (Fr , Fr ℎ ) are defined as with turbulence length scales in the energy-containing range. Since the horizontal turbulent integral length scale ( ℎ ) is not easy to compute in a spatially evolving flow, we use the TKEbased horizontal half wake width ( ) as a surrogate following CS20. at a downstream location is calculated by ( = , = 0, ) = 0.5 ( = 0, = 0, ). The mean vertical Froude number Fr (figure 13a) becomes (1) in both disk and spheroid wakes at the location at which the decay of slows down signaling the beginning of the NEQ regime, marked by ≈ 1 for the disk and ≈ 3 for the spheroid. For the Fr = 2 cases, Fr in the spheroid wake starts dropping faster beyond ≈ 40 − the streamwise location where the wake transitions from NEQ to Q2D regime.
The turbulent Froude numbers play an analogous role to those based on the mean. When their values become (1), turbulence starts being affected by buoyancy. Fr is defined using which has a significance to shear instability. Defined with , Fr is proportional to −1/2 , where is the Richardson number of the fluctuating shear ( Riley & DeBruynKops (2003), CS20). When Fr becomes (1) (figure 13b) is also when the spanwise and vertical components of the TKE start showing anisotropy of the turbulence stress tensor and > as shown in figure 8 for the Fr = 2 wakes. See also figure 8 of CS20. The avid reader will notice that the four Fr = ℎ / curves could collapse if plotted against instead of . Indeed, these curves collapse in the (Fr , ) plot for 1 < < 10 and the factor of 5 present between the curves in figure 13(b) is simply the ratio of buoyancy frequency between Fr = 2 and Fr = 10. This result points to the fact that the spheroid and disk wakes are at a sufficiently high Re so that there is a range of small-scale wake turbulence with dynamics relatively unaffected by the large scales and, therefore, by neither the shape of the wake generator nor buoyancy. ℎ / is a quantity that reflects such dynamics and thus, on dimensional grounds, evolves as ℎ / ∼ −1 so that Fr = ℎ / ∼ ( ) −1 . We note that buoyancy does affect large scale wake features, which are affected by the body particularities -hence the lack of collapse in Fr . The values of the buoyancy Reynolds number Re = / 2 and the Reynolds of the Taylor micro-scale Re = √ / where 2 = 15 2 / (not included here for brevity) are very similar in magnitude for both disk and spheroid wakes at the two levels of stratification. This finding is consistent with the initial collapse of Fr when plotted against .
In figure 13(c), the Fr = 10 wakes of both disk and spheroid reach Fr ℎ ∼ (1) at ≈ 10 − 20, the location at which the { } starts deviating from its unstratified counterpart ( figure 8(a,b)). In the Fr = 2 spheroid wake, Fr ℎ < 1 throughout the domain and { } deviates from the unstratified decay from the very beginning in both disk and spheroid wakes. Overall, we find that Fr ℎ ∼ (1) is a good indicator of the deviation of { } from the unstratified counterpart. In contrast, Fr ∼ (1) marks the location at which turbulence anisotropy between the vertical and the spanwise components starts growing. Figure 13(d) shows the evolution of turbulence Reynolds number, Re ℎ = ℎ / with ℎ being the intensity of horizontal turbulent fluctuations. The value of Re ℎ (figure 13d) changes slowly as the wake evolves, remaining at (10 4 ) for the disk wake and at (10 3 ) for the spheroid wake.
A consolidated view of the evolution of the state of fluctuations is provided in phase space (Brethouwer et al. 2007;Zhou & Diamessis 2019;de Bruyn Kops & Riley 2019;Chongsiripinyo & Sarkar 2020). In the phase-space portrait, one axis measures the buoyancy effect on the large scales through the turbulent Froude number (Fr ℎ ) and the other axis is a measure of Reynolds number that is not Re ℎ but one which accounts for buoyancy in addition to inertia. Ozmidov-scale eddies are the largest eddies unrestrained by buoyancy. The Ozmidov scales = ( / 3 ) 1/2 and = ( / ) 1/2 lead to the definition of Re = / = / 2 as the so-called buoyancy Reynolds number. Another convenient measure of Reynolds number, which does not require explicit computation of the turbulent dissipation rate, is the buoyancy-weighted Reynolds number, Re ℎ Fr 2 ℎ (Billant & Chomaz 2001;Riley & DeBruynKops 2003). The Reynolds numbers based on buoyancy tend to decrease with downstream distance as buoyancy progressively increases in importance and limits the range of scales which are susceptible to 3D turbulent motions. As long as Re ℎ Fr 2 ℎ > (1), viscous effects do not dominate.
Following CS20, figure 14 shows the evolution all four wakes in the Re ℎ Fr 2 ℎ − Fr ℎ phase space, where Re ℎ = ℎ / , ℎ being the intensity of horizontal turbulent fluctuations. The flow evolves in the direction of the arrow from a state where buoyancy effects are weak (almost negligible) to a region characterized by the presence of stratified turbulence. Within the state of stratified turbulence, three different regimes can be demarcated: weakly, intermediately, and strongly stratified turbulence (WST, IST, and SST). Unlike the case of freely decaying turbulence, the mean velocity is also important here. Therefore, CS20 elected to distinguish between WST (where buoyancy begins to affect the mean velocity) and IST (where the turbulence anisotropy begins to be affected) as we do here. The SST regime is one where buoyancy effects on the large scales is very strong (Fr ℎ 1) but nevertheless the value of Re ℎ Fr 2 ℎ is sufficiently large so that viscous effects are not dominant at the large scales.
The slope in the Re ℎ Fr 2 ℎ − Fr ℎ plane is similar for both disk and spheroid wakes in the weak buoyancy and WST stages. The slope almost follows a line of constant Re ℎ , marked by a dashed line, signaling the slow evolution of Re ℎ compared to that of Fr ℎ . However, there are significant differences in the way the spheroid wakes transverse the phase space. The spheroid wake starts out thinner than the disk wake, leading to larger local Froude numbers. Also, despite the body-based Re being larger in the spheroid wake by a factor of two, the turbulence intensity in the disk wake is higher, hence the Re ℎ difference observed in figure 13(d). This difference in Re ℎ can be seen in the horizontal offset between the disk and spheroid lines in the phase-space representation.
Despite the larger spheroid body-based Reynolds number -Re = 10 5 instead of Re = 5 × 10 4 , turbulence in the spheroid wake is unable to access either the the IST regime or the SST regime while the disk wake is able to access these regimes of stratified turbulence. Furthermore, the increase in TKE, which is not observed in the Fr = 2 disk wake, reverses the trajectory of the Fr = 2 spheroid wake. To the authors' best knowledge, a reversing trend of phase-space trajectory has not been seen before in the stratified turbulence literature. These differences reveal that the phase space evolution, at least for Fr 10, depends on the features of the wake generator, e.g., its aspect ratio or the type of BL separation.
Based on the phase-space portrait alone, one may hastily conclude that stratified slender body wakes always experience a weaker buoyancy effect relative to bluff bodies. This is true for the Fr = 10 case, at least for the limited simulation time, in terms of the relative amount of TPE (figure 10) and the vertical scale ( ) deviation from the unstratified case -note that both are smaller in the spheroid wake than in the disk. However, we have also shown that for the Fr = 2 spheroid wake, compared to the disk wake, there is a much earlier onset of the Q2D regime. The mean and turbulence quantities are highly intertwined and the strong modulation of the mean flow by steady lee waves of the high aspect-ratio spheroid in the Fr = 2 case ultimately leads to an early entry into the Q2D regime. Hence, it is not entirely accurate to say that the slender body wakes are weakly affected by stratification. The present work calls for a need to generalize the parameter space of {Re, Fr} in turbulent shear flows to account for the mean flow field (also possibly its instabilities) to build a more comprehensive understanding of buoyancy effects in shear flows. In the case of wakes, the shape of the body generator is brought into play through the mean flow as shown here.

Discussion and final remarks
The high-Reynolds number stratified wake of a slender body has been studied using a highresolution hybrid simulation. The wake generator is a 6:1 prolate spheroid with a tripped turbulent boundary layer, the diameter-based Reynolds number is Re = 10 5 and the Froude numbers, namely Fr = ∞ / = {2, 10, ∞}, take moderate to large values. By comparing the spheroid wake with the disk wake of Chongsiripinyo & Sarkar (2020) (referred to as CS20), we are able to study the influence of the wake generator -slender versus bluff -in the establishment and evolution of stratified wakes.
The near wake of the 6:1 prolate spheroid with a turbulent boundary layer is characterized by a small recirculation region (∼ 0.1 ). The recirculation region is surrounded by smallscale turbulence that emerges from the boundary layer and the flow does not show strong vortex shedding at the body (Jiménez et al. 2010;Posa & Balaras 2016;Kumar & Mahesh 2018;Ortiz-Tarin et al. 2021). As a result, the wake is much thinner and develops slower than the wake of a bluff body like the disk, which has a large recirculation region (∼ 2 ) and vortex shedding from the body (Nidhan et al. 2020). These body-dependent features of the near wake were recently shown to affect the decay of the far wake in environments without density stratification (Ortiz-Tarin et al. 2021). In the present stratified simulations also we find substantial differences in the decay of the disk and spheroid wake. Particularly, we find that the starting locations of the non-equilibrium (NEQ) and the following quasi-2D (Q2D) regions of wake deficit velocity depend on the wake generator.
At Fr = 2 ≈ ( / )/ , the wake of a 6:1 prolate spheroid is in a resonant state. The half wavelength of the lee waves is equal to the body length and, as a result, the flow separation and the wake are strongly modulated by the waves. Whereas previous works had described this regime in laminar-separation configurations of a sphere (Hanazaki 1988;Chomaz et al. 1992) and a 4:1 spheroid (Ortiz-Tarin et al. 2019), the present results show that the influence of the waves persists even at high Reynolds numbers and with the separation of a turbulent boundary layer. At Fr = 2, the flow and the turbulence in the spheroid wake evolve very differently from the disk wake. Both the lack of strong shedding in the near wake (Ortiz-Tarin et al. 2021) and the strong modulation of the mean flow by the lee waves, lead to a wake with vertical and horizontal profiles of mean velocity that depart strongly from Gaussian. These features are not observed in the disk wake at Fr = 2, which shows a vertically-squeezed Gaussian topology and a weak imprint of lee waves on the wake dimensions.
At Fr = 2, both disk and spheroid wakes transition to the NEQ regime at ≈ . However the transition to the Q2D regime -with enhanced wake decay relative to the NEQ regime -is very different; whereas the spheroid wake transitions at ≈ 15, the disk wake does not access the Q2D regime in a domain that spans ≈ 60. Other bluff bodies, e.g., the towed sphere (Spedding 1997) show transition to the Q2D regime at ≈ 50, a location which is also delayed with respect to the spheroid wake. The early transition to the Q2D regime of the spheroid wake is driven by its strong modulation -horizontal contraction and expansion of the wake width -in response to the vertical contraction and expansion by the lee waves. This modulation has a particularly strong effect on the slender wake of a spheroid where the horizontal contraction is a large fraction of the wake width. The early start of the Q2D regime in the spheroid wake is accompanied by a sustained increase of turbulent kinetic energy (TKE), driven by an increase of the horizontal mean shear which acts on the turbulence of the separated boundary layer. The TKE increase is limited to the horizontal velocity with the spanwise component being strongest, having almost an order of magnitude larger energy than the vertical. Although coherent vortical structures and spanwise flapping are seen in the horizontal motion, pancake eddies are incipient and not fully formed at the end of the domain, / = 80.
At Fr = 10 also, there are differences between the disk and the spheroid wakes. Particularly in the spheroid wake, the beginning of the NEQ stage occurs later, at ≈ 3 instead of ≈ 1 ( ≈ 30 instead of ≈ 10). The difference in the start of the NEQ can be attributed to the value of the local mean Froude number Fr = /2 . As noted previously, the spheroid wake is thinner than the disk wake, the mixing in the near wake is weaker, and as a result the defect velocity in the intermediate wake is larger. These features increase the value of the spheroid wake local Froude number and delay the onset of the buoyancy effect that gives rise to the NEQ regime. Additionally, the analysis of the mean kinetic energy (MKE) transport terms shows that the onset of buoyancy effect on the mean flow of both the disk and spheroid Fr = 10 wakes is associated with the decreased energy transfer from MKE to TKE.
Taking the unstratified case as a base line, the effect of buoyancy in the spheroid Fr = 10 wake is observed earlier (at ≈ 1) on the decay of the TKE than its effect (at ≈ 3) on the decay of . In the spheroid wake at Fr = 10, the transfer from TKE to TPE is responsible for the enhanced decay of TKE at ≈ 1. On the other hand, the decrease in turbulent production at a farther downstream distance (compared to the disk Fr = 10 wake) in the spheroid Fr = 10 wake is responsible for the slowed decay of the mean defect velocity at ≈ 3. The decrease in the production is caused by a reduction in the correlation (Jacobitz et al. 1997;Spedding 2002;Brucker & Sarkar 2010).
Meunier & Spedding (2004) compared the evolution far into the stratified wake, up to ≈ 8000, among several body shapes that also included a 6:1 prolate spheroid and a circular disk. The body Reynolds number was Re = 5000 and their diameter based Froude numbers were Fr = 4 and 16. When normalized using eff = √︁ /2 instead of , the evolution of the peak defect velocity of different wake generators exhibited approximate collapse for 50 (see their figure 5b) with a Q2D decay rate of ∼ −0.75 . Similarly, the wake width in the horizontal of different shapes approximately collapsed for / eff > 400 to exhibit a growth rate of ∼ 0.35 . Since their 6:1 spheroid data starts from / ≈ 100 (see their figure 5a), a direct comparison is not possible with our simulations that end at at / = 80.
In regard to the mean defect decay, a major difference between Meunier & Spedding (2004) and our results is that the transition to Q2D power-law behavior appears earlier, around ≈ 15, in the Fr = 2 spheroid wake relative to the the Fr = 2 disk wake which does not transition to the Q2D decay rate until the end of the domain at = 62.5. The value of Re = 10 5 in the spheroid wake is larger here and it is possible that the features that we have linked to the early onset of the Q2D stage for the spheroid Fr = (1) wake, i.e., the instability that leads to horizontal meanders and also the enhanced TKE production, are inhibited by viscous damping at the lower Re of the experiments.
The present simulations, both of the disk and the spheroid, do not extend into the very far wake regime reached by their experiments. Future hybrid simulations or experimental work at higher Re that probe the very far wake would clearly be useful. At any given Re and Fr, it will also be of interest to look at how tripping affects the wake evolution for the different body shapes.
The simulations show that the buoyancy timescale alone is not sufficient to determine the state of the wake decay for both generators. However, we find that the value of the local turbulent and mean Froude numbers can be a good proxy to describe some aspects of the wake state. For both disk and spheroid wakes, Fr = /2 becomes (1) at the location at which the decay of slows down; Fr ℎ = ℎ / ∼ (1) marks the location at which the area-integrated TKE of the stratified wake starts deviating from the unstratified case; and Fr = ℎ / ∼ (1) signals the location at which anisotropy between the different TKE components starts growing.
The buoyancy-weighted Reynolds number (Re ℎ Fr 2 ℎ ) has been used widely in stratified Ortiz, Nidhan and Sarkar flow as a convenient surrogate for the buoyancy Reynolds number (Re = / 2 ) since it displays similar trends during the flow evolution and the two quantities can be shown to be proportional using classical inviscid scaling of the turbulent dissipation rate. The surrogacy is true for the stratified wakes considered here except for the spheroid Fr = 2 wake after its entry into the stage of Q2D wake decay. The horizontal fluctuation energy, therefore Fr ℎ , increases owing to horizontal meanders and flapping of the flow. However, continues to decrease, albeit at a reduced rate relative to the NEQ regime. The value of Re = (10) is not high in the Q2D regime realized here at Fr = 2. It remains to be seen if, in the Q2D regime at even higher body-based Reynolds number, the equivalence between Re ℎ Fr 2 ℎ and Re is recovered and whether the unusual upward trajectory seen here in {Fr ℎ , Re ℎ Fr 2 ℎ } phase space is also seen in {Fr ℎ , Re } space. The duration of the upward trajectory in phase space until the eventual downward shift toward the viscous regime is also of interest.
The differences between bluff body (disk) and slender body (6:1 spheroid) wakes illustrate the difficulty of finding a universal scaling for the high-Re stratified wake. The initial magnitude of for different wake generators and levels of stratification can be roughly scaled with the global Fr and the body drag coefficient (Meunier & Spedding 2004). However, the start and the duration of the NEQ regime cannot be assumed to be independent of the wake generator. We find that rather than a particular , the local mean Froude number is a good proxy for the onset of the NEQ regime in the mean defect velocity and the values of local turbulent Froude number provide guidance for the behavior of TKE, e.g., the onset of buoyancy effect as well as the location at which the ratio of vertical to horizontal TKE starts decreasing. We are unable to connect Froude number to the Q2D regime transition of the wake. More numerical and experimental work spanning different wake generators, different sources of turbulence including freestream turbulence, and longer downstream distances will be instrumental in building a comprehensive picture of the effect of initial/boundary conditions on subsequent wake evolution.