Effect of Proximity to Bed on 30° and 45° Inclined Dense Jets: a Numerical Study

Employing inclined dense jets is a common way for the disposal of brine from coastal desalination plants. Reaching an optimal design to minimize the negative impacts of brine on marine ecosystems has been a popular environmental-related line of research over the last decades. This paper numerically analyzes the mixing and geometrical properties of 30° and 45° inclined dense jets when they discharge close to the bed. For this purpose, two series of numerical simulations were developed. In the first series, the nozzle was placed far enough from the lower boundary to act as a free jet. Meanwhile, in the second series, the distance between the nozzle tip and seabed was substantially reduced. By comparing these two series, the effect of proximity to bed on the behavior of dense jets has been investigated. The governing equations were solved by modifying a solver within the CFD package of OpenFOAM. The numerical results were presented in comparative figures and were compared to previous works. Comparisons indicated that the numerical model predicts the geometrical characteristics of dense jets in good agreement with the past experimental studies. However, the dilution predictions are conservative. It has been observed that proximity to the bed has almost no appreciable effects on the behavior of 45° jets. However, for 30° jets, when the bed proximity parameter (y0/LM) falls below 0.14, the normalized values of horizontal and vertical locations of centerline peak and return point dilution are slightly reduced while the terminal rise height remains untouched. • The properties of free and boundary-affected dense jets were numerically analyzed. • Some characteristics of 30° jets were found to be influenced by the lower boundary. • The 45° dense jets were insensitive to variations of the bed proximity parameter. • The model predicted the geometric characteristics with reasonable accuracy. • The dilution predictions were found to be conservative.


Introduction
In recent decades, seawater desalination has attracted a lot of attention as a new and unconventional resource of fresh water. The main by-product of the desalination process is a concentrated effluent called brine, which is commonly disposed of into the sea (Saeedi et al. 2012). The disposal of this effluent has raised serious concerns due to its potential impacts on marine biota, especially benthic communities. Submerged discharge in the form of inclined dense jets is a common way to mitigate brine impacts. By using the inclined dense jets, brine efficiently mixes with the surrounding water, and consequently, its concentration reduces to the accepted levels for the marine ecosystems (Abessi andRoberts 2015, 2018).
Predicting the behavior of flow is possible through the experimental, analytical, and numerical models. Zeitoun et al. (1970) conducted several experiments on the effect of nozzle orientation on dense discharges. They reported 60°angle as the most appropriate inclination since it produces the longest trajectory and the highest dilution compared to 30°, 45°, and 90°j ets. Afterward, 60°angle was accepted as the de facto standard for brine discharge, although this has been questioned in later studies. Based on the latter suggestion, Roberts et al. (1997) carried out detailed measurements on inclined dense jets at 60°angle using Laser-Induced Fluorescence (LIF) and a micro-conductivity probe. They derived flow properties along the near-field and proposed some experimental coefficients for flow parameters such as impact point and the end of the initial mixing zone. Kikkert et al. (2007) performed experiments using Light Attenuation (LA) and LIF techniques to validate their developed analytical solution. Their comparisons indicated that the analytical solution reasonably predicts the flow path and terminal rise height, but the minimum dilution at the location of the return point was underestimated. Lai and Lee (2012) reported a comprehensive series of experimental works on dense discharge in stationary water. They employed LIF and Particle Image Velocimetry (PIV) techniques to measure tracer concentration and velocity fields, respectively. According to their observations, the dimensionless maximum rise height y t /(d. Fr d ) is independent of source conditions for Fr d ≥ 25, and the impact point dilution is not sensitive to nozzle angle within the range of 38°-60°. Abessi and Roberts (2015) conducted experimental studies on single-port dense jets oriented at angles from 15°to 85°to the horizontal using the threedimensional LIF (3D-LIF) technique. The experimental results showed that the dilution at the impact point is almost insensitive to nozzle angle over the range of about 45°-65°, while the near-field dilution is more sensitive to nozzle inclination and is highest for 60°angle. Furthermore, they inquired about the bottom boundary effects on the impact point dilution. They observed that the time-averaged dilution along the jet centerline first increased and then decreased in a thin layer close to the wall. It was found to be due to an increase in turbulent intermittency and the accumulation of more saline patches of fluid at the bed. Papakonstantis and Tsatsara (2019) also carried out experiments on inclined dense jets utilizing a micro-scale conductivity probe for concentration measurements. They reported mixing characteristics of 35°, 50°, and 70°jets, which had been rarely investigated in previous studies, and reported data compatible with what was observed previously for other angles.
Although the 60°inclined dense jets have been widely accepted for brine discharge from desalination plants, the associated terminal rise height is relatively high for shallow coastlines; therefore, smaller nozzle angles and lower nozzles seem more practical in those regions (Abessi and Roberts 2016). Shao and Law (2010) carried out comprehensive experiments on the mixing behavior of dense jets at smaller angles of 30°and 45°using combined PIV and Planar LIF (PLIF) techniques. The effect of proximity to the bed on the flow behavior was also examined. They reported that in 30°dense jets, the flow is significantly affected by the bottom boundary for bed proximity parameter (y 0 /L M ) values less than 0.15. Abessi and Roberts (2016) have also performed experiments on inclined dense jets with nozzles oriented at 30°, 45°, and 60°in shallow water. They proposed the 30°jet for shallow water regime since it has fewer surface interactions and visual impacts and gives slightly better dilution at the impact point compared to other angles.
Besides experimental and analytical studies, various numerical studies have also been reported in the context of dense jets. Vafeiadou et al. (2005) used CFX-5 of ANSYS software with the k − ω Shear Stress Transport (SST) turbulence model to numerically simulate 60 o inclined dense jet in stagnant water. They found that the model underestimated the terminal rise height and the horizontal distance of the return point compared to the experimental data of Roberts et al. (1997). Oliver et al. (2008) investigated the initial mixing of negatively buoyant jets using the ANSYS CFX model. They tried to improve the predictions by adjusting the turbulent Schmidt number in the tracer transport equation for positively buoyant vertical jets and then applying it to inclined negatively buoyant jets. However, the overall effect of this approach on the quality of the results was small, and the dilution predictions were underestimated. Kheirkhah Gildeh et al. (2015) studied the geometrical and mixing characteristics of 30°and 45°inclined dense jets in stationary ambient using the OpenFOAM model. They evaluated the effect of several turbulence models on the accuracy of CFD predictions. By comparing the numerical results to previous experimental data, they indicated that the LRR and realizable k − ε turbulence models are able to predict the flow behavior more accurately. Zhang et al. (2017) performed a numerical study on 45°and 60°inclined dense jets in stagnant water with the Large Eddy Simulation (LES) approach. Their simulations could reproduce concentration build-up at the impact point reported by Abessi and Roberts (2015); however, dilutions at the impact point and within the spreading layer were underpredicted. In a recent study, Tahmooresi and Ahmadyar (2021) examined the effects of the turbulent Schmidt number on CFD predictions of 45°inclined dense jets. They reported that reducing the turbulent Schmidt number from 1.0 to 0.4 will improve the dilution predictions, whereas this change adversely affects the geometrical parameters and cross-sectional distributions of concentration.
The review of previous works highlights the absence of comprehensive numerical work for predicting the effect of bottom boundary on the behavior of the dense jets when they discharge close to the bed. There are generally three approaches to numerically predict turbulent flows, i.e., Direct Numerical Simulation (DNS), Large Eddy Simulation (LES), and Reynolds-Averaged Navier-Stokes (RANS). Although DNS and LES approaches have substantially more resolution than RANS, the latter is used considerably more widely due to its lower computational expenses. This paper numerically investigates the effects of bed proximity on mixing and geometrical characteristics of 30°and 45°inclined dense jets with the RANS approach. For this purpose, two series of numerical models are considered. In the first series, the nozzles are placed far enough from the floor to act as a free jet. Meanwhile, in the second series, the distance of nozzles from the bed is substantially reduced. The governing equations are implemented and solved within the OpenFOAM finite volume model, and the realizable k − ε turbulence model is employed for turbulent closure. The results of these two numerical series are compared to each other and also previous experimental studies. Through this comparison, assessing the probable effects of proximity to bed on the flow behavior would be possible.

Dimensional Analysis
The schematic side view of an inclined dense jet is shown in Fig. 1. The jet discharges upwardly from a round nozzle of diameter d, with jet velocity U 0 , jet density ρ 0 , and at the initial angle of θ to the horizontal. The jet rises and mixes with ambient water (ρ a < ρ 0 ) due to its initial momentum until it reaches a maximum height (y t ), and then it falls back to the floor because of its negative buoyancy. Finally, it impacts the bed at a horizontal distance of x i from the source and continues spreading as a density current.
The dimensional analysis of free jets in deep water is well-known and is described in Fischer et al. (1979), Roberts and Toms (1987) and Roberts et al. (1997), but for jets placed at a close distance to the bed, the height of the nozzle tip from the bed y 0 should be added. Based on the mentioned references, for fully turbulent jets in which the flow is independent of viscosity, all the dependent variables of the flow φ, could be characterized by the discharge angle θ, the jet kinematic fluxes of volume Q 0 , momentum M 0 , buoyancy B 0 , and the height of the nozzle tip from the bed y 0 : Þ=ρ a is the modified acceleration due to gravity. These fluxes can form the jet-to-plume transition length scale: This length scale represents the distance from the source of discharge in which the flow is dominated by initial momentum flux (jet-like). This characteristic length scale is also equal to (π/4) 1/4 d Fr d where Fr d is the jet densimetric Froude number: For high Froude numbers (Fr d > 20), the initial volume flux Q 0 is not dynamically important, and it can be neglected (Roberts et al. 1997). Therefore, for boundary-affected jets at a specific angle θ, the dependent variables are a function of M 0 , B 0 , and y 0 : Following a dimensional analysis, any characteristic length χ, can be expressed as: and dilution as: where dilution is defined as s = (C 0 − C a )/(C − C a ) where C 0 , C a , and C are jet discharge concentration, ambient concentration, and local time-averaged concentration, respectively.

Governing Equations
The governing equations consist of the continuity equation, Navier-Stokes equations with the Boussinesq approximation for buoyancy effects, and tracer advection-diffusion equation for incompressible three-dimensional flows. Based on the Boussinesq approximation, when the density variation is not large, the density can be treated as a constant in the unsteady and convection terms and as a variable only in the gravitational term (Ferziger and Perić 2002). After applying the Reynolds-averaging, the governing equations take the following form: Continuity equation: Navier-Stokes equations with the Boussinesq approximation for buoyancy: where U i and u 0 i are mean and fluctuating part of the fluid velocity, P is the mean flow hydrodynamic pressure, ρ is the fluid density, and ν is the kinematic viscosity.
Tracer advection-diffusion equation: where C and c ′ are the mean and fluctuating part of the tracer concentration, D = ν/Sc is the molecular diffusion, and Sc is the Schmidt number. The density of both the jet and the ambient water is calculated using the equation of state for seawater proposed by El-Dessouky and Ettouney (2002).

Turbulence Modeling
The above RANS equations have closure problem due to the Reynolds stress tensor τ ij ¼ − u 0 i u 0 j and turbulent scalar fluxes u 0 j c 0 . The Reynolds stress tensor is commonly modeled using the Boussinesq hypothesis, which assumes that the Reynolds stress is a linear function of the mean velocity gradients (Moukalled et al. 2016): where ν t is the turbulent viscosity, is the turbulence kinetic energy, and S ij is the mean strain-rate tensor: With this assumption, it is required to calculate the turbulent viscosity and turbulence kinetic energy rather than the Reynolds stresses. Several turbulence models have been developed to compute these two terms. The present paper employs the realizable k − ε model of Shih et al. (1995), which has been used and validated in several studies in the context of dense jets (Kheirkhah Gildeh et al. 2015; Ardalan and Vafaei 2019; Tahmooresi and Ahmadyar 2021). The formulation of the realizable k − ε turbulence model is given by: where G k ¼ 2 ν t S 2 ij represents the generation of turbulence kinetic energy due to the mean velocity gradients, D k = ν + (ν t /σ k ) is the effective diffusion coefficient for k, and D ε = ν + (ν t / σ ε ) is the effective diffusion coefficient for the rate of dissipation of turbulence energy ε. In the model, the turbulent viscosity is obtained using: and C μ is computed as follows: e where Ω ij is the mean rate of rotation tensor and ω k is the angular velocity. The parameters A 0 and A s are determined as follows: C 1ε is defined as follows: The model constants were assigned the following values: C 2ε = 1.9, σ k = 1.0, and σ ε = 1.
where Sc t is the turbulent Schmidt number.

Flow Configuration and Computational Setup
Based on the available experimental setups, a computational domain with dimensions of 1.25 m length, 0.4 m width, and 0.4 m depth was chosen, but only half of the domain (1.25 m × 0.2 m × 0.4 m) was considered since the problem is symmetrical. The depth and width of the considered domain satisfy the criteria for the deep flow regime (dFr d /H < 0.8) and acting as a single jet (s/(d • Fr d ) > 2, where s is port spacing) proposed by Abessi andRoberts (2014, 2016). For grid independency analysis, the dilution at the return point was considered as the monitoring point with a 2% tolerance. It was found that a grid of 650,000 cells with increasing grid spacing from the center of the nozzle to the boundaries, as shown in Fig. 2, provides independence of results from the grid cell sizes.
The boundary conditions for the nozzle inlet are as follows: The initial values for k and ε are calculated using the equation proposed by Huai et al. (2010) as follows: At the outlet section, a zero-gradient boundary condition perpendicular to the outlet plane is determined for variables U, k, ε, and C, while a Dirichlet boundary condition with a constant value of zero is applied to P to set a reference pressure. For the wall boundaries, the no-slip condition is imposed for velocity. Meanwhile, the OpenFOAM wall functions for low-and high-Reynolds number turbulent flow cases, namely kLowReWallFunction and epsilonWallFunction are used for k and ε, respectively. Moreover, at these sections, the fixedFluxPressure boundary condition is applied to P, which adjusts the pressure gradient so that the boundary flux matches the velocity boundary condition (Greenshields 2015). Finally, for the symmetry plane, the symmetry boundary condition is employed. This boundary condition is imposed by setting normal gradient of scalars and velocity vector as well as normal velocity to zero (Moukalled et al. 2016).
The simulations are performed using the open-source software package of OpenFOAM. The modified buoyantBoussinesqPimpleFoam solver, which is a transient solver for buoyant, turbulent, incompressible flows, is employed within OpenFOAM. This solver uses the Boussinesq approximation for buoyancy effects and the PIMPLE algorithm to solve the coupled pressure-momentum system. The PIMPLE algorithm is a combination of SIMPLE and PISO algorithms, which allows using of larger Courant numbers and time intervals. In the present study, however, the solver is run in PISO mode and also with the Courant number less than unity to ensure that any important transient information, which can significantly influence the flow, would not be missed. Further details about the PIMPLE algorithm can be found in Holzmann (2018).
The first-order implicit Euler scheme is used for temporal discretization. The advection (divergence) and diffusion (laplacian) terms are discretized using the standard Gaussian finitevolume integration in which the interpolation of values on cell faces from cell centers is required. The limitedLinear and linearUpwind interpolation schemes are employed for the advection terms of scalar variables and velocity, respectively, and also the linear scheme for the diffusion of all variables. The discretized equations form linear sets of equations as Ax = b, which have to be solved using a linear system solver method. For this purpose, the Preconditioned Conjugate Gradient (PCG) method with the Diagonal Incomplete Cholesky (DIC) preconditioner is used for solving the system of linear equations of the pressure field. Meanwhile, the Preconditioned Biconjugate Gradient (PBiCG) method with the Diagonal Incomplete LU (DILU) preconditioner is employed for the other fields. The convergence criterion of 10 −7 is set for U, k, ε, and C, while that for P is 10 −9 . The Sc and Sc t are considered constant with the values of 600 and 0.7, respectively (Rard and Miller 1979;Nimdeo et al. 2014;van Reeuwijk et al. 2016).

Results and Discussion
In total, twelve numerical simulations in two series have been performed to investigate the possible effects of bed proximity on flow behavior. The characteristics of these simulations are presented in Table 1. In the table, the far and proximate-to-bed jets have been denoted by F (stands for far) and N (stands for near), respectively.

General Observations
In most cases, the flow reached the steady-state condition about 30-40 s after the beginning of discharge. The simulations, however, were continued until 80 s after the start time. Unlike the DNS and LES approaches, all the turbulence fluctuations are modeled and represented in terms of the mean-flow characteristics in the RANS approach. Hence, after the flow reaches the steady-state condition, almost no appreciable difference is seen between the instantaneous and time-averaged values of a quantity at a specific point. While in real conditions, flow-dependent quantities at a specific point fluctuate widely; that is, the instantaneous values can be substantially higher or lower than the time-averaged value. Figure 3 shows the non-dimensional time-averaged concentration contour maps and lines of two cases in far and near series for both inclination angles. The selected cases from each angle have identical densimetric Froude numbers but different bed proximity parameter (y 0 /L M ). The figure does not show a significant difference in the flow behavior for these two cases. Like the free jet, the jet that is placed close to the bed rises and mixes with the ambient water until a maximum height and then falls back to the floor. As expected, due to the mixing process, the concentration of flow is considerably reduced when it impacts the seafloor.
The centerline trajectory is one of the simplest and most important characteristics of the flow, which is derived through the connection of the maximum velocity or concentration location at each cross-section perpendicular to the flow. It is worth noting that the trajectories obtained using the velocity variable almost coincide with those obtained through the concentration variable, but the latter descends slightly faster (Shao 2010). The centerline trajectories for various cases that were derived using the concentration variable have been illustrated in   Baum and Gibbes 2020). It indicates that the common RANS-based approaches are not able to sufficiently resolve the effects of turbulence along the plume-like region of the flow where the detrainment and buoyancy-driven instabilities are significant.
The slope of the ascending and descending parts of the trajectories for both series are almost equal, meaning that the trajectories are relatively symmetrical. The symmetry of centreline trajectories of free jets is in agreement with the observations of Shao and Law (2010). They also reported a slight asymmetry in boundary-affected cases, which is not seen herein. For 30°j ets, the major difference between the trajectories of two series F and N is that the trajectories of series N descend sooner in comparison to series F. This difference, which is also reported by Shao and Law (2010), happens due to the effect of Coanda on the flow behavior. The Coanda effect is the tendency of a fluid to adhere to a surface because of the reduced pressure due to flow acceleration around the surface (Constantin 2010). However, no notable difference is seen between the two series for 45°jets. This means that nozzle angle plays an important role in determining the bed influence on the flow behavior in addition to distance from the bottom boundary. So, it can be concluded that the steeper the nozzle angle, the less the effect of the boundary on the flow behavior is.

Centerline Peak
The centerline peak (x m , y m ) is the location where the centerline trajectory reaches its maximum rise height. The horizontal and vertical locations of this point are normalized with the jet diameter and plotted versus Fr d in Fig. 5. The experimental results from Cipollina et al. (2005), Kikkert (2006), and Shao and Law (2010), along with the analytical model of Kikkert (2006), are added for comparison. In Fig. 5, the numerical simulations slightly underestimate the location of the centerline peak in comparison to previous experimental data and the analytical model of Kikkert (2006). In free jets, the location of this point is controlled by the initial momentum flux and nozzle angle. It means that the centerline peak is situated in a region of flow, where the flow behavior is jet-like, not plume-like. It was reported by Ghayoor et al. (2019) after analyzing the turbulent fluctuations at different locations within the jet-and plume-like regions of dense jets. To perceive how proximity to the bed affects the location of the centerline peak, the numerical results are normalized with L M and plotted versus the bed proximity parameter in Fig. 6. No appreciable changes are observed in the centerline peak location with the variation of bed proximity parameter for 45°jets. For 30°jets, however, when y 0 /L M ≤ 0.14, the normalized values of both the horizontal and vertical location of the centerline peak are reduced compared to the free condition. It is worth mentioning that the reduction in the horizontal component is somewhat higher than the vertical component. The transition criterion of y 0 /L M ≤ 0.14 between free and boundary-affected 30°jets is in agreement with Shao and Law (2010), who reported y 0 /L M ≤ 0.15 for boundary-affected jets. As can be seen from Fig. 6a and c, the last case of series F, F4, which has the least bed proximity parameter in the series, has been placed on the line extending from series N. Meanwhile, all cases in series F have equal y 0 /d. Hence, it can be concluded that the bed proximity parameter y 0 /L M is a better deciding factor for the determination of the boundary effect than y 0 /d, which is also consistent with the presented dimensional analysis. This also means that both hydrodynamical and geometrical characteristics of discharge have roles in determining bed influence on the flow behavior.

Terminal Rise Height
As discussed earlier, the dense jet rises due to its initial momentum until it reaches a maximum height. At this point, which is known as the terminal rise height y t , the vertical component of momentum equals zero. This height plays a key role in the design of brine outfalls as it can lead to undesirable visual impacts and distinctive reduction in dilution (Jiang et al. 2014;Abessi 2018;Shrivastava and Adams 2019). There is a lack of consensus across past studies about the determination of terminal rise height. Lai and Lee (2012) considered the visual boundary as 0.25C max concentration contour, which corresponds to the radial position where turbulent intermittency (the fraction of time that a point was occupied by turbulent flow) is 0.5. The famous integral model of CorJet uses two cut-off levels of 3% and 25% for the visual boundary. In the present study, similar to Shao and Law (2010), the cut-off level of 3% is employed for the derivation of terminal rise height. Papakonstantis and Tsatsara (2018) considered a transient behavior for dense jets and reported two values for terminal rise height: an initial value for the developing stage of the flow and a final value for the quasi-steady stage. In the present study, however, we derived the terminal rise height based on the time-averaged concentration field after the flow reached steady-state. The terminal rise height of both series was normalized with the jet diameter and plotted against Fr d in Fig. 7. The numerical results are in acceptable agreement with previous experimental data and slightly lower than the analytical solution of Kikkert (2006). In order to investigate the possible effects of bed proximity on terminal rise height, the numerical results are non-dimensionalized with L M and plotted versus the bed proximity parameter in Fig. 8. As seen in the figure, the terminal rise height is almost insensible to the variations of bed proximity parameter for both angles over the range investigated herein, which is in agreement with the experimental study conducted by Shao and Law (2010).

Horizontal Distance of Return Point
After the initial rise, the dense jet falls back to the floor due to its negative buoyancy at a location so-called impact point. As discussed by Roberts et al. (1997), the minimum dilution along the lower boundary occurs at the impact point. Beyond this point, the dilution increases until it reaches a maximum value, and then entrainment collapses due to the decay of turbulent fluctuations and relaminarization of the flow. They defined the end of the active mixing zone as the location where turbulent intensity falls below 5%. The ultimate dilution at the end of the near field is roughly 60% more than the impact point dilution (Roberts et al. 1997;Abessi and Roberts 2015). Additional mixing and dilution beyond the near field will be mostly because of ambient turbulence, which is significantly lower than the initial mixing of the jet. As mentioned above, the minimum dilution along the bed occurs at the impact point; hence, the dilution at the impact point is commonly used as the indicator of the negative impacts of brine discharge on the benthic organisms. However, the location of the impact point depends on nozzle elevation and bed slope, so it is site-specific. Therefore, various past studies (Shao and Law 2010;Christodoulou et al. 2015;Crowe et al. 2016;Papakonstantis and Tsatsara 2018) have investigated the return pointwhich is independent of nozzle height and bed sloperather than the impact point for more generality. The return point is the location where the flow returns to the nozzle elevation. However, in most practical cases, the nozzle height and bed slope are typically small relative to the entire mixing zone. Thus, in these cases, the difference between the return and impact points would not be remarkable.
The horizontal distance of the return point after normalizing with the jet diameter is plotted versus Fr d in Fig. 9. Similar to the previously investigated geometrical characteristics, the numerical results are in good agreement with experimental data and slightly lower than the analytical model of Kikkert (2006). In Fig. 10, the derived dimensionless horizontal distances are plotted against the bed proximity parameter. No meaningful changes are observed in horizontal distances with the variations of bed proximity parameter for both angles. While Shao and Law (2010) reported that in boundary-affected 30°cases, the flow returns to the discharge source elevation at a slightly closer distance compared to the free condition.

Cross-Sectional Concentration Profile and Jet Spread
Mean centerline concentration profiles at three locations along the jet trajectories for one case of both series are demonstrated in Fig. 11. The profiles are perpendicular to the flow path and are normalized as C/C C against r/b C where C C , r, and b C are maximum local concentration at the cross-section, the radial distance, and the radius where C/C C = 1/e, respectively. These three sections are chosen so that the first section is in the jet-like region, the second one is at terminal rise height, and the last one is in the plume-like region of the flow. As seen in the figures, in the jet-like region where the flow is momentum-dominated, profiles are almost coincident with the Gaussian profile. By distancing from the nozzle along the trajectory, the inner (lower) half of the jet starts to deviate from the Gaussian profile, while the outer (upper) half remains Gaussian. This deviation from the Gaussian profile is caused by detrainment in the lower half, known as buoyant instability. The destruction of the axial symmetry of jets due to buoyant instabilities is also reported by Shao and Law (2010), Lai and Lee (2012), and Abessi and Roberts (2015). The 1/e width values of the upper half of the jets are extracted based on the normalized concentration profiles along the trajectory and plotted against the non-dimensional centerline length L C at the corresponding cross-section in Fig. 12a. The gradient of jet width to path length is typically used to characterize the spread of dense jets. It can be seen from the figure that the predicted growth rate is slightly higher than the reported value of Lai and Lee (2012) and lower than that of Jiang et al. (2014). However, as mentioned above, the cross-sectional profiles along the trajectory become asymmetric by distancing away from the source. Consequently, the inner and outer spread of dense jets is not equal to each other throughout the flow. The inner (lower) and outer (upper) widths of discharge and path length are normalized with d • Fr d and plotted against each other in Fig. 12b. The experimental results of Crowe et al. (2016) are also added for comparison. As can be observed in this figure, the inner and outer widths are almost identical near the source (from nozzle tip until about terminal rise height). Beyond this region, the inner and outer widths start to deviate from each other. The predicted lower widths are in good agreement with experimental data, but the upper width results are slightly bigger than experimental observations.

Minimum Dilution at Centerline Peak and Return Point
In order to investigate the possible effects of proximity to bed on mixing properties at centerline peak and return point, the dilution values of both 30°and 45°jets are normalized with corresponding densimetric Froude number and plotted versus y 0 /L M in Figs. 13 and 14, respectively. It is observed that the dilution predictions are close to the experimental values of Lai and Lee (2012) and Oliver et al. (2013), and considerably lower than those of Shao and Law (2010) as well as Abessi and Roberts (2015). The geometrical and mixing predictions of the present study, along with the results of various past experimental and analytical works, are quantified in Tables 2 and 3 for better comparison. It is seen in these tables that the dilution predictions of the present numerical model are superior than the well-known integral model of CorJet. The conservative dilution predictions of RANS models have also been reported by Zhang et al. (2017) and Robinson et al. (2016). This may be attributed to considering a constant turbulent Schmidt number throughout the fluid flow and the SGDH model in which the scalar flux vector is aligned with the mean scalar gradient vector (Pope 2000;Lai and Socolofsky 2019). As a result of the latter, the SGDH model is known to inaccurately predict the turbulent effects even in some simple turbulent flows (especially highly anisotropic and buoyant flows) (Combest et al. 2011). However, a noticeable variation in centerline peak or return point dilution is not observed with changes in bed proximity parameter for 45°jets. For 30°jets, however, it can be concluded that boundary-affected cases have a slightly lower dilution compared to free jets, which is in agreement with the observations of Shao and Law (2010). It is worth noting that the amount of reduction in their study is more than the present numerical predictions. The decline in dilution may be due to the reduced entrainment of ambient fluid in the lower boundary.  Using marine outfalls to form a dense jet is a common way to mitigate the environmental impacts of brine discharge from seawater desalination plants into the coastal waters. The main objective of the present study was to numerically investigate the mixing and geometrical characteristics of 30°and 45°inclined dense jets when they are placed at a close distance to the lower boundary. This condition is important when, in addition to reducing the nozzle angle, the designers need to reduce the nozzle height to prevent surface contact in shallow waters. For this purpose, two numerical series were performed by modifying a solver in the CFD package of OpenFOAM. In the first series, series F, the nozzle was relatively far from the bed to act as a free jet. In the second series, series N, the distance of the nozzle to the lower boundary was substantially reduced. Consequently, the effects of proximity to the bed became possible by comparing the results of these two series. The following conclusions can be drawn from the present study: 1. The bed proximity influences on the flow behavior depend not only on the nozzle height but also on its angle, and they are more significant in smaller angles. The proximity to the bed was found to have almost no appreciable effects on the flow behavior for 45°jets, but some changes in the geometrical and mixing characteristics of 30°jets were observed. 2. The bed proximity parameter, y 0 /L M , was found to be the controlling parameter for bed influence and not y 0 /d. For 30°jets, the influences of proximity to the lower boundary became evident when y 0 /L M ≤ 0.14. 3. Consistent with previous experimental studies, it was found that the geometrical and mixing characteristics of dense jets, including centerline peak, terminal rise height, the horizontal distance of return point, and the centerline peak and return point dilution, have a correlation with the densimetric Froude number. The proportionality coefficients were determined and presented in Tables 2 and 3. 4. In boundary-affected jets, the normalized values of horizontal and vertical locations of centerline peak were reduced compared to their counterparts in free jets. Simultaneously, the terminal rise height and the horizontal distance of the return point were almost  Kikkert et al. (2007) LA data 1.84 1.06 1.60 3.26 -- Kikkert et al. (2007) Model prediction 1.88 1.13 1.61 3.04 -- Papakonstantis et al. (2011a, b) Visual/Probe data --1.58 3.78 0.52 1.55 Abessi and Roberts (2015) 3D-LIF data -1.33 1.80 3.60 (i) 0.58 1.60 (i) Oliver et al. (2013) LIF data The denoted values by (i) are for the impact point. * The CorJet data are obtained from Palomar et al. (2012) untouched. The noted observations are in agreement with previous experimental data except for the horizontal distance of the return point. 5. The predicted geometrical properties matched well with previous experimental measurements; however, dilution predictions were conservative for both the free and boundaryaffected jets. 6. A slight reduction in the centerline peak and return point dilution was observed for boundary-affected jets, consistent with a past experimental study. However, they reported more reduction amount in comparison to the present numerical simulations.