Anisotropic RANS Turbulence Modeling for Wakes in an Active Ocean Environment

: The problem of simulating wakes in a stratiﬁed oceanic environment with active background turbulence is considered. Anisotropic RANS turbulence models are tested against laboratory and eddy-resolving models of the problem. An important aspect of our work is to acknowledge that the environment is not quiescent; therefore, additional sources are included in the models to provide a non-zero background turbulence. The RANS models are found to reproduce some key features from the eddy-resolving and laboratory descriptions of the problem. Tests using the freestream sources show the intuitive result that background turbulence causes more rapid wake growth and decay.


Introduction
Oceanographic flows include a broad variety of turbulence-generating phenomena, and the associated unsteady motions are in general inhomogeneous, non-stationary, and anisotropic. The thermohaline stratification of the ocean introduces a conservative body force which must be considered when examining flows in such an environment. Numerous effects including buoyancy, shear, near-free-surface damping, bubbles, and Langmuir circulations complicate any attempt to describe turbulent motions. The variety of production mechanisms includes (but is not limited to) wind shear, wave breaking, internal gravity waves, double diffusion, and overturning due to the alternating heating and cooling of the ocean surface.
The numerical simulation of engineering-related problems in such an environment is a daunting prospect. For the case of wakes generated by ships and other man-made objects, the associated Reynolds number can be O (10 9 ) in the near field, while approaching very small Froude number in the far field. The broad separation of scales makes use of eddy-resolving methods at full scale prohibitively computationally expensive for use in design and analysis problems. In previous work, we have tested a set of Reynolds stress models (RSMs) against laboratory representations of the oceanographic environment (see Wall and Paterson [1], publication pending). In this work, we then further develop the application of these models to the problem of wakes, testing against laboratory and eddy-resolving model descriptions of stratified wakes. The models are then modified with source terms to produce a finite background turbulence, intended to model the environmental turbulence in the ocean.
As has been noted, one of the key complications associated with the ocean is the density stratification, which causes anisotropy in the stress tensor. In the late wake the effects of the stratification inevitably dominate the character of turbulent motion. Models which account for buoyancy-induced stress tensor anisotropy in some way are therefore desirable for simulations of late-wake behavior, necessitating a second-moment closure. There is an extensive history of second-moment approaches in modeling geophysical problems (early examples are the hierarchy of Mellor and Yamada [2] or the summary of Rodi [3]). Similar approaches have also been adopted in dealing with stratified wakes. Algebraic closures paired with two-equation models have been used to good effect by, for example, Hassid [4], Voropaeva et al. [5], or Voropaeva et al. [6]. Such models reproduce key behaviors of stratified wakes, including the vertical collapse and production of internal waves. Full stress-transport models have also been applied to the problem, with the earliest example being Lewellen et al. [7]. More recently, Voropaeva [8] has even adopted algebraic and transport-equation based triple-moment closures.
Due to the expectation of strong stratification effects in the late wake, models which remain realizable in an approximately two-component turbulence state are desirable. The family of cubic tensorial expansion models developed at the University of Manchester were developed with the two-component limit (TCL) as an explicit constraint, and have been applied to a variety of flows with strong buoyant forcing (for example, see Craft et al. [9], Craft and Launder [10], Suga [11], or Craft et al. [12]) The so called TCL model has also been applied to doubly-stratified (simultaneous salinity and temperature stratifications) environs, as presented by Armitage [13].
In evaluating RANS closures, it is possible to employ a wealth of LES and DNS numerical experiments which have been conducted to complement previous experimental studies of stratified wakes. Temporal-model, or 3D+t simulations such as those conducted by Dommermuth et al. [14] or Brucker and Sarkar (Brucker [15], Brucker and Sarkar [16]) have helped to describe the distribution of energy within the wakes of both towed and self-propelled bodies. More recently, body-inclusive simulations such as those conducted by Chongsiripinyo and Sarkar [17] have done much to refine the understanding the scaling laws which can be applied a given stratified wake, and to qualify each stage encountered during its life. These stages were originally identified as the three-dimensional (3D), non-equilibrium (NEQ), quasi-two dimensional (Q2D), and viscous three-dimensional (V3D) stages by Spedding [18]. It is an interesting digression to note that these stages roughly align with the stages of a full-scale ship wake as defined by Somero et al. [19] (the near wake, the far wake, and the persistent wake), though the rigorous definitions are different in each case. Some recent experiments and LES/DNS studies have also concerned themselves with the effect of non-trivial free-stream turbulence. Studies such as that of Amoura et al. [20] and Rind and Castro [21] have shown that environmental turbulent motions can have dramatic effects on wake behavior. In the case of stratified turbulence, the simulations of Radko and Lewis [22] include consideration of pre-existing double-diffusive fluctuations. The authors of that study also establish fairly simple wake-detection criteria based on the centerline deficits of dissipation rate and turbulent scalar variance (θ 2 or s 2 ). It is fully understood that much of the physics described by these eddy-resolving simulations will be lost in adopting a RANS approach; the results of these studies must then be carefully applied to refining RANS models.
Having addressed the scope of the problem, we now recapitulate a set of general criteria which must be satisfied by a turbulence model which might be applied to full-scale, far-field ship wakes in an oceanic environment: • The model must be implementable as part of a general-use computational fluid dynamics package (in the case of this work, the finite-volume code OpenFOAM) • The model must be able to accommodate the anisotropy that arises in stratified turbulent flows.
• The model must gracefully handle free-stream environmental turbulence, with differing levels of turbulent variance in both the temperature and salinity fields. Ideally, the impact of pre-existing turbulence on wake similarity would be accurately accounted for as well. • The model must reproduce classic stratified wake experiments such those of Lin and Pao [23] • The model must reproduce key behaviors observed using eddy-resolving methods, including decay rates of turbulence kinetic energy (TKE) and turbulence potential energy (TPE), and the growth (or lack thereof) of the wake in the vertical and horizontal directions.
The simulations and evaluation described for this study were conducted primarily to address the third and fifth bullet points.

Simulation Methodology
The model system of equations was solved using extensions to the open-source finite-volume fluid dynamics package OpenFOAM. A "2D + t" approach was adopted; the same as that employed by, for example, Lewellen et al. [24] or Voropaeva [8]. Mean field transport and RANS model equations were solved on a two-dimensional grid, representative of a slice of fluid through which the wake progenitor has passed.

Mean Transport Equations
Momentum was transported according to the Reynolds-averaged incompressible Navier-Stokes equations under the Boussinesq approximation: where U i is the mean-velocity vector, u i is the fluctuating component of velocity, P is the mean kinematic pressure, g i is the gravitational vector, ν is the fluid viscosity, ρ is the fluid density, and u i u j is the Reynolds stress tensor. For the laboratory-scale simulations conducted in this work, a linear equation of state for the density ρ was deemed sufficient: where the relevant scalar values are the mean temperature Θ and the mean salinity S, the 0 subscript denotes a reference value, and the expansion coefficients are defined by: Note that all of the simulations presented in the Section 3 are singly-stratified. In keeping with the methods employed for experimental study of stratified wakes, a salinity stratification was employed. The equation of state reduces to: For the remainder of the work only the model for the transport of salinity will be provided. However, the same model can be applied as temperature field under certain conditions. The scalar quantities are transported according to an advection-diffusion equation: where the total mean scalar field S is assumed to be the sum of a background S B and a perturbation to that background δS, and su i is the turbulent flux of the scalar quantity. The quantities u i u j and su i are supplied by the turbulence model.

Stress/Flux Transport
In general for this work, the framework and nomenclature for second moment models laid out by Hanjalić and Launder [25] is adopted, where mean quantities are denoted by capital symbols (U, S, etc.), fluctuating quantities by lower case symbols (u, s) and averaging is denoted by an over-bar (u i u j , su i , θ 2 ). The stress tensor can be obtained by solving the associated transport equation: where the terms on the right-hand side are the dissipation tensor ij , the shear production: the production due to the buoyancy body force in a Boussinesq fluid: the re-distributive effects due to pressure interactions: and diffusive effects: A free-stream sustaining source P ij ∞ is also included, intended to maintain the TKE at some finite value (preferably associated with some background condition representative of the active ocean environment). The forms of the free-stream sources employed are given in Subsection 8 Similarly, the transport of the turbulent flux of a scalar quantity, such as the temperature or salinity in ocean water, was described according to the equation: where the physical interpretation of each term is the same as given for (7). The source terms are: The terms ij , Φ ij , D ij , Φ sj , D sj , and the quantity s 2 must be modeled in order to close the second-moment RSTM. The following sections describe the closure approaches employed; the different models so constructed are summarized in Table 1.

Diffusive Process Closure
For all of the RSMs used in this work, the generalized gradient-diffusion hypothesis (GGDH) model originally proposed by Daly and Harlow [26] was employed to approximate the diffusive effects: Other, more complex closures for these terms have been applied stratified problems. The most pertinent example is the set of models employed by Voropaeva et al. [6], using both complex empirical algebraic expressions and even transport equations for a subset of the triple correlations. Craft and Launder [10] also recommends using transport equations which account for buoyancy effects on a subset of the triple correlations for strongly stratified flows. However, the effects described in that work were primarily associated with sharp pycnoclines, in which inhomogeneity in the triple correlations became significant. Such approaches would then likely be unneccessary for the linear-stratification environment in this work. Under certain environmental conditions, the diffusion closure may need further revision.

Pressure Strain/Scrambling Closure
In modeling the pressure strain and scrambling terms, it is common to adopt the approach of Chou [27], in which the pressure fluctuations are eliminated from (11) by taking the divergence of the transport equation for u i and so obtaining a Poisson equation for p. The details of such a derivation are here elided. The resulting expressions for Φ ij and Φ si can be arranged into terms associated with different physical processes: the return to isotropy (Φ 1 ), isotropization of mean strain (Φ 2 ), isotropization of body forcing (Φ 3 ), and wall blocking effects (Φ w ). The wall blocking effects are neglected for this work, as the flow of interest is a free shear flow. The term (11) and analogous term from (13) can then be written as: Typically, models make use of the stress anisotropy tensor a ij and its invariants; the associated definitions are included for completeness: Lumley's flatness parameter is also employed by some models: Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 16 November 2020 doi:10.20944/preprints202011.0432.v1 which takes the value of unity in isotropic turbulence, and the value of zero in two-component turbulence. Additionally, some models are written in terms of the symmetric (S) and an antisymmetric (W) portions of the velocity gradient tensor: Two pressure-strain models were employed for the simulations presented in this work. The simpler was a linear model, and the other a cubic model based on the work of Craft et al. [9]. The linear model employed the return-to-isotropy model first proposed by Rotta [28], and the linear isotropization-of-production terms from Launder et al. [29]: The accompanying model of the pressure-scrambling processes in the scalar flux equations is detailed by Gibson and Launder [30]: The cubic model was that developed at the University of Manchester, and is detailed in the book by Hanjalić and Launder [25]. The model is designed to be realizable approaching the so-called two-component limit (TCL), at which one of the normal stresses is approximately zero. The pressure-strain process models were originally described by Craft et al. [9], and have here been re-cast in terms of the symmetric and antisymmetric portions of the velocity gradient tensor: The pressure-scrambling process models were further developed by Craft and Launder [10], and for a doubly-stratified system (such as the temperature/salinity stratification in the ocean) by Armitage [13]: Notably, the coefficients in (35), (37), and (38) are not emperical, and are determined only by the realizability constraints. The principle justification for adopting such a complex model is that, even for a wake with an initially high Re and Fr, the flow will eventually decay to the point at which the turbulent Froude number Fr T = /kN is small. The turbulence will then be dominated by stratification effects, and so approach the two-component limit. Recent LES/DNS studies such as those by Chongsiripinyo and Sarkar ([31], [17]) have lent some credence to this notion, showing that the vertical normal stress decreases much more quickly than the horizontal normal stresses.

Scale-Equation Closure
The scale-determining equation is typically the most empirical part of a given RANS turbulence model. In describing stratified flows, a number different quantities have been used to describe the scale of turbulent motion; the various transport equations for kL, , or ω each have virtues, but are ultimately somewhat interchangable (see Umlauf and Burchard [32]). For this work, the empirical model of the equation developed by Craft et al. [9] was adopted: As with the stress transport equation, a free-stream source P ∞ is included. For the stress transport models, the coefficients for the model transport equation (39) were taken to be: The model for given by (39) and (40) is designed for free-shear flows, and has been found to be ill-posed for homogenous turbulence. It is therefore likely to be less-suited to the far-wake than to near-wake regions. Further, Pereira and Rocha [33] has noted a general deficiency in models like (39) in the case of strongly-stratified turbulence. The empirical model equation is therefore perhaps the best target for model improvement in future work. Two different models for the dissipation rate tensor were tested. The first assumes that ij is isotropic: The second is the model of Hallbäck et al. [34], and adopts a nonlinear dependence on the anisotropy of the stress tensor:

Scalar-Variance Closure
Closure of (13) also requires the variance of the scalar fluctuations, s 2 . Per Radko and Lewis [22], s 2 is also a potentially useful quantity in its own right. A transport equation can be solved for s 2 : which includes a free-stream source P ss ∞ . Dissipation of s 2 was modeled using the algebraic expression of Craft et al. [9]: Some preliminary simulations conducted using a transport equation for ss , rather than (44), indicated that the algebraic expression was sufficient. However, if environmental conditions are expected to have significantly different scalar and mechanical time scales (fossilized turbulence, for example), this evaluation may need revision.

Eddy Viscosity Model
For the sake of comparison, an isotropic eddy-viscosity model was also applied the same set of wake conditions. A standard k − model was employed, with the body-force effects on the scale equation modeled after the approach of Henkes et al. [35]: where the remaining modeled values are closed with the Boussinesq approximation:

Environmental Turbulence Sources
Finally, in order to accommodate the existence of environmentally generated background turbulence, free-stream source terms were introduced to the turbulence quantity transport equations, as proposed by Spalart and Rumsey [36]. The terms maintain the background turbulence quantities at a specified value, and have the added benefit of improving numerical stability and convergence. For the u i u j , the simplest approach is to introduce an isotropic source (the ∞ subscript denotes a free-stream value): The sources for the , s 2 equations, respectively, are then: A second, nearly two-component anisotropic source was also implemented, to test the potential impact of free-stream anisotropy: Note that, as the dissipation of the scalar flux su i is typically included in the pressure-scrambling model, and any of the flux-vector components may potentially be negative, such source terms for the flux-transport equations were not applied. Unless otherwise stated, in the simulations conducted for this work, the free-stream values k ∞ ∞ were chosen such that ν t ∞ = C µ k 2 ∞ / ∞ ≈ 0.5ν For the simulations presented in this work, these forms were employed. Further study may be needed to select forms which properly account for the stratification and background anisotropy.

Simulation Approach
As indicated in the introduction, the models were implemented for the open-source finite-volume code OpenFOAM. For all of the simulations detailed here, a second-order backward numerical differentiation scheme was applied for the temporal derivatives. Second-order linear schemes were applied for all spatial derivatives, with limiters applied as necessary to ensure convergence. The momentum equation (1) and continuity equation (2) were coupled using the widely employed hybrid PISO-SIMPLE (PIMPLE) algorithm, modified to include a Boussinesq body force after the fashion of Issa and Oliveira [37]. Note that, under the 2D + t approach employed, the axial component of the velocity U 1 given by (1) is actually the velocity difference: where U B is the propagation velocity of the wake generator. The mean velocity and turbulence kinetic energy were initialized according to idealized models of drag and self-propelled (or net-zero momentum, NZM) wakes. The wakes were initially axially symmetric, and the Reynolds stress tensor was initially isotropic. For the drag wake the expressions for velocity and TKE, in terms of radial position r, were: where the initial centerline values are dependent on the case in question. The form of (57) produces a TKE distribution roughly in line with the sphere wake measurements of Uberoi and Freymuth [38]. For the self-propelled (NZM) wake the expressions for velocity and TKE were: The form of (59) produces a TKE distribution which roughly resembles the measurements in the wake of the disk/jet wake of Naudascher [39]. The expressions for both wake types are plotted as a function of radial position in Figure 1. Note that the expression for the self-propelled wake produces a smaller total amount of TKE than that of the drag wake. For all cases, the dissipation rate was set such that the turbulent Reynolds number Re T = k 2 /ν had a constant value of 10000 throughout the wake. In studies employing eddy-resolving methods (e.g. Dommermuth et al. [14] or Brucker and Sarkar [16]), it was found that initializing fluctuations in the scalar value did not substantially change the behavior of the wake. Likewise, employing algebraic expressions to initialize su i and s 2 was found to have little effect on the RSM predictions. The scalar-associated turbulence quantities were therefore simply initialized to zero for this study. The simulation domain consisted of a square two-dimensional grid, 120D in both vertical and horizontal extent. The pressure field p was given a Dirichlet boundary condition, with a fixed value of zero. The other flow variables were given mixed Dirichlet/Neumann boundary conditions, dependent on the flux of the quantity at the boundary.
For all of the simulations conducted in this study, the initial wake diameter was given a dimensional value of D = 1 m. The gravitational vector was aligned with the x 3 axis (g 3 = −9.81m/s). The body propagation velocity U B was chosen to obtain the desired Reynolds number Re = U B D/ν. The background salinity stratification ∂S B /∂x 3 was then set to provide the buoyancy frequency required to match a given internal Froude number Fr = U B /ND, where the buoyancy frequency is given by N = −(g/ρ 0 )(∂ρ/∂x 3 ).

Results and Discussion
The employed model variations and the associated equations are summarized in Table 1. RSM1 is a stress transport model employing the realizable, two-component limit pressure-strain model. The variant RSM1a employs an anisotropic model of the dissipation rate tensor ij . The variant RSM1b employs an anisotropic free-stream turbulence source. RSM2 is a a stress transport model employing the simpler linear pressure-strain model, and EVM1 is the eddy-viscosity model included for comparison.
The conditions simulated are summarized in Table 2. The eddy-resolving simulation studies of Brucker and Sarkar [16] and Dommermuth et al. [14] were selected for comparison primarily because of the use of the "temporal model", which is more analogous to the 2D + t approach than a body-inclusive simulation.

Turbulence Decay
The decay of the root-mean-square vertical-velocity fluctuations (the square-root of the vertical normal stress u 3 u 3 ) and root-mean-square scalar fluctuations (the square-root of the scalar variance s 2 ) along the wake centerline are depicted in Figure 2 and Figure 3, respectively. The stress transport models achieve the correct decay rate, matching the (Nt) −1 exponential decay measured for the experiments detailed by Lin and Pao [23]. The rate is captured for both u 3 u 3 and s 2 . The differences between RSM1 and RSM2 are trivial for this case. The eddy-viscosity model was not paired with a transport equation for s 2 , and so is omitted from Figure 3.
While reproducing the decay rate, the RSMs systematically under-predict the magnitude of the scalar variance (both the peak value and the value during the decay). A partial explanation of the deficiency may be the uncertainty in the initial conditions for the problem. As noted in the previous section, the scalar-associated turbulence quantities were initially uniformly zero. While computationally convenient this is clearly non-physical, as mixing of the scalar quantities begins at the onset of turbulence in the body boundary layer, not at a finite downstream distance.
The eddy-viscosity model (EVM1) predicts a much too rapid decay, likely indicating that the expression for the coefficient c 3 from Henkes et al. [35] is poorly tuned for this particular problem. The rapid extinction of turbulence quantities leads to an unrealistically high preserved wake momentum in the later stages of the wake, as will be discussed in Subsection 2.
Finally, the time-dependent behavior of the centerline value of is given for the self-propelled case under the same conditions as the LES of Brucker and Sarkar [16] as Figure 4. The RSMs predict a decay rate in keeping with the LES, while the EVM again predicts a much too rapid decay. The introduction of an anisotropic dissipation rate tensor for RSM1a dos not produce a significant difference in behavior, despite the dependence of equation coefficients on the stress anisotropy for the RSMs.

Wake Momentum Decay
Model predictions of the decay of the mean defect velocity are compared to the eddy-resolving simulations of Brucker and Sarkar [16]. Figure 5 shows the comparison for the drag wake. The models all correctly predict the prolonged duration of the momentum wake due to the suppression of turbulent mixing by stratification. As noted in Subsection 1, the rapid extinction of turbulence quantities for EVM1 results in that model predicting a much larger sustained centerline defect velocity.    The RSMs in general reproduce the overall velocity decay well. The RANs models under-predict relative to the eddy-resolving model in an approximate region between Nt ≈ 6 and Nt ≈ 70. This roughly corresponds to the NEQ region of the wake, according to the stage breakdown suggested by Spedding [18]. Further analysis is required to determine if the disagreement with the eddy-resolving simulation can be explained by some physical mechanism occurring in this stage of the wake. As with the model predictions of the turbulence decay, there is not a substantial difference between RSM1 and RSM2 for this metric. Figure 6 shows the comparison for a self-propelled (NZM) wake. The figure shows both the peak value of the momentum associated with the thrust portion of the wake (U + s ) and the value associated with the drag portion of the wake (U − s ). As with the drag wake, the preservation of the wake momentum to late Nt is reproduced by the RANs models. EVM1 again predicts a too-high preserved mean momentum.
The RSM predictions for both the thrust and drag portions of the wake are in fair agreement with the eddy-resolving model for the NZM condition as well. The under-prediction in the NEQ stage is not present for the self-propelled case. The differences between RSM1 and RSM2 are somewhat more pronounced. RSM1 predicts a slight increase in U − s near Nt = 5, however, this ultimately puts that model's prediction more in line with the eddy-resolved simulation predictions. Finally, it is notable that the adoption of the anisotropic dissipation rate expression in model RSM1a does not produce any qualitative differences in behavior, and appears to reduce the model's agreement with the eddy-resolving simulation.

Wake Dimensions
In evaluating model prediction of wake dimensions, the general definition of wake height/width suggested by Brucker and Sarkar [16] is adopted, which uses the second central spatial moment of the square of the mean streamwise velocity: The integrated expression for the momentum width or height allows for direct comparison between the drag and self-propelled cases. Figure 7 shows the model predicted wake dimensions for the pure drag case. The RANS model predictions of the wake growth rate in horizontal direction roughly agree with the eddy-resolving simulation (disregarding the eddy-viscosity model). However, after approximately Nt = 100, the RSMs predict a slow in horizontal growth which is not observed in the LES results. However, there is substantial disagreement in the predictions of the vertical growth of the wake. Both the RANs and LES approaches predict a local peak in R 3 shortly after Nt = 1. However, the eddy-resolving simulation of Brucker and Sarkar [16] predicts a wake which shrinks in the vertical axis over most of the wake's lifetime, while the RANS models predict a small but positive growth rate. The discrepancy is difficult to explain, and further study is needed to determine the cause of the qualitative difference in behavior. Figure 8 shows the model predicted wake dimensions for the self-propelled case. The RSM models systematically under-predict the wake width of the self-propelled wake in comparison with the LES, however, the growth rate is in approximate agreement over a portion of the wake lifetime. In contrast with the drag wake, the wake height predictions agree fairly well with the LES, with both the RSMs and the LES indicating a wake which maintains a roughly constant height R 3 /D ≈ 0.95 at late Nt. The stronger agreement with LES suggests that the discrepancy seen in Figure 7 may be due to poor initial conditions in the drag wake simulations, rather tha a deficiency in the RSMs themselves. For both the drag and self-propelled cases, the differences between the predictions of RSM1, RSM1a, and RSM2 are again mostly trivial.

Spatial Energy Distribution
The spatial distribution of the energy of a wake also be examined; again for comparison with the LES simulations of Brucker and Sarkar [16]. For the drag wake, Figure 9 shows a slice of the domain with the local mean kinetic energy (MKE), while Figure 10 supplies the same for the turbulent kinetic energy. Likewise, Figure 11 and Figure 12 show the MKE and TKE distributions for the self-propelled wake. The predicted MKE and TKE distributions are vertically symmetric.
The drag wake MKE shows a wake which has grown in the horizontal direction, while growth in the vertical is suppressed, in keeping with the thicknesses measured in Subsection 3. Figure 10 shows that at late Nt the TKE has separated into two peaks with a saddle point on the centerline, which is broadly in agreement with the behavior predicted by the LES of Brucker and Sarkar [16]. The primary difference between RSM1 and RSM2 is a slightly larger peak TKE value for the former. Examining Figure 11, the self-propelled wake possesses two distinct regions of mean kinetic energy. The thrust portion of the wake is still concentrated at the centerline, while the drag portion has been seperated into two "lobes" roughly one diameter above and below the centerline. The distribution of MKE compares favorably with the LES of Brucker and Sarkar [16], which predicted similar lobes in similar locations. Figure 12 illustrates perhaps the most significant difference in behavior between RSM1 and RSM2 observed in this study. The former predicts a core of TKE for the self-propelled case, while the latter predicts two separate peaks (similar to the behavior of the drag wake). The prediction of RSM1 is qualitatively more similar to the distribution observed by Brucker and Sarkar [16] for the self-propelled case.

Collapse-Induced Internal Gravity Waves
The internal gravity waves (IGWs) produced by the vertical collapse of the wake can be examined taking slices in the (x 2 − x 3 ) and (x 1 − x 3 ) planes (note that the x 1 direction is taken to be related to t by the body speed for the type of simulation conducted in this study, i.e. x 1 = U B t). Figure 13 shows the perturbation to the salinity field (which, in this case, is equivalent to showing the density perturbation) as a function of time and vertical location. The wake produces waves which propagate upward and downward through the linear stratification, with an oscillation period roughly corresponding to the buoyancy period 2π/N. Figures (14) and (15) show slices in the x 2 − x 3 plane, depicting waves which radiate from the wake centerline. The number of rays increases with increasing Nt, or equivalently, with downstream distance. It is important to note that both temporal model LES and 2D + t RANs models omit body-generated lee waves. Figure 16 shows the time variation of the turbulent kinetic and turbulent potential energy (TKE,TPE), integrated over the entire domain. The TKE is separated into vertical (VTKE) and horizontal (HTKE) components. The case is a drag wake under the same conditions as the LES of Dommermuth et al. [14], and the behavior depicted may be qualitatively compared with that study. The RSM predicts a decay rate roughly in line with the (−2/3) exponential decay measured in the LES over a portion of the wake's lifetime. Additionally, the vertical TKE and TPE oscillate with a period roughly equally to the buoyancy period 2π/N, exchanging energy as they do so. The oscillatory behavior is also in line with the behavior observed by Dommermuth et al. [14].

Free-stream Turbulence Effects
Finally, the effect of increased free-stream turbulence intensity may be briefly explored by increasing the strength of the sustaining sources included in the turbulence models. As indicated in Table 2, a second set of simulations was conducted with a high free-stream TI for the self-propelled case (case BS1a). The test was conducted for both an isotropic free-stream u i u j source, and the anisotropic source given by (54). Figure 17 shows the decay of the mean velocity. Strong background turbulence predictably increases the rate at which the mean velocity decays. The anisotropic source term appears to produce a stronger effect for the same free-stream TI. Figure 18 depicts the predicted wake dimensions under the same conditions. In this case, the sustaining sources are strong enough to overcome the buoyancy effects, and the wake grows in both the vertical and horizontal direction. The growth continues past the point at which the wake height ceases growth in quiescent conditions. By late Nt, the wake turbulence has been reduced to the background value, and turbulent transport of the wake momentum is exclusively by the background turbulence.
While the predicted behavior shown in Figure 17 and Figure 18 makes intuitive sense, further study is required to confirm the efficacy of the free-stream source approach employed.

CONCLUSION
In this work, we have demonstrated the use of a pair of anisotropic stress-transport RANS turbulence models intended to be used to simulate full-scale wakes in an active ocean environment. The models were found to reproduce a number of important stratified wake behaviors as observed in LES and laboratory studies. In particular, the models capture the decay rates of key turbulence quantities, the preservation of mean momentum to late Nt due to suppression of turbulence, and the wake collapse and internal gravity wave production. It was found that the more complex TCL pressure-strain based RSM did not differ significantly in behavior from the simpler linear model under the conditions simulated. Likewise, the use of an anisotropic dissipation-rate tensor for the stress transport equations did not substantially improve agreement with LES model predictions. Use of these more complex models required approximately 10%-20% more computing hours over the linear RSM for a given wake case. Further study at late Nt is needed to determine if the TCL model's additional cost is justified for low turbulence Froude numbers. Finally, the models were further modified with additional source terms to supply a nonzero background turbulence, which was found to increase the rate of wake decay and increase wake thickness growth. More tests are needed to carefully validate this capability.