Revisiting Mean Flow and Mixing Properties of Negatively Round Buoyant Jets Using the Escaping Mass Approach (EMA)

The flow formed by the discharge of inclined turbulent negatively round buoyant jets is common in environmental flow phenomena, especially in the case of brine disposal. The prediction of the mean flow and mixing properties of such flows is based on integral models, experimental results and, recently, on numerical modeling. This paper presents the results of mean flow and mixing characteristics using the escaping mass approach (EMA), a Gaussian model that simulates the escaping masses from the main buoyant jet flow. The EMA model was applied for dense discharge at a quiescent ambient of uniform density for initial discharge inclinations from 15◦ to 75◦, with respect to the horizontal plane. The variations of the dimensionless terminal centerline and the external edge’s height, the horizontal location of the centerline terminal height, the horizontal location of centerline and the external edge’s return point as a function of initial inclination angle are estimated via the EMA model, and compared to available experimental data and other integral or numerical models. Additionally, the same procedure was followed for axial dilutions at the centerline terminal height and return point. The performance of EMA is acceptable for research purposes, and the simplicity and speed of calculations makes it competitive for design and environmental assessment studies.


Introduction
The effluents of a wastewater treatment plant or power plant cooling waters are usually discharged in large water bodies like ocean, lakes, etc. Similarly, the gases from cooling towers and chimneys are emitted in the atmosphere. In the case of the latter, the effluents' density is less than the ambient density. Consequently, the flow has positive buoyancy, forming turbulent buoyant plumes that affect the receiver. These flows have been extensively investigated through the years [1][2][3][4][5][6][7][8][9]. In contrast, when the discharge fluid is denser than the ambient, the buoyant plumes initially have negative buoyancy, so they are called negatively buoyant jets. In the atmosphere, the emission of dense gases like chlorine, hydrogen fluoride or liquefied natural gas (LNG) occurs when they accidentally leak to the air environment (bounded or unbounded) [10,11]. These gases are usually toxic and/or flammable, causing serious adverse effects to human health, nonhuman biota and generally disasters. The accurate and fast prediction of the flow via mathematical models is crucial in evacuating affected areas and rescuing victims [12,13] and extremely important for the risk assessment of industries and storage areas that handle such hazardous materials [14][15][16].
inclined dense jets in shallow coastal waters, reporting the existence of three different mixing regimes. Papakonstantis et al. [53,54] experimentally studied inclined turbulent round jets with negative buoyancy discharging in a calm homogeneous fluid, providing geometrical and mixing characteristics for discharge angles from 45 • to 90 • to the horizontal. The same experimental apparatus was used by Nikiforakis et al. [55] to study negatively buoyant jets discharging at 30 • , 45 • and 60 • and, more recently, by Papakonstantis and Tsatsara [56,57]. The effect of nozzle inclination was studied experimentally by Oliver et al. [58,59], Abessi and Roberts [60] and Crowe et al. [61,62]. Oliver et al. [58,59] used the LIF technique, whil Crowe et al. [61,62] used the particle tracking velocimetry (PTV) flow visualization technique to study experimentally dense jets issuing from initial angle ranged from 15 • to 75 • in a homogenous calm ambient. Abessi and Roberts [60] used a three-dimensional LIF (3DLIF) to appraise the major flow characteristics from 15 • to 85 • dense jets in a homogenous still environment.
Integral models are used for the theoretical investigation of inclined buoyant flows. The integral models are based on either on Eulerian or Lagrangian approach [2,8,63]. The closure of integral models is based on either the entrainment concept [64] or the jet spreading rate assumption of the mean axial velocity and concentration fields [9,33]. A major disadvantage of the entrainment assumption is that entrainment does not have a global behavior for positive and negative buoyant jets. To overcome the disadvantage, Kaminski et al. [35] proposed a universal entrainment relationship, while Papanicolaou et al. [65] a reduced entrainment coefficient.
The two most widely used integral models are CorJet [7] and the Lagrangian model JETLAG/VISJET [8,66]. These models were initially developed to predict discharges of positive buoyancy. Jirka [67] applied CorJet to investigate inclined negatively buoyant discharges into stagnant homogenous ambient, and Lai and Lee [51] used JETLAG/VISJET incorporating a reduced entrainment coefficient for dense jets.
The simulation of negatively buoyant discharges appears to be more complicated than that of positively buoyant discharges. Lane-Serff et al. [42], Kikkert et al. [49], Ferrari and Querzoli [47] reported the phenomenon of fluid portions that detached from the main flow of the buoyant jet moving vertically upwards causing gravitational instabilities and consequently destroying the axial symmetry of transverse concentration and velocity profiles [49][50][51][52]54,60,[68][69][70][71][72][73]. This situation is favored by weak or sometimes negative values of transverse velocity in the concave region of negatively buoyant jets, as reported experimentally by Crowe [69] and theoretically by Yannopoulos and Bloutsos [74].
Although integral models are generally accepted as reliable tools to simulate inclined buoyant jets, the last years there is a notable trend of using computational fluid dynamics (CFD) to study the behavior of these flows. Vafeiadou et al. [75] studied inclined negatively buoyant jets using the commercial software ANSYS CFX and the shear stress transport (SST) model for turbulence closure. Comparing their findings to available experimental data, they concluded a slight underestimation of terminal rise height and a considerable underestimation of the return point. Additionally, Oliver et al. [76] used ANSYS CFX to solve Reynolds averaged Navier-Stokes (RANS) equations to investigate the geometrical and mixing characteristics of inclined negatively buoyant jets employing for turbulence closure the standard form of k-ε model and a calibrated k-ε model through adjustment of the turbulent Schmidt number in the tracer transport equation. Gildeh et al. [77] used the open-source CFD code called OpenFOAM to evaluate the accuracy of seven turbulence models in simulating wall thermal or saline buoyant jets, concluding that realizable k-ε and Launder-Reece-Rodi (LRR) closure turbulence models were found to be the most accurate and capable of accurately modeling thermal and saline wall jets discharged into stationary ambient water. Gildeh et al. [78,79]) studied the behavior of 30 • and 45 • inclined dense jets in stationary ambient using OpenFOAM. They applied several closure-turbulence-models to evaluate the accuracy of CFD predictions of geometrical and mixing characteristics. They reported that differences occur among the models with realizable k-ε and LRR tending to be more accurate among them. The realizable k-ε was used by Ardalan and Vafaei [80] too, to numerically model thermal-saline effluent discharging at an angle of 45 • in a homogenous and still environment. More recently, large eddy simulations (LES) were used by Zhang et al. [81,82], Cintolesi et al. [83] and Jiang et al. [72] to Fluids 2020, 5, 131 4 of 21 study in detail the geometrical mixing and turbulent characteristics of inclined dense jets, reporting that there are still challenges to simulate accurately the mixing behavior of inclined dense jets.
Briefly, the flow field of an inclined buoyant jet discharged into a stationary ambient fluid of uniform density is described by the initial flux of volume µ 0 = A 0 w 0 , the initial flux of momentum m 0 = A 0 w 0 2 and the initial flux of buoyancy β 0 = g 0 A 0 w 0 c 0 , where w 0 is the initial discharge velocity, A 0 = πd 0 2 /4 is the area of the exit nozzle, d 0 is the diameter of the exit nozzle, c 0 is the initial concentration of tracer, g 0 = g(ρ a − ρ 0 )/ρ 0 is the apparent gravitational acceleration, g is the gravitational acceleration, ρ 0 and ρ a are the densities of the fluid at the jet exit and the ambient, respectively. Yannopoulos and Bloutsos [74], following Shao and Law [50], stated that negative buoyancy occurs when the vertical component of the initial momentum of the buoyant jet is zero, or the jet discharges with purely horizontal momentum and ρ 0 > ρ a so g 0 < 0 or alternatively when the vertical component of the initial momentum is non-zero, and the resulting buoyancy force acts in the opposite direction as the initial motion of the buoyant jet. Additionally, dimensional analysis shows that the behavior of such a flow field is described by three independent non-dimensional variables: the initial densimetric Froude number F 0 = w 0 g 0 d 0 , the initial inclination angle θ 0 measured with respect to the horizontal plane and the dimensionless axial distance Ξ = ξ/l m = (ξ/d 0 ) F 0 −1 , where ξ is the axial geometrical distance and l m is the characteristic length scale l m = d 0 F 0 that is usually employed to make geometrical distances non-dimensional [74]. All other parameters are either constants or variables depending upon independent ones. The constant parameters have the same values both for vertical or inclined buoyant jets of either positive or negative buoyancy. The only extra parameter included in inclined buoyant jets with negative buoyancy is the parameter Λ, which is defined in Chapter 2. More details about the description of the flow field of an inclined buoyant jet discharged into a stationary ambient fluid of uniform density can be found at the relevant papers of Yannopoulos [9] and Yannopoulos and Bloutsos [74]. As the present paper revisits the previously published work regarding the negatively round buoyant jet using the EMA model [74], the authors effort is to incorporate in the current paper an overall review of all the experimental data and numerical results of such problems. The present paper determines the geometric and mixing characteristics of an inclined round negatively buoyant jet using the innovative mathematical model, called escaping mass approach (EMA) [74], which incorporates the exact description of the flow field through a curvilinear coordinate system, the application of a Gaussian integral model including second order terms [9]) and the consideration of escaping masses from the concave region of the inclined buoyant jet [74]. As the curvilinear coordinate system includes the main flow axis ξ and the transverse axis r, it enables the exact integration of the transverse Gaussian profiles of the axial velocity and concentration. Therefore, the calculation of the volume, momentum and buoyancy fluxes are computed exactly along the curvilinear path of the main flow of the buoyant jet.
The EMA model was applied for negatively buoyant jets of initial inclinations in the range 15 • ≤ θ 0 ≤ 75 • and initial Froude number range 4 < F 0 < 100 in a stationary and homogenous ambient environment. EMA's performance was checked by comparing geometric (centerline and external edge terminal height, horizontal location of centerline terminal height and return point, centerline length to terminal height and centerline and external edge length to return point) and mixing (minimum dilution at centerline terminal height and return point) characteristics to available experimental data, numerical and integral model predictions. The good agreement of the results of the EMA model with experimental data obtained from the international literature confirms the reliability of the model. In addition, the better performance in comparison to other integral or numerical model results makes the EMA model appropriate in the design of systems for the disposal of liquid, thermal and gas waste under these conditions. It could also be used for the verification of complex numerical models, but also for laboratory studies of related models.
An inclined negatively turbulent buoyant jet of density ρ 0 is discharged upwards with initial velocity w 0 from a round nozzle of size d 0 with an initial inclination θ 0 to the horizontal plane into a quiescent ambient fluid of uniform density ρ a (ρ 0 > ρ a ). At the first stage, the buoyant jet moves upwards, due to its initial momentum and the mixing with ambient fluid starts. Meanwhile, the negative buoyancy decreases the momentum flux causing its deceleration. The jet rises to a terminal (maximum) height, where the vertical component of momentum becomes zero, and then it falls as a positively buoyant jet to the initial discharge elevation. The general configuration of such a flow is shown in Figure 1, related to a Cartesian coordinate system (O, x, y, z). The flow is assumed to be steady state. The flow centerline (trajectory) is defined as the location of maximum velocity or concentration at each cross-section of the flow. In Figure 1, y and z are the horizontal and vertical distances and subscripts c, e, t and r denote the jet centerline, the external boundary, the terminal rise height and the return point to the source elevation, respectively. Following these definitions z ct and z et are the corresponding terminal heights of the centerline and the external boundary of the jet, which occur at a horizontal distance y ct from the source. The horizontal location of the return point of the jet axis is y cr and the corresponding location for the external boundary is y er . Additionally, S ct and S cr are the axial dilutions at terminal rise height and the return point, correspondingly. Previous studies [32,38,40,45,49,50,53,54,[56][57][58][59][60][61][62]81] have verified that the geometrical quantities non-dimensionalized by d 0 , and dilutions are proportional to F 0 for a specific angle θ 0 . negative buoyancy decreases the momentum flux causing its deceleration. The jet rises to a terminal (maximum) height, where the vertical component of momentum becomes zero, and then it falls as a positively buoyant jet to the initial discharge elevation. The general configuration of such a flow is shown in Figure 1, related to a Cartesian coordinate system (O, x, y, z). The flow is assumed to be steady state. The flow centerline (trajectory) is defined as the location of maximum velocity or concentration at each cross-section of the flow. In Figure 1, y and z are the horizontal and vertical distances and subscripts c, e, t and r denote the jet centerline, the external boundary, the terminal rise height and the return point to the source elevation, respectively. Following these definitions zct and zet are the corresponding terminal heights of the centerline and the external boundary of the jet, which occur at a horizontal distance yct from the source. The horizontal location of the return point of the jet axis is ycr and the corresponding location for the external boundary is yer. Additionally, Sct and Scr are the axial dilutions at terminal rise height and the return point, correspondingly. Previous studies [32,38,40,45,49,50,53,54,[56][57][58][59][60][61][62]81] have verified that the geometrical quantities non-dimensionalized by d0, and dilutions are proportional to F0 for a specific angle θ0.

Materials and Methods
The escaping mass approach (EMA) model [74] is an integral model that describes the flow and mixing fields of inclined round or plane buoyant jets issuing into stationary ambient environment of uniform density. For the case of round buoyant jets, the conservation partial differential equations (PDE) of continuity, momentum and tracer have been formulated in a curvilinear cylindrical coordinate system, which is relative to the corresponding curvilinear orthogonal Cartesian coordinate system for plane buoyant jets [84]. EMA simulates the reported phenomenon ( [42,47,[49][50][51][52]54,60,[68][69][70][71][72][73] of fluid portions that detach from the main flow of the buoyant jet, especially from the inner boundary, and move vertically upwards (or downwards for brines), causing gravitational instabilities, and consequently destroying the axial symmetry of transverse concentration and velocity profiles. EMA calculates the vertical velocity wB of the escaping masses from the buoyant jet field, employing the general definition of the local Richardson number for plume flows [9]. The concentration of these masses cB is estimated as a portion Λ of its centerline value cm. The value of Λ coefficient was adopted equal to 0.34 [74], where Λ includes the value π/4 for round buoyant jets) for all cases in the inclination range −75° ≤ θ0 ≤ −15° (or 15° ≤ θ0 ≤ 75° of negatively buoyant jets). Additionally, EMA is equipped with the second order approach (SOA), to include the variable second-order effect of turbulence to the mean flow properties [9].

Materials and Methods
The escaping mass approach (EMA) model [74] is an integral model that describes the flow and mixing fields of inclined round or plane buoyant jets issuing into stationary ambient environment of uniform density. For the case of round buoyant jets, the conservation partial differential equations (PDE) of continuity, momentum and tracer have been formulated in a curvilinear cylindrical coordinate system, which is relative to the corresponding curvilinear orthogonal Cartesian coordinate system for plane buoyant jets [84]. EMA simulates the reported phenomenon ( [42,47,[49][50][51][52]54,60,[68][69][70][71][72][73] of fluid portions that detach from the main flow of the buoyant jet, especially from the inner boundary, and move vertically upwards (or downwards for brines), causing gravitational instabilities, and consequently destroying the axial symmetry of transverse concentration and velocity profiles. EMA calculates the vertical velocity w B of the escaping masses from the buoyant jet field, employing the general definition of the local Richardson number for plume flows [9]. The concentration of these masses c B is estimated as a portion Λ of its centerline value c m . The value of Λ coefficient was adopted equal to Fluids 2020, 5, 131 6 of 21 0.34 [74], where Λ includes the value π/4 for round buoyant jets) for all cases in the inclination range −75 • ≤ θ 0 ≤ −15 • (or 15 • ≤ θ 0 ≤ 75 • of negatively buoyant jets). Additionally, EMA is equipped with the second order approach (SOA), to include the variable second-order effect of turbulence to the mean flow properties [9].
EMA's development [74] is based on PDEs of volume momentum and tracer's conservation, concerning the Reynolds approximation for mean flow and mixing properties and their fluctuations in steady-state conditions. Additionally, the usually accepted assumptions of (a) Boussinesq's approximation made for small initial density differences; (b) Prandtl's boundary-layer approximation; (c) negligible molecular viscosity terms and (d) no swirl are adopted. Regarding flow symmetry with respect to the centerline ξ axis of the curvilinear coordinate system O l (r, ϕ, ξ), the governing equations of the model are: Continuity where w, u are the mean velocity components in the directions ξ and r, respectively and w , u are their corresponding turbulent fluctuations; r = (n 2 + ψ 2 ) 1/2 is the transverse (radial) distance; h = 1 + r sin ϕ dθ/dξ is the scale factor of the coordinate system, θ is the local inclination angle of the ξ axis, c = (ρ a − ρ)/(ρ a − ρ 0 ) is the local mean concentration and c the corresponding turbulent fluctuation; ρ is the local mean fluid density; w 2 , w c , u c are the local fluxes due to turbulent fluctuations of w, u, c; τ rξ is the mean turbulent shear stress; p d is the dynamic pressure; S 1 = ∂(hq r )/∂r is a source term for the escaping masses from the flow field of the buoyant jet; and q r is the escaping load of masses per unit area in the r direction. The flow field is further described when Equations (1)-(4) satisfy the following boundary conditions: On the jet axis (r = 0) On the jet boundary (r = B w ) where U B = −u B − w B cosθ sinϕ. Based on experimental measurements, the region between the nominal half-width of a buoyant jet, b 0.5 , and the total half-width, B w , is characterized by intermittencies of the velocity field (see Yannopoulos & Bloutsos 2012 [74]). Therefore, our assumption is that all this region could contribute to escaping masses. The velocity w B becomes zero at r = b 0.5 , while has a positive value at the boundary of the buoyant jet (r = B w ). The escaping masses from the convex boundary enter the buoyant jet and are thus swept by the main flow. However, only masses from the concave boundary manage to escape when the entrainment of the buoyant jet becomes very weak or negative. Consequently, in addition to Equation (6), on the jet boundary (r = B w ) further boundary conditions are applied: where w B and c B are the vertical velocity and concentration of escaping masses; u B is the entrainment velocity of the inclined buoyant jet at the actual boundary B w = d 0 /2 + n w b w and n w = 2 −1/2 [7,8]). The main flow of the buoyant jet is assumed that remains symmetrical to the jet nξ plane. Therefore, the dimensionless mean axial velocity and concentration profiles are assumed Gaussian in the zone of established flow (ZEF): where b φ = K φ ξ is the nominal width of the buoyant jet where the property φ reduces at the e −1 of the maximum value φ m and K φ is the spreading rate coefficient. Thus, φ stands for either the mean axial velocity w or the mean concentration c.
In the zone of flow establishment (ZFE), the EMA model incorporates the advanced integral model (AIM; [85]). Thus, the actual similarity profiles of velocity and concentration are yielded by composing the flow and mixing fields of a group of infinite vertical buoyant jets issued from closely spaced sources on the basis of flux conservation of momentum, buoyancy and kinetic energy for the mean motion and enhancing the dynamic adaptation of the individual buoyant jet axes to the group centreline.
The entrainment velocity u is calculated using the transverse velocity profile hu/w m that is produced by integrating the equation of continuity (1) with respect to distance r [9,74]: According to [74], the vertical velocity w B of the escaping masses at the inner boundary is estimated as: where p = 0.144, K w = 0.11, λ Bp = 1.15, λ p = 1.04, R p = 0.3521 and Y p = 1.001 [74]. Additionally, the concentration c B of the escaping fluid by the inner boundary is assumed proportional to the corresponding centerline concentration as: where Λ is a constant that is determined approximately to 0.34 [74].

Variation of Basic Parameters
Integration of the PDE of tracer mass (4), under the assumptions (5)-(8) yields the ordinary differential equation (ODE) of tracer mass d dξ dξ Λc m w B cos θ, which gives the longitudinal variation of the slop of the buoyancy flux [74]. Figure 2 shows the variation of the normalized buoyancy flux with respect to the dimensionless distance Ξ. It is apparent that the buoyancy flux reduces along centerline path until Ξ 10 and then stabilizes to a constant value. This reduction increases via initial inclination angle θ 0 until θ 0 40 • -45 • and then decreases up to 75 • . Thus, as seen in Section 4, the distance y cr increases until θ 0 40 • -45 • and then decreases up to a minimum value approaching zero for vertical discharges. After the flow bends over, beyond the cross section where the maximum zct occurs, the buoyancy flux becomes positive. The major part of this region is beyond Ξ  10 ( Figure 2).
Therefore, the stability of the buoyancy flux means that the assumption made regarding the escaping mass has no effect in the region of positive buoyancy flux. The longitudinal variation of the local angle θ of the buoyant jet axis to horizontal with respect to the dimensionless distance Ξ = ξ/(d0F0) for various initial inclination angles θ0 (θ0 = 15°, 30°, 45°, 60°, 75°) and initial Froude numbers F0 (F0 = 5, 20) is shown in Figure 3. It is apparent that the effect of F0 to θ is rather insignificant for practical applications.

Results and Discussion
PDEs (1)-(4) are integrated with respect to r and ϕ on the actual cross-section of a round buoyant jet under the aforementioned boundary conditions, the spreading concept [9,85] and all the appropriate approximations [74]. The integration produces a system of ordinary differential equations (ODE) that is solved using a 4th order Runge-Kutta method [86]. More details about the procedure may be found at the relative paper of Yannopoulos and Bloutsos [74].
In this paper, EMA is applied to predict the mean flow geometrical quantities and mixing characteristics for inclined negatively turbulent buoyant jets for initial inclination angles 15° ≤ θ0 ≤ 75° and initial densimetric Froude numbers 4 ≤ F0 ≤ 100. For each angle, in the range 15° ≤ θ0 ≤ 75°, the variation of the non-dimensional geometric quantities of zct/d0, zet/d0, yct/d0, ycr/d0 and yer/d0 and mixing characteristics Sct and Scr are calculated for every initial Froude numbers 4 ≤ F0 ≤ 100. For After the flow bends over, beyond the cross section where the maximum z ct occurs, the buoyancy flux becomes positive. The major part of this region is beyond Ξ 10 ( Figure 2). Therefore, the stability of the buoyancy flux means that the assumption made regarding the escaping mass has no effect in the region of positive buoyancy flux.
The longitudinal variation of the local angle θ of the buoyant jet axis to horizontal with respect to the dimensionless distance Ξ = ξ/(d 0 F 0 ) for various initial inclination angles θ 0 (θ 0 = 15 • , 30 • , 45 • , 60 • , 75 • ) and initial Froude numbers F 0 (F 0 = 5, 20) is shown in Figure 3. It is apparent that the effect of F 0 to θ is rather insignificant for practical applications. After the flow bends over, beyond the cross section where the maximum zct occurs, the buoyancy flux becomes positive. The major part of this region is beyond Ξ  10 ( Figure 2).
Therefore, the stability of the buoyancy flux means that the assumption made regarding the escaping mass has no effect in the region of positive buoyancy flux. The longitudinal variation of the local angle θ of the buoyant jet axis to horizontal with respect to the dimensionless distance Ξ = ξ/(d0F0) for various initial inclination angles θ0 (θ0 = 15°, 30°, 45°, 60°, 75°) and initial Froude numbers F0 (F0 = 5, 20) is shown in Figure 3. It is apparent that the effect of F0 to θ is rather insignificant for practical applications.

Results and Discussion
PDEs (1)-(4) are integrated with respect to r and ϕ on the actual cross-section of a round buoyant jet under the aforementioned boundary conditions, the spreading concept [9,85] and all the appropriate approximations [74]. The integration produces a system of ordinary differential equations (ODE) that is solved using a 4th order Runge-Kutta method [86]. More details about the procedure may be found at the relative paper of Yannopoulos and Bloutsos [74].
In this paper, EMA is applied to predict the mean flow geometrical quantities and mixing characteristics for inclined negatively turbulent buoyant jets for initial inclination angles 15° ≤ θ0 ≤ 75° and initial densimetric Froude numbers 4 ≤ F0 ≤ 100. For each angle, in the range 15° ≤ θ0 ≤ 75°, the variation of the non-dimensional geometric quantities of zct/d0, zet/d0, yct/d0, ycr/d0 and yer/d0 and mixing characteristics Sct and Scr are calculated for every initial Froude numbers 4 ≤ F0 ≤ 100. For

Results and Discussion
PDEs (1)-(4) are integrated with respect to r and φ on the actual cross-section of a round buoyant jet under the aforementioned boundary conditions, the spreading concept [9,85] and all the appropriate approximations [74]. The integration produces a system of ordinary differential equations (ODE) that is solved using a 4th order Runge-Kutta method [86]. More details about the procedure may be found at the relative paper of Yannopoulos and Bloutsos [74].
In this paper, EMA is applied to predict the mean flow geometrical quantities and mixing characteristics for inclined negatively turbulent buoyant jets for initial inclination angles 15 • ≤ θ 0 ≤ 75 • and initial densimetric Froude numbers 4 ≤ F 0 ≤ 100. For each angle, in the range 15 • ≤ θ 0 ≤ 75 • , the variation of the non-dimensional geometric quantities of z ct /d 0 , z et /d 0 , y ct /d 0 , y cr /d 0 and y er /d 0 and mixing characteristics S ct and S cr are calculated for every initial Froude numbers 4 ≤ F 0 ≤ 100. For every θ 0 , these quantities vary proportionally to the Froude number by a constant k i (among others: [49,53,54]) as: where k i (i = 1, . . . , 7) is the proportionallity constant. The results are compared to experimental data available in the literature published by several researchers, as well as to other numerical models predictions such as CorJet, VisJet, the reduced buoyancy flux (RBF) model by Oliver et al. [58], modified RBF by Crowe et al. [61] and analytical solutions proposed by Kikkert et al. [49]. Available data from CFD analysis [76,78,81,82] are also used for comparisons. Figure 4 shows the variation of dimensionless centerline terminal rise height z ct to the initial inclination angle. EMA's prediction is compared to available experimental data. The terminal rise height EMA passes through the scatter of experimental data following their trend. The latter is more obvious if a polynomial fit to experimental data [45,47,[49][50][51]54,[57][58][59]61,73,81,87,88] is made. For this purpose, a third degree polynomial is fitted to experimental data using the least square method. The equation of the polynomial is presented in Table 1. In the range 15 • ≤ θ 0 ≤ 75 • EMA's prediction is very close to the polynomial fit. It must be noted that EMA predicts the centerline as the midpoint of the round cross-section. The same procedure is followed by Bosanquet et al. [87] and Bashitialshaaer et al. [88]. Instead, Cipollina et al. [45], Kikkert et al. [49], Ferrari and Querzoli [47], Shao and Law [50], Lai and Lee [51], Oliver et al. [58], Nikiforakis et al. [89], Crowe et al. [61], Ramakanth [73], and Zhang et al. [81] determine the position of centerline as the point where the transverse profile of velocity or concentration is maximized. In particular, this point is shifted towards the external edge, moving away from the midpoint of cross-section. As the initial inclination angle is rising, the deviation of actual transverse profile is deviates intensively from the Gaussian presenting this difference, while, when the initial inclination reaches the vertical, this deviation is reduced, and EMA's prediction coincides with the experimental data. Thus, the deviation between EMA's prediction and the experimental data in the range 20 • < θ 0 < 70 • is reasonably expected. In Figure 4b, EMA's prediction of centerline terminal rise height is compared to other models results. For θ 0 < 30 • all models almost coincide. For θ 0 > 30 • , Kikkert's analytical solution and the RBF model are diverging upwards calculating higher values than EMA's and modified RBF's predictions that approximately coincide up to 75 • . The commercial packages CorJet and VisJet are underestimating the centerline terminal rise height for θ 0 > 30 • . In Figure 4b, the results of CFD analysis from Gildeh et al. [78,79] Oliver et al. [76] and Zhang et al. [81,82]    The geometric parameter of the final rising height of the external edge (boundary) zet of the dense jets is of great environmental importance, because it indicates whether the mixing processes are taking place under the water surface or not. The dimensionless values of zet that are calculated by the EMA model through the initial inclination angles are shown on Figure 5a,b. In Figure 5a, the EMA's results are compared to previous published experimental data, and in Figure 5b, the The geometric parameter of the final rising height of the external edge (boundary) z et of the dense jets is of great environmental importance, because it indicates whether the mixing processes are taking place under the water surface or not. The dimensionless values of z et that are calculated by the EMA model through the initial inclination angles are shown on Figure 5a,b. In Figure 5a, the EMA's results are compared to previous published experimental data, and in Figure 5b, the performance of EMA is compared to other models. Again EMA, without intense deviation, follows the trend of a third degree polynomial fitted to experimental data [38,[40][41][42][43][44][45][48][49][50][51]53,56,58,60,69,73,81,[88][89][90][91][92][93] (Table 1). In Figure 5a, the experimental data of Nemlioglou and Roberts [48] and Bloomfield and Kerr [44] are the upper and lower data of the scatter, respectively. The wide scatter of experimental results is due to the different definitions of the external edge of the buoyant jet among the investigators. Indicatively, Kikert et al. [49] and Oliver et al. [58,78] define the external boundary at a distance twice the concentration spread, while Papakonstantis et al. [54] define this length as 1.5 times the concentration spread. Lai and Lee [51] and Abessi and Roberts [60] define the boundary at a locus where concentration values c are 25% or 10% the maximum concentration (c max ) at centerline, respectively. Similarly, Jiang et al. [52] and Zhang et al. [81] define the external boundary at c/c max = 5% and Shao and Law [50] and Gildeh et al. [78] define the external boundary at c/c max = 3%. According to EMA, the outer boundary of the jet is defined at a distance B c = d 0 /2 + n c λ b w , where b w is the nominal spread width of the velocity field of the buoyant jet, λ is the concentration-to-velocity spreading rate coefficient ratio and n c is the non-dimensional total spread of concentrations of a buoyant that according to profile measurements by Ramaprian and Chandrasekhara [94] and Shao and Law [50] is approximately 1.9. More details can be found at the relative paper of Yannopoulos and Bloutsos [74]. Among EMA's prediction and Kikkerts's solution, the RBF and modified RBF models there are slight differences up to approximately 70 • . A little deflection occurs between EMA's and CorJet's prediction for the cutoff level of c/c max = 3%. As VisJet's boundary is defined at 0.25 c max , its results are like CorJet's cutoff level of 25%, and obviously differ to the EMA relative values.
In Figure 6a,b, the dimensionless horizontal distance to centerline terminal rise height y ct predicted by EMA is shown. The distance y ct is increasing gradually up to initial angle of 45 • , approximately, and then decreases smoothly up to 75 • . Similar behavior is observed at experimental data of various previous works demonstrating the good performance of EMA. Excluding the experiments of Bosanquet et al. [87] and Lindberg [43], the experimental scatter is quite narrow, and EMA predictions are within the experimental data for the range 15 • ≤ θ 0 ≤ 75 • . CFD analyses provide data only for initial angles of 15 • , 30 • , 45 • and 60 • . The scatter of CFD data seems to be wider than experiments. Both CorJet and VisJet diverge from experimental data, underestimating the horizontal distance of centerline peak. EMA's predictions underestimate the corresponding experimental data, following the trend of a third degree polynomial fitted to experimental data [43,45,[49][50][51]53,56,58,61,73,81,87,89] (Table 1).
The dimensionless path length of the centerline up to terminal height sct ( Figure 11) almost coincides to corresponding experimental data of Crowe [69], Oliver [68] and Ramakanth [73], where the scatter of experimental data is negligible.
The only existing experimental data of the path length of external edge of negatively buoyant plumes until the return point are provided by Christodoulou and Papakonstantis [95], and are compared to the corresponding predictions of EMA model in Figure 12. It is obvious that the EMA and experimental results almost coincide.    The dimensionless path length of the centerline up to terminal height s ct (Figure 11) almost coincides to corresponding experimental data of Crowe [69], Oliver [68] and Ramakanth [73], where the scatter of experimental data is negligible.    The only existing experimental data of the path length of external edge of negatively buoyant plumes until the return point are provided by Christodoulou and Papakonstantis [95], and are compared to the corresponding predictions of EMA model in Figure 12. It is obvious that the EMA and experimental results almost coincide.

Conclusions
A second order Gaussian integral model that incorporates the mechanism of escaping masses from an inclined buoyant jet (EMA) was applied to predict the mean flow and mixing properties of inclined turbulent negatively round buoyant jets in stationary and uniform ambient. The predictions were made for a range of initial inclination angle from 15 • to 75 • to the horizontal and initial Froude number from 4 to 100. EMA's dimensionless results were compared to experimental data and other integral or numerical predictions available in the literature. For comparison reasons, a third degree polynomial is fitted to experimental data for every parameter under study. EMA's prediction of dimensionless terminal height of external edge almost coincides to the experimental fit, while the corresponding prediction of centerline terminal height differs slightly from the polynomial fit, for initial angles θ 0 from 35 • to 40 • . This difference is attributed to the way that centerline is defined by each researcher. EMA's predictions of horizontal distance of centerline terminal height are lower than the experimental fit polynomial, but within the range of experimental values. Regarding the horizontal location of return point of centerline and external edge, EMA's predictions follow the corresponding polynomial fits presenting a slight difference for the case of centerline variation. As an overall observation, EMA's predictions agree well to the experimental data, which assures that EMA performs better than other models. The definition of buoyant jet's centerline affects the estimated variation of axial dilution at centerline terminal height regarding the available experimental data, leading to overestimations for initial inclinations less than 45 • . Additionally, EMA's prediction of the axial dilution values at centerline return point is rather conservative compared to the corresponding experimental fit, but much closer than other models' predictions. Finally, EMA's of centerline length to terminal height and external edge length to return point coincides excellently to available experimental data. Overall, it seems that the application of EMA model gives reliable predictions without heavy computational cost or adopting controversial assumptions.