CFD Modeling of E ﬄ uent Discharges: A Review of Past Numerical Studies

: E ﬄ uent discharge mixing and dispersion have been studied for many decades. Studies began with experimental investigations of geometrical and concentration characteristics of the jets in the near-ﬁeld zone. More robust experiments were performed using Laser-Induced Fluorescence (LIF) and Particle Image Velocimetry (PIV) systems starting in the 20th century, which led to more accurate measurement and analysis of jet behavior. The advancement of computing systems over the past two decades has led to the development of various numerical methods, which have been implemented in Computational Fluid Dynamics (CFD) codes to predict ﬂuid motion and characteristics. Numerical modeling of mixing and dispersion is increasingly preferred over laboratory experiments of e ﬄ uent discharges in both academia and industry. More computational resources and e ﬃ cient numerical schemes have helped increase the popularity of using CFD models in jet and plume modeling. Numerous models have been developed over time, each with di ﬀ erent capabilities to facilitate the investigation of all aspects of e ﬄ uent discharges. Among these, Reynolds-averaged Navier-Stokes (RANS) and Large Eddy Simulations (LES) are at present the most popular CFD models employing e ﬄ uent discharge modeling. This paper reviews state-of-the-art numerical modeling studies for di ﬀ erent types and conﬁgurations of discharges, including positively and negatively buoyant discharges, which have mostly been completed over the past two decades. The numerical results of these studies are summarized and critically discussed in this review. Various aspects related to the impact of turbulence models, such as k- ε and Launder-Reece-Rodi (LRR) models, are reviewed herein. RANS and LES models are reviewed, and implications for the simulation of jet and plume mixing are discussed to develop a reference for future researchers performing numerical investigations on jet mixing and dispersion. aimed to discover whether the RANS models were suitable for predicting the velocity and concentration properties of wall jets. They also aimed to identify the best performing turbulence model. The jets were modeled with OpenFOAM, and those turbulence models were compared to those from the numerical studies of [58,59] ( Fr d = 11.61 to 42.33). The linear (standard k- ε , RNG k- ε , realizable k- ε and SST k- ω ) and Reynolds stress (Launder-Gibson and LRR) turbulence models were tested in their study. This was the ﬁrst study for wall jet modeling that compared several turbulence models.


Introduction
Increasing population and industrial growth have meant a global increase in effluent discharge in water bodies. Industries located near shore that discharge significant volumes of thermal and saline effluent include desalination plants, mining operations, water treatment plants, and nuclear power plants.
The direct discharge of wastewater into lakes, rivers, and seas can increase turbidity and change the ambient temperature [1]. Salinity is also a major public and scientific concern [2]. Coastal waters receive concentrated salt brine as discharge from seawater desalination plants (Figure 1), chemical waste from biofouling (e.g., chlorine), and fertilizers. The water bodies that receive the industrial discharge are often very sensitive environments [3,4], and the design of the discharge facility to disperse effluent and reduce the concentration of effluents may help protect the receiving water body. Ref. [5] noted that certain coastal ecological zones are particularly vulnerable to effluent discharges, including salt marshes, mangrove forests, coral reefs, and other low-energy intertidal areas. The Persian Gulf and the Red Sea are particularly sensitive to effluent due to their low hydro-dynamism.
Local fisheries, tourism industries, and other economic consequences are affected by the health of coastal environments. The mixing and dispersion qualities of jets were first noted in the 1950s. Length-scale and integral models were among the first numerical studies on jet and effluent mixing; however, these models are only briefly described here, as the focus of this paper is to review CFD studies on jet mixing.
Length-scale models use relationships calibrated based on experimental data to forecast the steady-state behavior of the effluent discharge. These models use non-dimensional numbers to categorize discharge regimes (see Section 2.1 for details) and are fast to process. However, they are sensitive to user error and Monte Carlo testing, and thus an experienced user may be required to interpret the data when the result is near the boundary of two regimes, as small changes induce large differences in these models [6]. As well, these models are often inaccurate when used on different parameters than those for which they were validated [7]. NRFIELD [8][9][10] and CORMIX [11] are two commonly used length-scale models.
Jet integral models, according to [6], solve mass and momentum conservation equations based on the assumptions that the velocity profiles of jets have no radial variation, and that the jet profile is axisymmetric and Gaussian. In the 1950s and 1960s, first-order jet integral models were proposed by [12,13] based on the jet entrainment closure approach and by [14,15] based on the jet diffusion approach. Refs. [16][17][18] developed second-order jet integral models. Visual Plumes [19] and CORMIX (CorJet, based on the integral model [18]) are two popular mixing tools based on jet integral models.
According to [6], the integral models are less reliable when there is any of the following: (i) the discharge's initial momentum and buoyancy acting in opposite directions, resulting in instabilities on the edge, as observed in the inner half (lower half) of inclined dense jets [20,21]; (ii) noticeable interaction between the mean flow and the jet, (iii) an unsteady mean ambient flow; (iv) a great effect due to horizontal or lateral boundaries [18]; (v) an unstable near-field area, with a re-entrainment of concentrated effluent into jet [11]; or, (vi) a large re-entrainment of concentrated effluent from mid-and far-fields into the near-field jet due to tidal cycles [22].
Jet flow modeling by the CFD tools is not perfect, but it is an improvement over the parameter-based jet integral and length scale models. Issues that remain with CFD tools include the following: (i) Water 2020, 12, 856 3 of 31 accuracy, (ii) stability, (iii) computational time, (iv) complicated codes that require expert knowledge to use them efficiently and accurately, and (v) simulations that need calibrating and validating.
Turbulent length-scale models are often resolved with a turbulence model to parameterize unresolved mixing and dispersion scales. One should apply turbulence models with caution, as they sometimes provide stable but unrealistic solutions, such as when they are applied to physical scenarios for which they have not been validated.
When using a CFD model, it can be a challenge to create and resolve the mesh and to define appropriate boundary conditions (e.g., intensity and the turbulence dissipation rate). A high mesh resolution is often needed for a stable solution, even when the turbulence model is a good match. This means that CFD modeling is numerically expensive. Even with current computing systems, accurate CFD models for near-field dispersion and mixing might need simulation times of several days or weeks. This is much more expensive compared to the parametric-based models that can produce results on the order of minutes and seconds [6]. There is a balance between model stability, numerical diffusion, mass and momentum conservation, boundedness, and computational cost. These choices can significantly influence the estimation of modeled concentration.
However, once built, calibrated, and validated, CFD models can produce high-resolution three-dimensional images of jet mixing and dynamics. CFD models are free from some of the assumptions that restrict integral models. Since CFD models do not require the assumption of a steady-state condition or self-similarity in the jet profile, they can include a variety of external effects such as the presence of surface waves and encompass a wide range of boundary conditions to allow users to directly simulate the boundary interaction.
CFD modeling of jet discharges has been approached in a variety of ways, including both hydrostatic and non-hydrostatic approaches to the Reynolds-averaged Navier-Stokes (RANS) and the Large Eddy Simulations (LES) models. Both models have functioned well over the past decade to simulate effluent discharges [23][24][25][26][27]. RANS models are based on a time-averaging method and result in a time-averaged mean velocity field, which is averaged over a longer time period than the time constant of the velocity fluctuations and results in a constant mean velocity without nuance for time-dependent variations. LES is based on filtering instead of averaging. A filter size is identified, and flow scales equal to or larger than this size are calculated exactly, and scales smaller than the filter size will be modeled. The smaller the filter size, the more concise the calculated time variation resolution of the velocity vectors. RANS models are more numerically efficient than LES models, while providing enough detail for engineering applications. Thus, they have become the most prevalent CFD models used for the design of outfall systems.
The Direct Numerical Simulation (DNS) method is less applicable to engineering problems, functioning more as a research tool. It is CPU-intensive, as it attempts to resolve Navier-Stokes equations with no approximation of the turbulence and requires a very fine numerical resolution to capture all the turbulence details. It basically resolves entire turbulence scales temporally and spatially. Mesh systems should be very fine to resolve all the spatial scales [28].
The hydrostatic pressure assumption in CFD models has been widely used to model shallow water systems in which the horizontal length scale of the problem is much larger than the vertical length scale. Ref. [29] stated that brine discharges from a desalination plant in the near field is a case in which the non-hydrostatic pressure approach should be adopted as the eddy viscosity and the acceleration terms in the momentum equation are comparable to that of the gravitational acceleration term for the vertical velocity component. Table 1 (after [7]) summarizes the existing modeling packages (commonly used in the industry) for the simulation of jet and plume mixing. The objectives of this review paper are to (i) provide an overview of the theoretical advancements on jet mixing modeling, (ii) review the current state-of-the-art research on CFD modeling of jets and plumes, (iii) investigate the most common turbulence models used in the jet mixing modeling, and (iv) indicate areas of research needs. These objectives are accomplished by completing a thorough literature review, discussing research gaps, and outlining the methodological solutions to address current and future research needs. This paper is organized in the following order: Section 2 presents the numerical details of jet studies, including dimensional analysis, that govern numerical equations. Section 3 reviews numerical studies on inclined dense jets. Numerical studies of vertical jets, horizontal jets, and surface discharges are reviewed in Sections 4-6, respectively. The discharge port configuration is discussed in Section 7. Finally, the conclusions complete the study.

Jet Studies: Details of Numerical Analysis
Most CFD codes for numerical jet studies employ the following theories.

Dimensional Analysis
This section introduces the dimensional analysis and fluid governing equations used in the numerical analysis of effluent discharges.
Water 2020, 12, 856 5 of 31 Figure 2 shows a sketch of an inclined dense jet with negative buoyancy in stagnant ambient water. The effluent discharges with an initial jet velocity U 0 , jet density ρ 0 , jet nozzle diameter D, angle θ to the horizontal, and ambient water density ρ a (with ρ 0 > ρ a ). As it is discharged, the jet reaches a maximum or terminal rise height, y t , and then falls due to the negative buoyancy, mixing with the ambient water. The jet lands on the seabed and then it spreads horizontally as a density current.
Most CFD codes for numerical jet studies employ the following theories.

Dimensional Analysis
This section introduces the dimensional analysis and fluid governing equations used in the numerical analysis of effluent discharges. Figure 2 shows a sketch of an inclined dense jet with negative buoyancy in stagnant ambient water. The effluent discharges with an initial jet velocity U0, jet density , jet nozzle diameter D, angle ϴ to the horizontal, and ambient water density (with > ). As it is discharged, the jet reaches a maximum or terminal rise height, yt, and then falls due to the negative buoyancy, mixing with the ambient water. The jet lands on the seabed and then it spreads horizontally as a density current. The dispersion of concentration relies on both ambient and jet characteristics, including U0, D, and ϴ, described above, plus the initial density difference ∆ = − , the jet discharge concentration, C0, the turbulence intensity of the jet, ITJ, and the ambient water depth Ha. The jet densimetric Froude number (Frd) is a key parameter for dense jet analysis. Frd is the ratio of inertia to buoyancy, and it is calculated with the reduced gravitational acceleration ( ) as follows: Inclined negatively buoyant jets are described in the equations that follow, using the buoyancy flux (B0), the jet discharge volume flux (Q0), and the kinematic momentum flux (M0): A characteristic length, such as the maximum terminal rise height (yt), for a jet may be written with the momentum and source length scales (LM and LQ), using dimensional analysis as follows: The return point dilution, Sr, and centerline peak dilution, Sm, can be expressed as: The dispersion of concentration relies on both ambient and jet characteristics, including U 0 , D, and θ, described above, plus the initial density difference ∆ρ 0 = ρ 0 − ρ a , the jet discharge concentration, C 0 , the turbulence intensity of the jet, I TJ , and the ambient water depth H a . The jet densimetric Froude number (Fr d ) is a key parameter for dense jet analysis. Fr d is the ratio of inertia to buoyancy, and it is calculated with the reduced gravitational acceleration (g 0 ) as follows: Inclined negatively buoyant jets are described in the equations that follow, using the buoyancy flux (B 0 ), the jet discharge volume flux (Q 0 ), and the kinematic momentum flux (M 0 ): A characteristic length, such as the maximum terminal rise height (y t ), for a jet may be written with the momentum and source length scales (L M and L Q ), using dimensional analysis as follows: The return point dilution, S r , and centerline peak dilution, S m , can be expressed as: Water 2020, 12, 856 6 of 31

CFD Governing Equations
The Navier-Stokes equations, as follows, govern the CFD mechanisms. The continuity equation governs the conservation of mass: Momentum is also conserved, as shown in the following equations: where t is time; u, v, and w are the mean velocity components in the x, y, and z directions, respectively; υ eff is effective kinematic viscosity (υ eff = υ t + υ); υ t is turbulent kinematic viscosity; P is fluid pressure; g is gravitational acceleration; ρ is fluid density; and ρ 0 is the reference fluid density. Dividing by density (ρ) and adding the term of buoyancy to the vertical momentum equation accounts for variable density effects. Then, the seawater equation of the state is used to calculate the density for both the jet and the ambient water [30].
The advection-diffusion equations are used to model the concentration and temperature changes over time.
where T is the fluid temperature, C is the fluid concentration (salinity), d is the diffusion coefficient, k eff is the heat transfer coefficient, Pr is the Prandtl number, and Pr t is the turbulent Prandtl number.

Discharge through Inclined Dense Jets
An inclined jet is discharged at an angle of 0 • to 90 • from the horizontal and may be positively or negatively buoyant. Most of the inclined jets used in engineering applications are negatively buoyant jets, which are also known as inclined dense jets. Most high-density discharge outfalls are designed as inclined jets due to their proven near-field dilution rate and efficiency.
Many new Computational Fluid Dynamics (CFD) techniques were developed in the 1970s, but detailed and reliable experimental data were difficult to come by. In the 1990s, Particle Image Velocimetry (PIV) and Laser-Induced Fluorescence (LIF) experimental techniques provided more reliable data, and these studies became benchmarks for CFD model validation. [31] and others performed some numerical studies on jet mixing in the 1990s, but [32] were the first to fully apply CFD tools to estimate the effluent discharge jet flow behavior. Vafeiadou simulated dense effluent discharges at 45 • , 60 • , 70 • , 80 • , and 90 • inclinations, by using the ANSYS CFX with a mesh of 400,000 elements and adopting the k-ω Shear Stress Transport (SST) turbulence model. Comparing Vafeiadou's results with experimental data by [33,34], it appeared that their numerical results agreed, supporting the further use of such numerical models to estimate the characteristics of dense effluent discharges. Vafeiadou's study pioneered the use of CFD for the mixing of effluent discharges into water bodies.
Ref. [35] also used ANSYS CFX and applied the standard k-ε turbulence model for the RANS equations closure. Experimental data of dense effluent discharges were simulated by the k-ε turbulence model using both a standard and adjusted turbulent Schmidt number of S ct = 0.9. Both of these simulations with the k-ε turbulence model were more accurate than the integral models and the analytical solutions. In the experimental studies, the buoyancy-induced instabilities were observed on the lower (inner) half of the jet. However, in the study of [35], the k-ε simulations (both the standard and the calibrated) underpredicted the spread of the jet and the integrated centerline dilution at the top rise of the jet, because they overestimated the influence of the stabilizing density gradients. This study only achieved qualitative comparisons and did not arrive at any conclusions about which model performed the best.
Ref. [24] used OpenFOAM to explore the numerical modeling of 30 • and 45 • inclined dense jets in calm ambient water conditions. OpenFOAM has a solid structure of various solvers, including a wide range of turbulence models that could be used to the jet mixing and dispersion analysis. This study focused on how the choice of the turbulence model can affect the jet trajectory and predicted dilution and the importance of turbulence closure in the Navier-Stokes equations. They applied five RANS turbulence models to investigate the accuracy of the CFD predictions: the Launder-Reece-Rodi (LRR) Reynolds stress model, the Launder-Gibson Reynolds stress model, the RNG k-ε linear eddy viscosity model, the realizable k-ε linear eddy viscosity model, and a non-linear k-ε eddy viscosity model. A summary of these turbulence models is discussed in [36].
Ref. [24] discretized the temporal term with a first-order implicit Euler scheme. The standard finite volume method with a Gaussian integration was used to discretize the advection-diffusion terms. The preconditional conjugate gradient (PCG) scheme was used to solve the pressure field, and the preconditioned biconjugate gradient (PBiCG) scheme was used for the other fields: U, T, C, k, ε and ω. This study focused on analyzing the inclined dense discharge, including the jet trajectories, jet terminal rise height, jet centerline peak, jet centerline, and jet horizontal return point (i.e., impact point). The modeled values for velocity and concentration profiles were compared to experimental data. The effects of the increased jet velocities of the high Froude numbers observed in field applications were investigated with a sensitivity analysis. Ref. [24] found that of the turbulence models tested in their study, the realizable k-ε and LRR models were more accurate than the other three tested. Tables 2 and 3 (after [24]) summarize the results for these two models (only presented 30 • jets for the sake of brevity).
The geometry, concentration, and velocity fields in the near field of effluent discharges can be characterized when the jet flow pattern is understood. Figure 3 presents the normalized concentration maps (C/C 0 , where C is the computed concentration and C 0 is the discharge concentration) and the dilution isolines for inclinations of 30 • and 45 • ; realizable k-ε and LRR turbulence models were used for each case.
Effluent discharge trajectory is an important factor in the design of ocean outfall systems. Discharge trajectory would basically identify the flow path that a jet would travel through until it impacts the bed. A dilution equation (Equation (15); C a is ambient concentration) was used to calculate the dilution values for plotting in Figure 3. The discharge trajectory and flow growth are also illustrated in Figure 3 as the jets travel downstream with clear differences between the two turbulence models. LRR is an anisotropic turbulence model and performs more reliably for calculating the shear forces on jet edges and providing more realistic predictions. Tables 2 and 3 summarize the discharge trajectory and flow growth compared to other turbulence models.
Water 2020, 12, 856 8 of 31  Table 3. Evaluation of performance of two turbulence models (realizable k-ε and LRR models) for 30 • inclined jet (Source: [24]).    The jet centerline often follows the maximum cross-sectional velocity or concentration along the flow path perpendicular to the discharge trajectory. The best way to extract the centerline is to start at the nozzle and create a velocity vector map. Ref. [40] noted that the maximum concentration and velocity profiles almost coincide, but the maximum concentration profile usually decreases more quickly than the velocity. The momentum and the buoyancy-induced instabilities also affect the trajectories of the inclined dense discharges. While the discharge rises close to the nozzle, negative buoyancy forces affect the upward momentum and decrease it until they dominate the concentration transport (somewhere after the discharge rise peak) and the maximum concentration profile sinks down toward the bed faster than the velocity maximum profile. The centerline results in [24] were extracted following the maximum cross-sectional velocity similar to [40].
Ref. [25] used the buoyancy modified standard k-ε model and the two best performing models in [24] (realizable k-ε and LRR) for inclined dense jets. They modified the standard k-ε turbulence scheme to include the standard Boussinesq gradient diffusion hypothesis (SGDH) and the general gradient diffusion hypothesis (GGDH). The governing equations for the k-ε turbulence model that includes the buoyancy term are as follows: where the shear production term P and the buoyancy production term G in the k-ε turbulence model in Equation (16) are calculated as follows: The way the term ρ u i in Equation (19) is solved defines the SGDH and GGDH methods. The corresponding equations for SGDH and GGDH are as follows, respectively.
The modified SGDH and GGDH k-ε turbulence models were used in OpenFOAM to explore the influence of the buoyancy term, using turbulence coefficients C µ = 0.09, C 1ε = 1.44, C 2ε = 1.92, and a calibrated coefficient, C 3ε for sensitivity tests with C 3ε = 0.9, 0.6, and 0.4. Figure 4 shows numerical and experimental results for a 45 • inclined dense jet on the central plane. The standard k-ε results are relatively close to the data by experiments, both with and without modifications ( Figure 4a). Adding the buoyancy term to the turbulence model causes the jet to spread more along the inner (lower) half of jet, where stronger buoyancy-induced forces would be experienced. When the buoyancy effects are set to be stronger (i.e., with a smaller C 3ε ), the discharge growth rate is larger, and the numerical results align more closely to the data by experiments (not presented in this paper). The SGDH and GGDH methods showed similar results when the computed and experimental discharge trajectories were compared. This could be attributed to the small difference in density between the discharge and the receiving water (smaller than 1%). The buoyancy effect on the effluent discharge mixing characteristics should be investigated further, especially for larger density differences between the discharge and ambient water (larger than 1% in this study) as well as C 3ε calibration for mixing applications. Water 2020, 12, x FOR PEER REVIEW 8 of 30 Ref. [25] further studied the concentration and velocity effects on the dispersion properties and the discharge growth rate for the 30° and 45° inclinations ( Figure 5). The figure presents the normalized cross-sectional concentrations for the 30° and 45° discharges using the realizable k-ε and LRR turbulence models. The realizable k-ε model predicts a narrower jet width than the LRR model does. The jets predicted by Linear Eddy Viscosity Models (LEVMs) are usually thinner than the ones predicted using the Reynolds Stress Models (RSMs), due to the assumption of isotropic turbulence in the LEVMs. The improved behavior of the LRR model could be due to stress anisotropy.
LEVMs are less expensive numerical models than the RSMs. There are reasons that justify extra terms in the RSMs over the two-equation models such as k-ε for the mixing application. To begin with, Reynolds stresses are affected by the pressure strain term through both the turbulent fluctuations (which is related to the gradient in mean velocity) and to the turbulence-turbulence interaction (eddies effects in the fluctuation of pressure), which do not relate directly to the mean flow changes. The first process is called rapid pressure, and the second is known as slow pressure.
This affects the energy redistribution between stress components, which is important in estimating the degree of the stresses' anisotropy. Inclined dense discharges disperse according to the interaction of eddies. Therefore, both the rapid pressure and slow pressure terms may be important. It is possible that the reason the LRR model works well is that it includes both the slow and rapid pressure terms. Ref. [25] further studied the concentration and velocity effects on the dispersion properties and the discharge growth rate for the 30 • and 45 • inclinations ( Figure 5). The figure presents the normalized cross-sectional concentrations for the 30 • and 45 • discharges using the realizable k-ε and LRR turbulence models. The realizable k-ε model predicts a narrower jet width than the LRR model does. The jets predicted by Linear Eddy Viscosity Models (LEVMs) are usually thinner than the ones predicted using the Reynolds Stress Models (RSMs), due to the assumption of isotropic turbulence in the LEVMs. The improved behavior of the LRR model could be due to stress anisotropy.
LEVMs are less expensive numerical models than the RSMs. There are reasons that justify extra terms in the RSMs over the two-equation models such as k-ε for the mixing application. To begin with, Reynolds stresses are affected by the pressure strain term through both the turbulent fluctuations (which is related to the gradient in mean velocity) and to the turbulence-turbulence interaction (eddies effects in the fluctuation of pressure), which do not relate directly to the mean flow changes. The first process is called rapid pressure, and the second is known as slow pressure.
This affects the energy redistribution between stress components, which is important in estimating the degree of the stresses' anisotropy. Inclined dense discharges disperse according to the interaction of eddies. Therefore, both the rapid pressure and slow pressure terms may be important. It is possible that the reason the LRR model works well is that it includes both the slow and rapid pressure terms. OpenFOAM was used to simulate 45° inclined discharges using the LES method and the Smagorinsky and Dynamic Smagorinsky sub-grid scale (SGS) by [26]. This resulted in numerical predictions including geometrical characteristics, the jet trajectory, jet spread, and eddy structures. Eddy sizes in the LES method are based on the grid spacing. In LES models, large eddies are simulated directly by computing the Navier-Stokes equations, but small eddies are simulated using, for instance, Boussinesq hypothesis assumptions. The simulations were run with the twoLiquidMixingFoam solver in OpenFOAM. Most of the parameters in Table 2 were used to compare [26]'s study results with other experimental and numerical studies. Eddy structures and turbulence characteristics were computed utilizing the experimental data of [41][42][43] for comparison. The coherent structure of the inclined dense discharge plays an important role in the development of flow in the jet and can help elucidate the mixing properties and the intensity of turbulence.
LES seemed to underestimate the turbulence intensity in the area of buoyancy-induced instabilities (i.e., lower half of the jet), as reported by [26]. For turbulent flows, the effect of the small scales on the simulated resolved scales are accounted for with an SGS approach. To complete the closure, for the Smagorinsky model, the parameter (Cs) in the eddy viscosity, = ( △) , where Δ is the LES filter size and Sij is the strain rate tensor, needs to be specified. The accuracy in choosing this parameter is directly related to the kinetic energy dissipation.
In 1991, [44] made an important addition to turbulence theory and modeling when they proposed a model for the dynamic determination of Cs. Rather than obtaining the minimal coefficient from predetermined expressions, the model analyzed large-scale turbulences during the numerical simulation to deduce them. This dynamic model approach was proved to be more applicable in the cases with complex interactions.
Ref. [6] predicted the mixing and dispersion of inclined dense discharges in stagnant surroundings by developing a three-dimensional CFD model that uses the CFD tool named Fluidity [45]. Fluidity uses a range of control volumes and finite element discretization to solve the 3D Navier-Stokes equations using unstructured mesh; it is capable of including Coriolis effects, density variations, turbulence models, tidal forcing, and associated buoyancy forces. Fluidity includes an anisotropic adaptive mesh capability that controls the solution accuracy locally throughout the model domain. The standard k-ε and V-LES (very-large eddy simulation) models in the dense jet simulations were used to close the Navier-Stokes equations. Although they only studied the 60˚ jet for the centerline rise height, impact point, and its dilution, and though their study was preliminary and not detailed in terms of model validation and comparisons, it seems that their model underpredicted the terminal rise height, the distance from the nozzle to the impact point, and the OpenFOAM was used to simulate 45 • inclined discharges using the LES method and the Smagorinsky and Dynamic Smagorinsky sub-grid scale (SGS) by [26]. This resulted in numerical predictions including geometrical characteristics, the jet trajectory, jet spread, and eddy structures. Eddy sizes in the LES method are based on the grid spacing. In LES models, large eddies are simulated directly by computing the Navier-Stokes equations, but small eddies are simulated using, for instance, Boussinesq hypothesis assumptions. The simulations were run with the twoLiquidMixingFoam solver in OpenFOAM. Most of the parameters in Table 2 were used to compare [26]'s study results with other experimental and numerical studies. Eddy structures and turbulence characteristics were computed utilizing the experimental data of [41][42][43] for comparison. The coherent structure of the inclined dense discharge plays an important role in the development of flow in the jet and can help elucidate the mixing properties and the intensity of turbulence.
LES seemed to underestimate the turbulence intensity in the area of buoyancy-induced instabilities (i.e., lower half of the jet), as reported by [26]. For turbulent flows, the effect of the small scales on the simulated resolved scales are accounted for with an SGS approach. To complete the closure, for the Smagorinsky model, the parameter (C s ) in the eddy viscosity, υ t = ρ(C s ) 2 S ij , where ∆ is the LES filter size and S ij is the strain rate tensor, needs to be specified. The accuracy in choosing this parameter is directly related to the kinetic energy dissipation.
In 1991, [44] made an important addition to turbulence theory and modeling when they proposed a model for the dynamic determination of C s . Rather than obtaining the minimal coefficient from predetermined expressions, the model analyzed large-scale turbulences during the numerical simulation to deduce them. This dynamic model approach was proved to be more applicable in the cases with complex interactions.
Ref. [6] predicted the mixing and dispersion of inclined dense discharges in stagnant surroundings by developing a three-dimensional CFD model that uses the CFD tool named Fluidity [45]. Fluidity uses a range of control volumes and finite element discretization to solve the 3D Navier-Stokes equations using unstructured mesh; it is capable of including Coriolis effects, density variations, turbulence models, tidal forcing, and associated buoyancy forces. Fluidity includes an anisotropic adaptive mesh capability that controls the solution accuracy locally throughout the model domain. The standard k-ε and V-LES (very-large eddy simulation) models in the dense jet simulations were used to close the Navier-Stokes equations. Although they only studied the 60 • jet for the centerline rise height, impact point, and its dilution, and though their study was preliminary and not detailed in terms of model validation and comparisons, it seems that their model underpredicted the terminal rise height, the distance from the nozzle to the impact point, and the minimum dilution at impact point compared to the experimental data. This could be due to little water entrainment toward the jet centerline; therefore, the simulated jet is not as diluted as it should be. This could be attributed to the mesh used in their preliminary study.
Ref. [27] built upon their previous study [26] to include the tail of the return point density current analysis with the LES model in OpenFOAM. The standard k-ε RANS turbulence model was run to compare the results to the LES model, and the 45 • and 60 • inclinations were selected. The jet spreading layer thickness was defined as the lateral distance where the concentration would reach 5% of the maximum concentration at the jet centerline. The 60 • jet spreading layer results were compared to those from experiments by [37]. Similar to the study of [26], the eddy structures and concentration profiles were plotted along the jet trajectory, and the results showed that LES was well able to predict the return point location and the coordinates of the centerline peak. However, their model underestimated the dilution at the return point by approximately 20% compared to data collected by previous experiments. This could be due to the mesh resolution that was adopted in this study, as they argued. Comparing to the experimental data, LES was able to reproduce the concentration density in the density current region, while k-ε results could not do that. It was observed that for both geometrical and flow/dilution characteristics, k-ε results were generally lower in values compared to the LES results.
Ref. [46] built on the work of [26,27] by studying the effect of swirls at the outfall nozzles. Experimentally, they found that the addition of swirls at the nozzles lowered the terminal rise height of inclined dense discharges and increased the dilution rate at the return point. This could potentially affect the outfall design for coastal waters. Their numerical study is the first of its kind for the CFD modeling of swirls in effluent discharges. Similar to [27], two approaches were taken in consideration for this study: (i) LES using the dynamic Smagorinsky sub-grid approach and (ii) RANS using the standard k-ε turbulence model. In the case with swirls at the nozzle, RANS had a better performance than LES when the results were compared to the experimental data, which could be due to the reduction in turbulence anisotropy and more axisymmetric distribution of turbulence. On the other hand, in the case of no swirls at the nozzle, LES performance was superior, as it could account for the anisotropy in a more accurate manner. When the swirls number (maximum tangential velocity divided by the discharge velocity at the nozzle) was larger than 0.33, it was observed that both RANS and LES models underestimated the dilution of discharge.

Discussion on Differences in RANS and LES Models for Effluent Mixing Problems
RANS and LES models are very popular models in the modeling of effluent discharges in ambient waters, as reviewed in the studies above. Navier-Stokes equations were described before in Equations (8) to (11). Reynolds decomposition involves splitting any instantaneous quantity in mean and fluctuating components by time averaging for both velocity and momentum equations. RANS equations resemble basic governing equations, except for the turbulent momentum and density fluxes, which results from the Reynolds decomposition and averaging. Six unknown terms of the Reynolds stress tensor and three unknown terms of the density flux term imply that the number of unknowns is larger than the available equations. Therefore, it leads to an undetermined system of equations commonly referred to as the closure problem. To resolve this severe shortcoming, several hypotheses and methods are prescribed. The turbulent viscosity and gradient-diffusion hypotheses are the most widely used concepts to deal with the closure problem. The turbulent viscosity hypothesis (TVH) assumes that the deviatoric Reynolds stress is proportional to the mean shear strain rate as [47]: where k = 1 2 u i 2 = 0.5 u 2 + v 2 + w 2 is the turbulent kinetic energy and υ t is the turbulent eddy viscosity. The gradient-diffusion hypothesis (GDH) assumes that the turbulent density (scalar) flux is aligned with the mean density (scalar) gradient as [47]: where κ t is a positive scalar that is named the turbulent eddy diffusivity. These hypotheses resolve the closure problem by decreasing the number of unknowns, but they still require a correct proposition for the turbulent viscosity and diffusivity. Different closure schemes are introduced widely to define the turbulent viscosity (υ t ). Regarding how to solve for υ t , in terms of additional equations, these closure schemes are classified as zero-equation, one-equation or two-equation models. Zero-equation or algebraic models do not require additional partial differential equations (PDEs) for transport equations and provide a prediction of the turbulent viscosity (υ t ) directly from the mean flow variables. One-equation models involve the use of one additional transport equation (usually turbulent kinetic energy) and assess the turbulent viscosity (υ t ) based on the estimated turbulent quantity. Two-equation RANS closure schemes such as the standard k-ε model make use of two more transport equations for turbulence quantities to define the turbulent viscosity (υ t ). However, to provide closure for the turbulent flux term in the density transport equation, most turbulence schemes make use of a turbulent Prandtl number (Pr t ) instead of defining the turbulent diffusivity (κ t ) explicitly. The turbulent Prandtl number is defined as: Due to averaging, the accuracy of RANS models is always a concern in mixing problems, and the development of more robust numerical schemes is an ongoing effort. On the other hand, LES models are now getting more attention, especially with the advancement in computational resources. LES models would resolve the unsteady nature of large eddies and use turbulence models for small-scale eddies. To handle this, LES uses a filtering technique to split the velocity field into mean (U) and residual (u') values. The filter size can either be determined implicitly by the numerical domain grid size or by introducing filter functions [48]. The LES equations for a 3D unsteady-state flow using Boussinesq approximation can be written as follows: The filtered density transport equation is: In Equations (25) and (26), τ SGS ij and χ SGS j are the SGS tensor and sub-grid scale flux vector respectively and are defined as: With the existence of residuals, similar to RANS models, LES suffers from the closure problem and requires SGS models to resolve closure. The accuracy of LES is related to the efficiency of SGS models used to define the SGS motions. The simplest and most well-known closure method is introduced by [49]. A linear turbulent viscosity (υ t ) is used to model the residual motions as: The computational studies on the mixing of effluent (review of previous studies in this paper) in the near-field region showed that, especially for LES models, the influence of small changes in the discharge momentum (i.e., inflow momentum) can result in large variations in the concentration field [50]. As seen in previous studies, for effluent mixing problems, RANS models perform relatively well in replicating the geometrical and flow properties of effluent discharges and thus need more exploration in terms of the sensitivity of these models to various input parameters. RANS models are very popular in mixing applications between researchers and professionals, as it has a low CPU demand and is still an appropriate methodology. There are some applications such as the mixing of toxic agents that the models should go beyond the estimate of mean concentration and should resolve smaller scales of the fluid motions. In these cases, LES models are a suitable approach, and they have been an approved method based on studies in the past few years.

Vertical Jets
Vertical discharges have traditionally been preferred over inclined jets when there are horizontal confinements, and they are usually more suitable for deep waters where there is enough height in the water column to diffuse the concentration. This section reviews the previous CFD studies completed on vertical jets. Unlike inclined dense discharges, this section includes both positively and negatively buoyant jets. Vertical jet modeling is attractive for researchers due to its simplicity and practicality. Ref.
[51] modeled a confined turbulent buoyant discharge using [52]'s realizable k-ε model and modified it to account for turbulence due to buoyancy. They also validated the numerical model with an experimental study, which included the jet temperature, axial velocities, radial velocity, dynamic pressure, centerline velocity, mass, and momentum fluxes. This study is limited because it only performs 2D modeling of the 3D experiment. Numerically, they compared the realizable k-ε model with other k-ε models (standard k-ε, RNG k-ε, and k-ω) for temperature propagation in the jet. The realizable k-ε model performed the best among all models, especially closer to the nozzle, and it satisfies some constraints in the Reynolds' stresses, which make it more suitable for this specific physics of flow. Ref. [24] made a similar recommendation but for inclined dense discharges. They expressed an essential coefficient of the realizable k-ε model, C µ , as a function of mean flow and turbulence properties, instead of assuming that it was a constant as in the standard model. This expression (i.e., realizability) satisfies certain constrains on the Reynolds stresses, which makes this model stronger compared to the other k-ε models for mixing problems.
Ref. [53] also studied vertical dense discharges, with both experimental and numerical investigations. A single round outfall was modeled experimentally and numerically. Several numerical simulations were completed with the ANSYS Fluent CFD package, and the numerical results were compared to the experimental data. They considered a homogeneous ambient water with no current and reported the geometrical and concentration properties of the discharge along its trajectory. A wide range of conditions was covered by studying various combinations of port diameters and concentrations of effluent salinities. Their experimental results were compared to the experimental data in the literature. However, this study only compared the numerical results with previous semi-empirical formulas, so it lacks comparison to a good experimental dataset, even to their own data. The numerical model also identified a penetration depth and the existence of multiple peaks for the brine concentration.
Ref. [54] investigated laterally confined vertical buoyant jets with the OpenFOAM model. They modified the standard k-ε turbulence model with the GGDH buoyancy approach and linked the Prandtl and turbulent Prandtl numbers to the Froude number. They verified the model using experimental results from [55] and extracted normalized concentration Gaussian profiles along the jet centerline, perpendicular to the trajectory. Ref. [56] calculated the free jet concentration and comparing this to the modeled confined jet concentration in this study shows that using the same Froude number gives a higher dilution at a certain point along the jet centerline in confined jets. The study had the following results: • An assumption was made that improves the determinacy and aids in the calibration of CFD modeling, which is: the Prandtl number Pr and turbulent Prandtl number Pt r would be related to the densimetric Froude number Fr d : • The influence of water surface is smaller than the influence of confinement, although the distribution of concentration along a cross-section in a confined discharge could be impacted by boundaries. • A conclusion was made to enable engineers and researchers to quickly estimate the evolution of a laterally confined vertical buoyant discharge, which is that the rate of discharge concentration growth is almost equal to b gc /s = 0.0938, where b gc is the concentration 1/e width (e is the Napier's constant) and that is where the jet concentration reaches to 1/e of the centerline maximum concentration.
Ref. [57] used the OpenFOAM model to simulate the mixing properties of vertical buoyant jets discharged from multiport diffusers. Four turbulence models were investigated: standard k-ε, RNG k-ε, k-ω, and SST k-ω. By comparing cross-sectional variations of mean axial velocities at different elevations, they observed that the RNG k-ε model performed the best of the four models. By varying the ratio of the port spacing over the port diameter, they demonstrated that the port spacing affects the dilution characteristics of the jets (i.e., the concentration decreases with increasing port spacing). This study found that when the merging point is located above the nozzle, the centerline concentration drops to approximately 20% of the initial concentration, and the port spacing effect becomes less important. Therefore, a logarithmic term was added to derive a new empirical formula, which includes the port spacing effect to describe the centerline dilution: where d p is the port diameter, p s is the port spacing, and α, β, and γ are regression coefficients.

Horizontal Jets
There are two categories of horizontal jets: offset jets and wall jets. Wall jets are attached to a solid boundary (usually at the bottom), and offset jets are elevated away from the bottom and experience no effects from the wall. This section covers numerical studies on both types to date.
Ref. [58] simulated horizontal buoyant wall jets in three dimensions by applying a realizable k-ε model to analyze results for temperature dilution, centerline trajectory, and cling length at various sections, as well as jet velocity. They compared their results to those of [59,60]. A linear relationship was found for the cling length as a function of the densimetric Froude number: L/D=3.2Fr d . An exponential relationship was found for the jet centerline velocity.
The conclusions from [58] included the following. (i) Temperature dilution is related to the nozzle diameter D, the distance of the jet trajectory to the nozzle x, and the densimetric Froude number Fr d . The temperature dilution in the wall jet region decays according to Equation (32). (ii) The central surface velocity profiles resemble that of a turbulent wall jet and generally agree with the classical curve of wall jets. (iii) The velocity profile shows strong similarity in the vertical plane beyond an x > 5D distance from the nozzle, fitting to a Gaussian shape. (iv) The centerline velocity decay is related to the densimetric Froude number Fr d and the nozzle diameter D and fits Equation (33), where U 0 is the velocity at the source and U m0 is the centerline maximum velocity.
Ref. [61] used the LES model to study, for the first time, the interaction between a parallel offset jet and a plane wall jet; the large eddies were simulated directly, and the small eddies were obtained with the Dynamic Kinetic energy Sub-grid-scale Model (DKSM) and the Dynamic Smagorinsky-Lily Model (DSLM). In predicting the turbulent intensity and the mean stream-wise velocity, the agreement was good between the numerical results and the experimental data, especially for DKSM.
Ref. [61] also analyzed some aspects of the turbulence mechanism, such as the function of correlation, the velocity Probability Density Function (PDF), and the coherent structures. In some locations, the mean stream-wise velocity profiles were similar a certain distance after merging two jets. Near the jet exit, the velocities along the centerline of the two jets were negative. The stream-wise velocities and turbulent intensity both became positive and rapidly increased to the maximum after a stagnation point, and then they gradually decreased downstream; this shows that the interaction of jets is dominant near the exit. Finally, tracer characteristics were used to demonstrate the dilution by the dual jet. The constant concentration was maintained in both the stream-wise and wall normal directions close to the exit because of the interaction between the two jets and the presence of the wall. After the two jets merge completely, the concentrations (C/C 0 ) are distributed in the wall normal direction parabolically, while the maximum value decreases along the stream-wise direction linearly. The concentration (C/C m ) profiles were similar in the region with a steady C 0 value.
Horizontal wall jets ( Figure 6) were examined by [23], using a finite-volume model (FVM) and several turbulence models. Their study aimed to discover whether the RANS models were suitable for predicting the velocity and concentration properties of wall jets. They also aimed to identify the best performing turbulence model. The jets were modeled with OpenFOAM, and those turbulence models were compared to those from the numerical studies of [58,59] (Fr d = 11.61 to 42.33). The linear (standard k-ε, RNG k-ε, realizable k-ε and SST k-ω) and Reynolds stress (Launder-Gibson and LRR) turbulence models were tested in their study. This was the first study for wall jet modeling that compared several turbulence models. velocities and turbulent intensity both became positive and rapidly increased to the maximum after a stagnation point, and then they gradually decreased downstream; this shows that the interaction of jets is dominant near the exit. Finally, tracer characteristics were used to demonstrate the dilution by the dual jet. The constant concentration was maintained in both the stream-wise and wall normal directions close to the exit because of the interaction between the two jets and the presence of the wall. After the two jets merge completely, the concentrations (C/C0) are distributed in the wall normal direction parabolically, while the maximum value decreases along the stream-wise direction linearly.
The concentration (C/Cm) profiles were similar in the region with a steady C0 value. Horizontal wall jets ( Figure 6) were examined by [23], using a finite-volume model (FVM) and several turbulence models. Their study aimed to discover whether the RANS models were suitable for predicting the velocity and concentration properties of wall jets. They also aimed to identify the best performing turbulence model. The jets were modeled with OpenFOAM, and those turbulence models were compared to those from the numerical studies of [58,59] (Frd = 11.61 to 42.33). The linear (standard k-ε, RNG k-ε, realizable k-ε and SST k-ω) and Reynolds stress (Launder-Gibson and LRR) turbulence models were tested in their study. This was the first study for wall jet modeling that compared several turbulence models. Figure 6. Schematic view of the model by [23].
As the fluid leaves the nozzle in a wall jet attached to the horizontal wall and discharging effluent into a water body, water entrains the jet from all directions except for the wall region. Therefore, there is a higher pressure exerted on the top of the jet than on the wall, and the jet stays on the wall up to the point where the suction pressure from above reduces, and the buoyancy force exceeds the pressure difference. There are four regions for wall buoyant jets: (i) Initial Jet Region: Figure 6. Schematic view of the model by [23].
As the fluid leaves the nozzle in a wall jet attached to the horizontal wall and discharging effluent into a water body, water entrains the jet from all directions except for the wall region. Therefore, there is a higher pressure exerted on the top of the jet than on the wall, and the jet stays on the wall up to the point where the suction pressure from above reduces, and the buoyancy force exceeds the pressure difference. There are four regions for wall buoyant jets: (i) Initial Jet Region: from the inlet to the point where the velocity profile is equal to the maximum initial velocity, and nearly uniform; (ii) Wall Jet Region I: from the end of the Initial Jet Region to the point where the jet centerline departs from the horizontal and starts rising; (iii) Wall Jet Region II: from the end of Wall Jet Region I to the point where the outer layer of the jet rises up off the floor; and (iv) Free Jet Region: starts after the Wall Jet Region. These regions are shown in Figure 6.
The cling length for thermal wall jets is the distance from the nozzle to the location where the floor temperature reaches (T − T a )/ (T 0 − T a ) = 3% [58]. Figure 7 shows the numerical cling length results from [23] compared to the experimental data and numerical results of [58]. The results agreed with both the other researchers' numerical data and the experimental data. For larger Froude numbers, the cling length values reported by [59] were smaller than those from [58], though they align well with the numerical results obtained by [23]. Table 4 shows the relationships between L/D and Frd for each turbulence model evaluated in [23].
All of the turbulence models estimated a smaller L/D than the experimental data. RNG k-ε has the closest cling length value to the experimental result, and SST k-ω has the smallest.
Predicting jet trajectories is a key factor in the outfall design to predict the path that a jet travels from the exit point until it reaches the surface. This is particularly critical in locations where the receiving water is not deep enough to dilute the effluent completely. Figure 8 shows trajectories from several studies, including that of [23]. It seems that the trajectory results from the family of k-ε turbulence models are much more accurate than the SST k-ω model. The results agreed with both the other researchers' numerical data and the experimental data. For larger Froude numbers, the cling length values reported by [59] were smaller than those from [58], though they align well with the numerical results obtained by [23]. Table 4 shows the relationships between L/D and Fr d for each turbulence model evaluated in [23]. All of the turbulence models estimated a smaller L/D than the experimental data. RNG k-ε has the closest cling length value to the experimental result, and SST k-ω has the smallest.
Predicting jet trajectories is a key factor in the outfall design to predict the path that a jet travels from the exit point until it reaches the surface. This is particularly critical in locations where the receiving water is not deep enough to dilute the effluent completely. Figure 8 shows trajectories from several studies, including that of [23]. It seems that the trajectory results from the family of k-ε turbulence models are much more accurate than the SST k-ω model. the closest cling length value to the experimental result, and SST k-ω has the smallest.
Predicting jet trajectories is a key factor in the outfall design to predict the path that a jet travels from the exit point until it reaches the surface. This is particularly critical in locations where the receiving water is not deep enough to dilute the effluent completely. Figure 8 shows trajectories from several studies, including that of [23]. It seems that the trajectory results from the family of k-ε turbulence models are much more accurate than the SST k-ω model. Ref. [23] obtained detailed results on stream-wise velocity profiles ( Figure 9). The stream-wise (x-y) velocity profiles of the buoyant wall jet centerline were extracted from various simulations. The velocity field results were obtained for different jet cross-sections along the x-direction (various values of x/D) at the symmetry plane. In Figure 9, the variables are as follows: U m is the x-direction velocity (along y at the central plane), U m0 is the maximum of U m , with the ordinate of y, and y m/2 is the velocity-half-height, which is the height where U m =U m0 /2. All profiles for stream-wise velocity are self-similar and in good agreement with the data by [63], as plotted in Figure 9. Ref. [60]'s equation, which is suitable for 2D wall jets, is also plotted along with other data for comparison. Ref. [ Profiles of various values of x/D showed a good similarity as seen in Figure 9. Buoyancy-induced instabilities justify the larger deviations in the farther x/D locations as they mostly happen at higher elevations (i.e., y/y m/2 > 1) where the momentum forces dissipate, and buoyancy forces get stronger. The literature often reports the velocity self-similarity profiles at the central plane for both experimental and numerical studies, although usually without presenting results for the offset measurement from the centerline. the velocity-half-height, which is the height where Um=Um0/2. All profiles for stream-wise velocity are self-similar and in good agreement with the data by [63], as plotted in Figure 9. Ref. [60]'s equation, which is suitable for 2D wall jets, is also plotted along with other data for comparison. Ref. [60]'s equation reads: . Self-similarity of stream-wise velocity profiles for various turbulence models (Source: [23]).
Profiles of various values of x/D showed a good similarity as seen in Figure 9. Buoyancy-induced instabilities justify the larger deviations in the farther x/D locations as they mostly happen at higher elevations (i.e., y/ym/2 > 1) where the momentum forces dissipate, and buoyancy forces get stronger. The literature often reports the velocity self-similarity profiles at the central plane for both experimental and numerical studies, although usually without presenting results for the offset measurement from the centerline.
Experimental data for offset velocity profiles were first reported by [63] for two offset sections, z/D = 1.818 and z/D = 3.636. The current study compared numerical results with [63] as well as [60]'s equation, as shown in Figure 10, in which ym/2 is the local length scale and Ums is the maximum velocity for the offset sections. As seen in the figure, the numerical results (current study) for z/D = 3.636 do not agree well with Verhoff's curve in the area close to the nozzle (x/D = 5 and x/D = 10). This is Experimental data for offset velocity profiles were first reported by [63] for two offset sections, z/D = 1.818 and z/D = 3.636. The current study compared numerical results with [63] as well as [60]'s equation, as shown in Figure 10, in which y m/2 is the local length scale and U ms is the maximum velocity for the offset sections. As seen in the figure, the numerical results (current study) for z/D = 3.636 do not agree well with Verhoff's curve in the area close to the nozzle (x/D = 5 and x/D = 10). This is primarily due to the development of the jet in the tank width ( Figure 11). As z/D increases, the jet may not develop at the values along the width of the tank yet, and thus the scatters show less self-similarity. z/D = 1.818 and z/D = 3.636. The current study compared numerical results with [63] as well as [60]'s equation, as shown in Figure 10, in which ym/2 is the local length scale and Ums is the maximum velocity for the offset sections. As seen in the figure, the numerical results (current study) for z/D = 3.636 do not agree well with Verhoff's curve in the area close to the nozzle (x/D = 5 and x/D = 10). This is primarily due to the development of the jet in the tank width ( Figure 11). As z/D increases, the jet may not develop at the values along the width of the tank yet, and thus the scatters show less selfsimilarity.  The current study investigated the temperature results from [23] further. The dilution contours for temperature at the plane of symmetry (z = 0) are plotted for realizable k-ε in Figure 12 for dilution rates of 12, 15, 20, 30, and 60. The innermost contour line reads a dilution of S = 12, and the outermost contour line corresponds to dilution of S = 60. Dilution would obviously increase with distance from the nozzle, and it depends on both discharge and receiving water properties (such as discharge diameter D, densimetric Froude number Frd, and ambient water depth Ha, etc.). It was observed that the effect of distance (the path that the jet travels through) on the dilution factor was larger than the effect of the Froude number at the nozzle.
Temperature profiles (self-similar profiles) would behave the same (i.e., Gaussian shape) for various Froude numbers. Figure 13 presents the temperature Gaussian distribution at different crosssections for three different cases (i.e., three different Froude numbers), while bulked together. As shown, the temperature profiles, similar to the velocity profiles, seem to be independent of the Froude number. The current study investigated the temperature results from [23] further. The dilution contours for temperature at the plane of symmetry (z = 0) are plotted for realizable k-ε in Figure 12 for dilution rates of 12, 15, 20, 30, and 60. The innermost contour line reads a dilution of S = 12, and the outermost contour line corresponds to dilution of S = 60. Dilution would obviously increase with distance from the nozzle, and it depends on both discharge and receiving water properties (such as discharge diameter D, densimetric Froude number F rd , and ambient water depth H a , etc.). It was observed that the effect of distance (the path that the jet travels through) on the dilution factor was larger than the effect of the Froude number at the nozzle.
Temperature profiles (self-similar profiles) would behave the same (i.e., Gaussian shape) for various Froude numbers. Figure 13 presents the temperature Gaussian distribution at different cross-sections for three different cases (i.e., three different Froude numbers), while bulked together. As shown, the temperature profiles, similar to the velocity profiles, seem to be independent of the Froude number. the effect of distance (the path that the jet travels through) on the dilution factor was larger than the effect of the Froude number at the nozzle.
Temperature profiles (self-similar profiles) would behave the same (i.e., Gaussian shape) for various Froude numbers. Figure 13 presents the temperature Gaussian distribution at different crosssections for three different cases (i.e., three different Froude numbers), while bulked together. As shown, the temperature profiles, similar to the velocity profiles, seem to be independent of the Froude number.  Ref. [64] used a sigma SGS eddy viscosity model to perform an LES model of turbulent horizontal buoyant and non-buoyant jets with several Reynolds (Re) and Richardson (Ri) numbers in stratified ambient water. Then, the effects of varying Ri (to characterize density differences) and Re (to characterize injection momentum) was studied. Turbulent production and dissipation rates were found to be asymmetric in the mid-vertical plane but symmetric in the horizontal plane. Ref. [64]'s results on jet centerline velocity decay, mean velocity self-similarity, radial spread, and turbulent fluctuations were in good agreement with previous experimental results. Using the results from their LES model and looking at the instantaneous velocities in the jet studies, they realized that the jet close to the nozzle would behave the same as the developed jet farther from the discharge point in the stratified ambient water. They also identified the stable and unstable stratification regions and how the turbulent vortex rings and Ri are related to them. Circular, square, and rectangular nozzles were considered in [65]'s CFD study of the effects of nozzle geometry on turbulent offset jet development in the near-field region. Turbulence models that were studied include the standard k-ε model [66], the realizable k-ε model [67], the Launder-Sharma k-ε model [68], and the Yang-Shih k-ε model [69]. The Yang-Shih k-ε turbulence model came out on top for predicting the jet properties in a comparison of these turbulence models (with a Re number of approximately 8500) with previous experimental studies. Square-shaped offset jets seem to spread more in the wall normal and lateral directions and result in more efficient mixing with the Ref. [64] used a sigma SGS eddy viscosity model to perform an LES model of turbulent horizontal buoyant and non-buoyant jets with several Reynolds (Re) and Richardson (Ri) numbers in stratified ambient water. Then, the effects of varying Ri (to characterize density differences) and Re (to characterize injection momentum) was studied. Turbulent production and dissipation rates were found to be asymmetric in the mid-vertical plane but symmetric in the horizontal plane. Ref. [64]'s results on jet centerline velocity decay, mean velocity self-similarity, radial spread, and turbulent fluctuations were in good agreement with previous experimental results. Using the results from their LES model and looking at the instantaneous velocities in the jet studies, they realized that the jet close to the nozzle would behave the same as the developed jet farther from the discharge point in the stratified ambient water. They also identified the stable and unstable stratification regions and how the turbulent vortex rings and Ri are related to them. Circular, square, and rectangular nozzles were considered in [65]'s CFD study of the effects of nozzle geometry on turbulent offset jet development in the near-field region. Turbulence models that were studied include the standard k-ε model [66], the realizable k-ε model [67], the Launder-Sharma k-ε model [68], and the Yang-Shih k-ε model [69]. The Yang-Shih k-ε turbulence model came out on top for predicting the jet properties in a comparison of these turbulence models (with a Re number of approximately 8500) with previous experimental studies. Square-shaped offset jets seem to spread more in the wall normal and lateral directions and result in more efficient mixing with the surrounding fluid than circular offset jets. However, the maximum shear stress on the adjacent wall in the case of the square-shaped nozzle was slightly higher than that in the circular nozzle case.
Ref. [70] studied turbulent circular wall jets both experimentally and by the LES model; then, they compared those results with the numerical results from two RANS models: the standard k-ε and standard k-ω models, with enhanced wall functions. Then, [70] tested the grid convergence by the Grid Convergence Index (GCI), with low-resolution (0.9 million) followed by high-resolution (3.1 million) cells and ending with both the velocity and concentration properties of the jets. Ref. [70]'s results showed that LES is better than both RANS models for reproducing the scalar mixing and kinematic characteristics. Their LES results generally agreed better with the experimental data, although they did underpredict the span-wise velocity profiles away from the jet centerline. Vorticity distribution and turbulence intensities (u'/U, with u' being the root-mean-square of turbulent velocity fluctuations and U being the mean velocity) were also extracted and compared to the past experimental data, which showed a better agreement of the LES model rather than RANS models. In the region far from the centerline, the span-wise turbulence intensity was observed to decrease faster, which was likely due to the inadequacy of the Smagorinsky SGS model for low-turbulence intensity situations in those regions. In the plane of symmetry, the y-direction (lateral-perpendicular to the plane of symmetry) vorticity, caused by wall presence, dominated the vorticity distribution. In the horizontal plane, the primary contributor was observed to be the z-direction (vertical-perpendicular to the horizontal plane) vorticity due to the jet-flow shear layer in the ambient water.
In the most recent study, [71] studied the offset buoyant jets with various properties for the discharges and ambient water, experimentally and numerically. Discharges were set to be both thermal and non-thermal, but positively buoyant all the time. The ambient water was stagnant, and they used a PIV system to collect the experimental data. All comparative experiments were conducted with the same densimetric Froude numbers (F rd , ranging from 9.9 to 29.8) and density differences (∆ρ, ranging from 5.1 to 17.41). Three RANS turbulence models were adopted for their numerical study: standard k-ε, realizable k-ε, and buoyancy-modified k-ε. They concluded that the realizable k-ε model was more successful in predicting the discharge trajectories. The main finding of this study was that while using different combinations of parameters in discharge (salinity versus temperature) for keeping the same properties of the jet (the same values of F rd and same ∆ρ), the trajectory and mixing characteristics of the jets would be different in the same ambient water. Therefore, it is important not to only look at the relative buoyancy between discharge and ambient water, but also the properties of discharge such as salinity and temperature, which could be very important in the overall mixing efficiency of the jets.

Surface Discharges
Surface discharges of effluent into water bodies are less common due to lower mixing efficiency with the ambient water (i.e., the top portion of the jet/plume does not get ambient water entrainment especially close to the nozzle where the momentum effect is higher). This could be the reason that surface discharges have been less studied both experimentally and numerically. There is no CFD modeling of surface discharges to the best knowledge of the authors of this paper. The CORMIX3 empirical-based model was developed for surface discharge outfalls, which is not considered a CFD tool. This gap in the literature could be bridged with the following suggestions for both experimental and numerical modeling of surface discharges: • Surface discharges using different channel geometries in calm ambient water • Surface discharges using different channel geometries in co-flow ambient water • Surface discharges using different channel geometries in cross-flow ambient water

Discharge Port Configuration
All of the cited publications above focused on the numerical modeling of single port outfalls. However, multiport diffusers are commonly used in ocean outfall designs due to their efficiency. The configuration of nozzles could be varied depending on the ambient condition and design considerations. Figure 14 shows five different configurations of multiport diffusers that are commonly used. Another emerging multiport diffuser configuration is the rosette outfall configuration, as shown in Figure 15. Similar to the surface discharges, multiport effluent discharges have been less studied, both experimentally and numerically.
surface discharges have been less studied both experimentally and numerically. There is no CFD modeling of surface discharges to the best knowledge of the authors of this paper. The CORMIX3 empirical-based model was developed for surface discharge outfalls, which is not considered a CFD tool. This gap in the literature could be bridged with the following suggestions for both experimental and numerical modeling of surface discharges: • Surface discharges using different channel geometries in calm ambient water • Surface discharges using different channel geometries in co-flow ambient water • Surface discharges using different channel geometries in cross-flow ambient water

Discharge Port Configuration
All of the cited publications above focused on the numerical modeling of single port outfalls. However, multiport diffusers are commonly used in ocean outfall designs due to their efficiency. The configuration of nozzles could be varied depending on the ambient condition and design considerations. Figure 14 shows five different configurations of multiport diffusers that are commonly used. Another emerging multiport diffuser configuration is the rosette outfall configuration, as shown in Figure 15. Similar to the surface discharges, multiport effluent discharges have been less studied, both experimentally and numerically.  Ref.
[72] studied a submerged multiport diffuser discharging thermal water into ambient water with co-flow conditions (similar configuration as shown in Figure 14a). They developed their own  Ref.
[72] studied a submerged multiport diffuser discharging thermal water into ambient water with co-flow conditions (similar configuration as shown in Figure 14a). They developed their own 3D model and used a numerical slot diffuser concept to implement the momentum discharge flux properly in the model. Thermal jet trajectory and temperature distribution along the trajectory were Ref. [72] studied a submerged multiport diffuser discharging thermal water into ambient water with co-flow conditions (similar configuration as shown in Figure 14a). They developed their own 3D model and used a numerical slot diffuser concept to implement the momentum discharge flux properly in the model. Thermal jet trajectory and temperature distribution along the trajectory were compared between their model and the experimental data. The qualitative comparison of the results showed a relatively good agreement between the numerical results and experimental data.
Very recently, [73] published a paper on 3D numerical simulations of rosette multiport diffusers using two RANS models: standard k-ε and RNG k-ε turbulence models. They used the popular OpenFOAM model for this study. They introduced an initial dilution region (y < 5D), and this is the region in which jets do not interact with each other after discharge and concentration is higher with a lower dilution rate. As the jets grow and travel toward the surface, their lateral widths increase, and they start interacting with each other. It was reported that for the regions around y = 25D, the concentration was reduced compared to the initial dilution region, and thus the dilution rate increased. As the jets traveled farther, their interaction with each other increased (e.g., y = 65D). The RNG k-ε turbulence model performed better than the standard k-ε turbulence model when the numerical results were compared to the experimental data. As most of the terms in both turbulence models are the same (i.e., the same transport equations for k and ε), the difference in the results could be attributed to the improvements made in the RNG k-ε turbulence model as listed below: • The RNG k-ε model accounts for the influence of the Reynolds number on the effective turbulence transport.

•
The RNG k-ε model calculates the inverse effective Prandtl numbers, using a more advanced equation.

•
The RNG k-ε model includes a new term in the ε transport equation that improves the calculation of the turbulent viscosity.
Ref. [73] also appreciated that the computational cost of using the RNG k-ε turbulence model is almost the same as that of standard k-ε turbulence model, while models such as LES and DNS will increase the computational cost dramatically.
Other than the above-mentioned studies, no CFD modeling study was found for the modeling of multiport diffusers. The main reason that the present available CFD models have not been applied to the modeling of multiport diffusers could be attributed to the difficulties in the mesh generation and stability of such models. However, it seems that with the work of [73], there will be more CFD modeling of such cases in the near future.

Critical Review and Future Research Needs
The above presented review and discussion outline the state-of-the-art knowledge on the CFD modeling of effluent discharge in the context of near-field mixing and dispersion. The cited literature above covers fundamental processes of discharge mixing modeling, which would help the efficient design of outfalls for engineering applications.
In advancing the numerical modeling of effluent discharges, special attention should be given to collecting reliable field data for real-sized projects in real project conditions. Coupling near-field and far-field models, which are the most practical in engineering projects, should be studied in more detail to investigate coupling techniques. Most of the CFD studies focused on fine-tuning small-scale parameters in the numerical schemes to calibrate their models, and little attention was spent on the applicability of such models. This is specifically true when it comes to the practice of outfall designs, where wave and wind-induced currents have inevitable effects on jet mixing and dispersion in the near field. Now that the numerical models have shown an excellent performance in replicating the experimental cases, it is time to collect more field data to simulate the large-scale practical cases to further investigate the challenges in complex geometries, environmental forcing, and boundary conditions. Additionally, there are a few jet configurations that have been studied less numerically or have not been studied at all. Although surface discharges are less popular due to poor dilution efficiency, there are several scenarios in which they are more suitable compared to inclined, submerged jets. Examples are available for the shallow places, or when the maintenance of submerged outfalls are more difficult compared to surface channels. Moreover, most of the previous numerical studies have focused on stagnant ambient water in their CFD models. More realistic ambient conditions with the presence of co-flow and cross-flow will advance the applicability of CFD models in effluent discharge problems. In reality, and especially for coastal outfalls, the discharge point is often laid on a steep bed, which will have a significant impact on the jet dispersion and run-out distance, especially beyond the impact point (i.e., when the density current on the bed is formed). Ref. [27] have studied the dilution characteristics of density currents built after the impact point of an inclined dense jet. However, their study was on a flat bed with no inclination. CFD models could be used to efficiently model these practical scenarios in more detail. The effect of bed roughness and vegetation could be critical for turbulent wall jets and/or density currents beyond the impact point.
Multiport and rosette diffusers could potentially increase the efficiency of the effluent discharges by increasing the number of jets discharging at the same time and multiplying the momentum-length scales before merging the jets. CFD modeling of complex multiport and rosette jet configurations would provide valuable insight into designing these jets in practice for professional engineers. Ambient stratification is another aspect that could be investigated using CFD models in more detail for both stagnant and non-stagnant ambient conditions.

Conclusions
Advances in computational resources and strength have led to an increase in CFD modeling of jet mixing and dispersion over the past two decades. Some jet configurations have received more attention than others due to various reasons. These reasons include but are not limited to easier compilation in numerical models, less complicated boundary conditions, simpler geometries for mesh generation, more practical efficiency in practice, and the availability of data. However, recent advances in experimental setups in jet studies (e.g., PIV and LIF systems), sampling techniques in field studies, and the availability of open-source CFD tools have opened the doors to more realistic CFD modeling in jet studies. Based on an extensive literature review of CFD modeling of effluent discharges in the near-field region, the following conclusions can be drawn: • Numerically, the most studied effluent discharge configuration has been in inclined dense jets, due to their applicability in industry. Most studies focused on lab-size experiments to calibrate their models. Details on jet trajectories and dilution and velocity characteristics have been investigated and compared to experimental data. RANS and LES turbulence models are popular for such studies.

•
Vertical jets are also popular in CFD studies, and the new trend in studying these jets involves considering the ambient conditions that may affect these jets such as lateral confinement and water shallowness, where the jet is attached to the top boundary. Cross-flow in vertical jets could have a significant influence in terms of the trajectory and dilution, both of which are getting more attention from researchers using CFD models.

•
Horizontal jets could be either positively buoyant or negatively buoyant with attachment to the bed (i.e., wall jets) or elevated (i.e., offset jets). Single jets have been studied experimentally and numerically during the past years, and more attention is now given to interactions of multiple horizontal jets when they merge after a certain distance from the discharge point.

•
Surface discharges have not been studied yet using a CFD approach, even though they are used in the industry and experimental data are available on such jets. This literature gap also exists for more complex jet configurations such as multiport and rosette diffusers.
• Previous studies have mainly focused on the stagnant ambient water due to the simplicity of internal and boundary conditions. However, to replicate real-life conditions more precisely, there is a need to move toward more complex ambient conditions to study the effects of wave, wind, co-flow and crossflow, density stratification, etc. on jet mixing and dispersion in CFD models. • A wide range of turbulence models are available and have already been implemented in different CFD platforms that could be used for discharge mixing studies. Modifications of turbulence models, such as implementing the buoyancy terms, has been shown to be effective in improving the prediction of jet characteristics. Acknowledgments: The second author, H.K.G., is a recipient of the NSERC Postgraduate Scholarship-Doctoral (PGSD).

Conflicts of Interest:
The authors declare no conflict of interest.