Zonation of Positively Buoyant Jets Interacting with the Water-Free Surface Quantified by Physical and Numerical Modelling

The evolution of positively buoyant jets was studied with non-intrusive techniques—Particle Image Velocimetry (PIV) and Laser Induce Fluorescence (LIF)—by analyzing four physical tests in their four characteristic zones: momentum dominant zone (jet-like), momentum to buoyancy transition zone (jet to plume), buoyancy dominant zone (plume-like), and lateral dispersion dominant zone. Four configurations were tested modifying the momentum and the buoyancy of the effluent through variations of flow discharge and the thermal gradient with the receiving water body, respectively. The physical model results were used to evaluate the performance of numerical models to describe such flows. Furthermore, a new method to delimitate the four characteristic zones of positively buoyant jets interacting with the water-free surface was proposed using the angle (α) shaped by the tangent of the centerline trajectory and the longitudinal axis. Physical model results showed that the dispersion of mass (concentrations) was always greater than the dispersion of energy (velocity) during the evolution of positively buoyant jets. The semiempirical models (CORJET and VISJET) underestimated the trajectory and overestimated the dilution of positively buoyant jets close to the impact zone with the water-free surface. The computational fluid dynamics (CFD) model (Open Field Operation And Manipulation model (OpenFOAM)) is able to reproduce the behavior of positively buoyant jets for all the proposed zones according to the physical results.


Introduction
The economic growth throughout the 20th century fostered the establishment of industries near estuarine and coastal zones. Their activities demanded suitable areas for the disposal of their discharges, being buoyant jets. The presence of these discharges together with the high environmental requirements of these water bodies supposes a management challenge [1,2]. These wastewaters can contain physical, chemical and/or biological pollutants and cannot be adequately quantified without actual measuring and testing. Buoyant jets evolution is dominated by momentum and buoyancy. Buoyancy is a density-dependent property. Accordingly, a jet with a lower density than the receiving medium will have a positive buoyancy and will tend to rise to the surface due to gravity forces, known as positively buoyant jets or forced plumes. In the case of negative buoyancy, the discharge will have a higher density than the receiving environment so the jet will tend to descend to the bottom, known as negatively buoyant jets.
Industrial and urban wastewaters are typically freshwaters. They are characterized by a positive buoyancy and an initial momentum when discharging into brackish or salt waters such as estuarine differences between different tests to calculate the length scales for a given buoyant jet. Furthermore, precise information of vertical discretization is needed to analyze the evolution of a jet when its behavior is like a plume evolution and these data are not always available. For instance, this discretization introduces important errors of up to 30% in the jet centerline trajectory and up of 50% in jet centerline dilution [38].
Therefore, this study aims to assess the zonation of positively buoyant jets from its discharge into the receiving water body until its advance in the far-field by proposing a new methodology based on the angle (α) shaped by the tangent of the centerline trajectory and the longitudinal axis. α is a value from which it is possible to define which type of force is dominant over the jet during its advance in a fluid environment. To accomplish this work, the specific objectives are: (1) collecting physical data on the evolution (trajectory and dilution) of positively buoyant jets in a stagnant medium by physical modelling using a coupled PIV-LIF technique, (2) calibrating two semiempirical models (CORJET and VISJET) and one CFD model (OpenFOAM) with the collected data, and (3) applying the proposed method to assess the accuracy of numerical modelling to calculate the zone limits of positively buoyant jets.

Description of the PIV and LIF Techniques
PIV is a laser optical technique used for the measurement of instantaneous velocity flow fields. The flow is seeded with suitable tracer particles which are assumed to follow the flow faithfully. These particles scatter light when illuminated by means of a laser. A Charge Coupled Device (CCD) camera, synchronized to the laser, captures several images of the light scattered by each particle within a short time interval. The unit displacement of the particles can then be computed from successive images by means of cross-correlation. Finally, the local velocity is obtained by dividing the displacement over the time interval between the laser pulses [29].
LIF is a technique used in visualization experiments. However, if the experiments are performed carefully, it is possible to attain quantitative concentration measurements. This technique is based on the physical property of certain organic molecules, that allows an electron in the outer orbit to be exited. It moves from one of the vibrational levels of the ground state into the so-called singlet state. A laser is used to excite a fluorescent species within the flow. Typically, the tracer is an organic fluorescent dye such as fluorescein or rhodamine. The dye absorbs a portion of the excitation energy and spontaneously reemits a portion of the absorbed energy as fluorescence. The fluorescence is measured optically and used to infer the local concentration of the dye [39]. Accordingly, LIF system requires an appropriately paired CCD camera, laser and fluorescent dye combination, with the laser frequency within the absorption band of the dye.
The use of PIV and LIF techniques have an advantage over the classical methods when carrying out measurements of velocities and concentrations in a fluid environment because this measurement techniques have an optical nature and they are not intrusive on the flow. Thus, the magnitudes are not disturbed. On the contrary, classic measurements introduce an instrument that obviously influences the flow patterns, such as hot wire anemometry.

Setup of the Physical Modelling
The physical tests were performed at IHCantabria by means of simultaneous PIV and LIF techniques [40,41]. The physical setup comprised (see in Figure 1a): (1) a steel test tank with 3 m length, 3 m width, and 1 m depth, acting as the receiving water body, (2) a 100 L steel constant level tank containing the effluent, (3) a 1000 L plastic effluent storage tank, and (4) a piping system connecting the three tanks. The LIF camera has a band-pass filter, which allows only the fluorescent light emitted from the tracer dye excited by the laser to pass through. The filter ensures an excellent blocking efficiency of the excitation wavelength with a steep edge at 540 nm providing maximum transmission in its working range (545-800 nm). The measured images were realized with a total window size of 71 × 71 cm 2 and an acquisition rate of 5 Hz. For the PIV measurements, an interrogation window size of 32 × 32 pixel was used. The separation time (∆t) was 20 µs for the jet centerline path around the nozzle and 7 ms for the spreading layer where the jet presents a marked upward trajectory, according to [43].
Polyamide particles with a 50 µm diameter and density of 1030 kg/m 3 were used as seeding to follow the flow velocity variations. The standard cyclic Fourier Transform (FT) correlation function was applied to eliminate the PIV-generated noise and detect the peak that shows the velocity for each time step. For the LIF measurements, the exposure time was 20 ms. The absence of surfactants in the water-free surface was checked by liquid chromatography, since small amounts of them may alter the turbulence phenomena [44]. As passive tracer, Rhodamine WT (C29H29ClN2Na2O5) with a concentration of 100 µg/L was used. This concentration was selected to assure the right measurement in the highest dilution zones and to avoid any significant laser attenuation caused by self-absorption The test tank was made with two side glass windows. One window was used for illuminating the flow with a laser sheet and, the other window, for imaging the flow. Moreover, walls, bottom and the aluminum discharge tube were painted in black to avoid reflections of the laser during the tests. The nozzle of the discharge tube had a 5 mm inner diameter (D) and was located 50 cm above and parallel to the bottom of the test tank. The discharged flow was measured by an ultrasonic flow meter. The constant level tank contained a computer-controlled thermostat as a heating device (DIGITERM 2000, Selecta) and was located on top of a 5 m structure to obtain a gravity-driven flow as shown in Figure 1a. Furthermore, a recirculating pump was placed inside of the constant level tank to ease the fluid homogenization and to avoid the formation of thermoclines. Finally, the piping system and the walls of the constant level tank were insulated with neoprene to reduce heat losses.
Regarding the PIV and LIF techniques, a Q-switched double pulse Nd-Yag laser was used to illuminate the flow, using a Quantel model Twins BSL140 dual-cavity laser to create a bidimensional sheet, with a thickness between 0.5 and 2.5 mm. The energy level was 130 mJ per pulse and the pulse duration about 8 ns. Light emitted was green with a 532 nm wavelength and with a 30 Hz pulse repetition rate of each cavity. The laser was equipped with a telescopic arm for vertically placing the laser sheet and passing through the center of the discharge tube.
To characterize the trajectory and dilution of the buoyant jet, three cameras were available to record images: two cameras from a stereoscopic PIV system and another camera from a LIF system (LaVision Imager 90 ProX 4M). The stereoscopic PIV system was used as a conventional PIV system, since the system was only focused on the central axis of the jet and the variation in the perpendicular axis to the jet's advance was not swept away. The cameras were synchronized to take simultaneous measures of velocity and concentration.
The CCD array is a 2048 × 2048 pixel sensor consisting of square detectors 7.4 × 7.4 µm 2 with physical dimensions of 15.15 × 15.15 mm 2 . Each CCD camera was mounted with a Nikon AF Nikkor 50 mm 1:1.8D objective with an aperture varying between 1.8 and 22. The CCD cameras were placed parallel to each other and approximately perpendicular to the laser sheet and the jet centerline. Both cameras had the same horizontal and vertical field-of-view (FOV) of 71 cm with a total overlap of both FOV. One camera took the image at time t and the other camera at t + ∆t. The velocity composition was carried out by an inter-correlation algorithm [42], which divides the two successive images of the recorded particles at time t and t + ∆t into small interrogation areas. The position of the correlation peak between the two images corresponds to the most probable displacement of the particles for each of these interrogation areas [42]. The effects of geometric distortion, resulting from small angles in camera alignment, were automatically corrected during image processing and post-processing using LaVision DaVis software. The cameras were moved forward or backward depending on the test configurations in order to zoom in or zoom out on the flow area. Firstly, the cameras were relocated to define their results in the plane that marks the central axis of the jet advance. Secondly, the cameras were placed in such a way that the plane was wide enough to contain from the start of the jet until it was in the so-called lateral dispersion dominant zone, taking into account the number of pixels available for the camera to define accurately the jet evolution.
The LIF camera has a band-pass filter, which allows only the fluorescent light emitted from the tracer dye excited by the laser to pass through. The filter ensures an excellent blocking efficiency of the excitation wavelength with a steep edge at 540 nm providing maximum transmission in its working range (545-800 nm). The measured images were realized with a total window size of 71 × 71 cm 2 and an acquisition rate of 5 Hz. For the PIV measurements, an interrogation window size of 32 × 32 pixel was used. The separation time (∆t) was 20 µs for the jet centerline path around the nozzle and 7 ms for the spreading layer where the jet presents a marked upward trajectory, according to [43].
Polyamide particles with a 50 µm diameter and density of 1030 kg/m 3 were used as seeding to follow the flow velocity variations. The standard cyclic Fourier Transform (FT) correlation function was applied to eliminate the PIV-generated noise and detect the peak that shows the velocity for each time step. For the LIF measurements, the exposure time was 20 ms. The absence of surfactants in the water-free surface was checked by liquid chromatography, since small amounts of them may alter the turbulence phenomena [44]. As passive tracer, Rhodamine WT (C 29 H 29 ClN 2 Na 2 O 5 ) with a concentration of 100 µg/L was used. This concentration was selected to assure the right measurement in the highest dilution zones and to avoid any significant laser attenuation caused by self-absorption of the Rhodamine [39,45]. The data acquisition system presents uncertainties of around 10 −4 m/s for velocities and 1 µg/L for concentrations.
The relation between laser counts and concentration (see Equation (1)) and the relation between thermostat temperature and nozzle temperature (see Equation (2)) were calibrated using a thermometer (Pt-100 Testo 176 T2). In addition, the conductivity was measured by a Crison CM 35+ conductivity probe and translated into salinity by using the Reference Salinity transformation according to the UNESCO thermodynamic equation of seawater [46].
Rhodamine concentration (µg/L) = 0.0176 × Light intensity (counts) Nozzle temperature ( • C) = 0.9635 × Thermostat temperature ( • C) Water 2020, 12, 1324 6 of 22 The characteristics of the four physical tests (V1 to V4) are summarized in Table 1. S a , T a , and ρ a are the ambient salinity, temperature and density, respectively. S o , T o , ρ o , and Q o are the effluent salinity, temperature, density and flow, respectively. Re, Fr, and H are the Reynolds number (see in Equation (3)), the Froude number (see in Equation (4)), and the distance from the jet nozzle axis to the surface of every test, respectively.
where U o is the jet velocity, D the nozzle diameter, g the gravity acceleration, and ∆ρ = ρ a − ρ o the difference between the ambient and effluent density. Subscript "a" refers to the receiving tank and subscript "o" to the discharge effluent.
For every physical test, 1000 images were approximately taken. The images corresponding to the time necessary to achieve the steady state for the mixing regions were eliminated. According to [35][36][37], physical tests reached the steady state when the behavior of the buoyant jet showed variations less than the 5%, being 400 images in our physical tests. PIV measurements were processed calculating the time-averaged velocity (U) and the time-averaged turbulent kinetic energy (TKE), nondimensionalized with the inlet velocity at the nozzle (U MAX ) and the TKE at the nozzle (TKE MAX ), respectively. LIF measurements were processed calculating the time-averaged concentration (C) and the time-averaged fluctuation of the concentration (C'), nondimensionalized with the maximum concentration, i.e., the initial concentration (C MAX ) discharged through the nozzle.

Description of Numerical Models
The numerical simulations were carried out by two semiempirical models (CORJET and VISJET) and one CFD model (OpenFOAM). All numerical simulations were performed by the Neptuno supercomputer located at the IHCantabria facilities, featuring 80 nodes (64 GB RAM per node) of 16 cores (2x Intel Xeon E5-2670).
Both semiempirical models consider three types of forces acting on the jet: (1) an intrusive force in the direction of the longitudinal advance of the jet, (2) a buoyant force in the vertical direction, and (3) a drag force perpendicular to the jet due to the normal component to the path. It is important to note that both semiempirical models cannot simulate the interaction of jets with solid contours or the water-free surface, ending their calculation when the trajectory reaches any of them.
CORJET is an eulerian model that allows the characterization of the detailed geometry of the jets, both with positive and negative buoyancy. For this purpose, it solves the three-dimensional equations integrated in the jet, obtaining the trajectory and the dilution along the axis, considering a Gaussian distribution of the different magnitudes of the flow (speed, concentration, etc.). Within this hypothesis, CORJET considers that three types of forces act on the jet: intrusion force in the X direction due to the current, buoyancy force in the Z direction and dragging force perpendicular to the jet due to the normal component to the jet trajectory in each of the planes that cut it. The CORJET model is not able to represent the jet evolution once it reaches the contours of the system (at the bottom or on top). Due to this restriction, the model does not consider the jet evolution in the lateral dispersion dominant zone. More details of the theory and computation can be found in [47][48][49].
The VISJET program is a flow visualization tool that allows the description of the evolution and interaction of single-point discharges with different angles in a fluid environment. The Lagrangian calculation model of VISJET is called JETLAG. This model, which has been contrasted with experimental and field data, is capable of predicting the initial mixing of discharges in a fluid environment. The model estimates the evolution of the average properties of the discharge in each time step, by integrating the vertical and horizontal momentum equations, the conservation of mass equation and the conservation of heat equation. The vortices at the mass input to the jet are accurately determined, while the drag pressure is ignored. VISJET calculation requires the definition of the receiving environment (environmental parameters), the single-point characteristics and the nature of the discharge. The model takes into account the domain as a regular parallelepiped and once the jet reaches any part of the domain (e.g., the water-free surface) the simulation concludes. More details of the theory and computation can be found in [12].
OpenFOAM is a free and open source CFD model, based on finite volume discretization, to solve complex fluid mechanics problems. OpenFOAM is released and developed mainly by OpenFOAM Ltd. since 2004, being later distributed by the OpenFOAM Foundation. It has a wide user base in most areas of engineering and science, both in commercial and academic organizations. OpenFOAM is organized in a set of C++ modules that make it possible to solve problems from complex fluid flows, involving chemical reactions, turbulence, heat transfer, acoustics, solid mechanics and electromagnetics. It includes tools for meshing in and around complex geometries (e.g., a vehicle), and for processing and visualizing data. A detailed description of the OpenFOAM with the Reynolds equation, as used in this study, can be found in [50].

Setup of the Computational Fluid Dynamics Model
The 3D mesh grid was nondimensionalized by the nozzle diameter with a minimum element size of D/5 near the nozzle and increasing the element size a 10% with a cross-flow size-function until a maximum element size up to D (Figure 2a,b). For the longitudinal advance of the buoyant jet, a constant discretization of D/2 was taken ( Figure 2b). This discretization was chosen based on the biunivocal relationship between the mesh discretization and the jet centerline position and dilution to simulate a discharge by means of CFD [32]. The relationship reported that the required grid size close to the single-port has to be at least 1/5 the nozzle diameter to simulate properly the jet entrance. Moreover, the size of the grid elements can be increased using a size function to achieve a maximum element size of 1 diameter close to the contours of the domain. The numerical discretization was conducted by the SIMPLE algorithm (Figure 2c).
In order to assess the model sensitivity to the mesh grid size, three mesh grids were used: (1) coarse grid with 3 × 10 5 elements, (2) intermediate grid with 6 × 10 5 elements, and (3) fine grid with 2 × 10 6 elements. Independently of the simulation case and/or the mesh grid, the numerical discretization was a Crank-Nicholson second order scheme in time and a Van Leer second order scheme for the advection terms. In order to ensure the numerical model stability, a variable time step depending on the turbulent flow was used. The time step started with 10 −5 s at the beginning of any simulation and finished with 5 × 10 −4 s at the end of any simulation. The simulation time was 200 s, recording model results every 10 −2 s for every numerical test.
The water-free surface (upper boundary) was modelled as a surface with zero shear stress because the variation of water level in the tank was negligible during the duration of the tests. The bottom, wall and nozzle were defined as non-slip closed boundaries. The boundaries in front and parallel to the nozzle were considered as open boundaries under hydrostatic pressure in order to comply with the continuity equation and take out the exceedance mass of the domain. Finally, the initial condition for

Setup of the Semiempirical Models
The semiempirical modelling (CORJET and VISJET) was carried out with the same geometry of the CFD modelling (see Figure 2) and taking into account the same physical conditions for both the discharge and the receiving water body. In addition, it should be noted that these models do not interact with the boundaries, so the simulations conclude when the buoyant jet hits the water-free surface.

Performance Metrics of Numerical Modelling
To evaluate the models' performance, firstly, the coefficient of determination (R 2 ) was calculated as the difference between modelled results and observed values, as expressed in Equation (5). Secondly, the error between both series was assessed using the root-mean square error (RMSE) and the normalized RMSE (NRMSE), displayed in Equations (6) and (7), respectively. Finally, the error between both series was also calculated using the model efficiency (CE), developed by [51] and showed in Equation (8).
where R i is the i-data of the physical tests, S i is the i-data of the numerical tests, R is the average data of the physical tests, and i is the ith value from 1 to N data of the physical tests. CE ranges between −∞ and 1.0 (1.0 inclusive), with CE = 1 being the optimal value. Values between 0.0 and 1.0 are generally viewed as acceptable levels of performance, whereas values <0.0 indicate that the mean observed value is a better predictor than the simulated value, which indicates unacceptable performance. Depending on the CE value, the comparison is considered acceptable (poor) if CE < 0.4, acceptable (-) if 0.4 ≤ CE < 0.6, acceptable (convenient or good) if 0.6 ≤ CE < 0.8, and acceptable (excellent) if CE ≥ 0.8. CE was employed for two major reasons: (1) it is very commonly used in river and estuary model applications [52][53][54][55][56], and (2) [57] CE is also found to be the best objective function for reflecting the overall fit of a model output.

Zonation of Positively Buoyant Jets
As it was mentioned, the evolution of positively buoyant jets is generally separated into four zones (see Figure 1b,c) according to the competitive relation between momentum and buoyancy [11]. In this work, we proposed a new method to delimitate these zones using the angle (α) shaped by the tangent of the centerline trajectory and the longitudinal axis or X-axis (see Figure 3). α can be seen as an indicator of the relation between momentum and buoyancy.
noted that 45 degrees marks the point where the forces in the horizontal direction (momentum forces) are smaller than the forces in the vertical direction (buoyancy forces) according to the force decomposition.    When introducing a positively buoyant jet in a fluid environment, its behavior is governed by its momentum forces so the jet centerline is parallel to the single-port projection. Thus, the momentum dominant zone can practically be assumed to be up to α equals to 10 degrees because density-driven mechanisms are weak. When the positively buoyant jet goes beyond this angle, it reaches a transition area, where the momentum forces act with a similar magnitude as buoyancy forces. Next, if the angle reaches a value greater than 45 degrees, the force decomposition indicates that the buoyancy forces are greater than momentum forces and, therefore, the jet is in the buoyancy dominant zone. Lastly, when the jet angle is again less than 10 degrees, the positively buoyant jet is considered to advance only by lateral dispersion, entering in the so-called lateral dispersion dominant zone.
The Limit 1 is defined when α is equal or greater than 10 degrees, i.e., buoyancy begins to be significant. It should be noted that lateral mixing of the jet with the receiving water starts to decrease (shear layer mixing). This effect is an indicator of the transition between the dominance of momentum forces with respect to buoyancy. Next, the Limit 2 is established when α is equal or greater than 45 degrees because it is the turning point where buoyancy is stronger than momentum. It should be noted that 45 degrees marks the point where the forces in the horizontal direction (momentum forces) are smaller than the forces in the vertical direction (buoyancy forces) according to the force decomposition. All physical tests presented an evolution dominated by momentum for x/D ranging between 29 (V2) and 41.5 (V1). The transition zone was located between 29 ≤ x/D < 65 for V2 and 41.5 ≤ x/D < 90 for V1. Moreover, the buoyant jet evolution was dominated by buoyancy in the region 65 ≤ x/D < 91 for V2 and 90 ≤ x/D < 104 for V1. Independently of the physical test, the buoyant jet evolution shows, for x/D ≥ 104, a movement parallel to the free surface, i.e., entering in the lateral dispersion dominant zone. From these results, we can highlight (1) the higher the flow discharge, the greater the extension of Zones 1 and 2; (2) the higher the density gradient between the discharge and the receiving water body, the greater the extension of Zone 3; and (3) the higher the water depth of the receiving water body, the greater the extension on all zones. the centerline trajectory and the longitudinal axis (X-axis).   All physical tests presented an evolution dominated by momentum for x/D ranging between 29 (V2) and 41.5 (V1). The transition zone was located between 29 ≤ x/D < 65 for V2 and 41.5 ≤ x/D < 90 for V1. Moreover, the buoyant jet evolution was dominated by buoyancy in the region 65 ≤ x/D < 91 for V2 and 90 ≤ x/D < 104 for V1. Independently of the physical test, the buoyant jet evolution shows, for x/D ≥ 104, a movement parallel to the free surface, i.e., entering in the lateral dispersion dominant zone. From these results, we can highlight (1) the higher the flow discharge, the greater the extension of Zones 1 and 2; (2) the higher the density gradient between the discharge and the receiving water body, the greater the extension of Zone 3; and (3) the higher the water depth of the receiving water body, the greater the extension on all zones.    Figure 5 shows, in all cases, the concentration was reduced to 15%-10% of the initial concentration (Co) at Limit 1. Between X/D = 0 and X/D = 20, the test results presented outliers since the rhodamine concentrations used showed a saturation in their emission values and could not provide reliable results in this area. Next, the concentration was reduced to values of the order or less than 5% of Co at the Limit 2. However, in the physical test V2, the development of the buoyant jet generates a higher concentration due to the relationship between the jet momentum and the water depth until it hits the water-free surface. As can be seen in Figure 6, the fluctuation values are less than 1% of CMAX in all the zones of the buoyant jet. Nonetheless, these fluctuations can present values corresponding to 2% of CMAX in the central area of the buoyant jet before reaching the Zone 3, independently of the physical test.

Setup of the Computational Fluid Dynamics Model
In order to compare the results with the three grid resolutions for OpenFOAM, Figure 7 shows a scatter plot of the centerline trajectory (Figure 7a Figure 5 shows, in all cases, the concentration was reduced to 15%-10% of the initial concentration (C o ) at Limit 1. Between X/D = 0 and X/D = 20, the test results presented outliers since the rhodamine concentrations used showed a saturation in their emission values and could not provide reliable results in this area. Next, the concentration was reduced to values of the order or less than 5% of C o at the Limit 2. However, in the physical test V2, the development of the buoyant jet generates a higher concentration due to the relationship between the jet momentum and the water depth until it hits the water-free surface. As can be seen in Figure 6, the fluctuation values are less than 1% of C MAX in all the zones of the buoyant jet. Nonetheless, these fluctuations can present values corresponding to 2% of C MAX in the central area of the buoyant jet before reaching the Zone 3, independently of the physical test.

Setup of the Computational Fluid Dynamics Model
In order to compare the results with the three grid resolutions for OpenFOAM, Figure 7 shows a scatter plot of the centerline trajectory (Figure 7a   As can be seen in Figure 8, the two semiempirical models shown a similar result. It is noteworthy to point out that these models did not need any boundary conditions but an initial condition for the ambient density (calculated with the UNESCO equation from temperature and salinity). From X/D = 20, semiempirical models represented a tendency to decrease in concentration, as seen in the physical tests, but overestimated the dilution by 50% when the models reached the water-free surface and finished their calculation. Regarding OpenFOAM results, it can be seen that CFD simulation provided better results than semiempirical models.
In order to assess the performance of the numerical models to reproduce the physical tests, R 2 , RMSE, NRMSE, and CE were calculated ( Table 2). Figure 9 displays a scatter plot of the centerline trajectory (Figure 9a) in dimensionless positions (X/D and Y/D), and a scatter plot of the C/CMAX along the centerline (Figure 9b) for the physical results versus the numerical results obtained with CORJET (green markers), VISJET (blue markers) and OpenFOAM (red markers) for the momentum dominant zone (square markers), the momentum to buoyancy transition zone (circle markers), the buoyancy dominant zone (triangle markers), and the lateral dispersion dominant zone (diamond markers), respectively. As can be seen in Figure 8, the two semiempirical models shown a similar result. It is noteworthy to point out that these models did not need any boundary conditions but an initial condition for the ambient density (calculated with the UNESCO equation from temperature and salinity). From X/D = 20, semiempirical models represented a tendency to decrease in concentration, as seen in the physical tests, but overestimated the dilution by 50% when the models reached the water-free surface and finished their calculation. Regarding OpenFOAM results, it can be seen that CFD simulation provided better results than semiempirical models.
In order to assess the performance of the numerical models to reproduce the physical tests, R 2 , RMSE, NRMSE, and CE were calculated (Table 2). Figure 9 displays a scatter plot of the centerline trajectory (Figure 9a) in dimensionless positions (X/D and Y/D), and a scatter plot of the C/C MAX along the centerline (Figure 9b) for the physical results versus the numerical results obtained with CORJET (green markers), VISJET (blue markers) and OpenFOAM (red markers) for the momentum dominant zone (square markers), the momentum to buoyancy transition zone (circle markers), the buoyancy dominant zone (triangle markers), and the lateral dispersion dominant zone (diamond markers), respectively. Water 2020, 12, x FOR PEER REVIEW 16 of 23     For both the trajectory and C/C MAX , R 2 values were above 0.8 for OpenFOAM. Moreover, CE values for OpenFOAM were considered as excellent (CE > 0.8) and good (CE > 0.6) for centerline position and dilution along the centerline, respectively. On the other hand, this metric shows semiempirical models provided good results (CE > 0.6) for the centerline position but their results were poor (CE < 0.4) in defining the dilution of a positively buoyant jet.
The three considered models provided good agreement with physical tests for the centerline position. Semiempirical models CORJET and VISJET obtained a CE value of 0.78 and 0.79 (classified as good), respectively, (Table 2) and the CE value for OpenFOAM was classed as excellent, being CE = 0.98. The lower NRMSE was obtained by OpenFOAM with a value of 4.62%, better than the approximately 12% obtained by semiempirical models. In terms of R 2 , the three models obtained a very high value (>0.87), being especially high for OpenFOAM (0.98).
The analysis of the dilution shown the best results were obtained by OpenFOAM, with a CE of 0.64, a NRMSE of 12% and a R 2 greater than 0.8. However, in the case of semiempirical models, the CORJET obtained a CE of 0.45 (acceptable result) and the VISJET of 0.37 (poor result), a NRMSE of 16% and 20%, respectively, and a R 2 less than 0.8 for both semiempirical models.

Zonation of Positively Buoyant Jets
The zone limits of the positively buoyant jet for the four physical tests (V1 to V4) were determined from the experimental, OpenFOAM, CORJET, and VISJET data and summarized in Table 3. In order to assess the performance of the numerical models to define the limits of positively buoyant jet zones, the RMSE, and NRMSE were calculated in Table 4. It should be noted Limits 1, 2 and 3 are defined by α value of 10 degrees (limit of separation between the momentum and transition zone), 45 degrees (limit between the transition zone and the buoyancy dominant zone) and again 10 degrees (limit of the zone where the flow pattern is dominated by lateral dispersion), respectively.  According to the experimental data, Limit 2 was reached at a distance 2.14 times the Limit 1. In the case of numerical models, this distance was 2.15, 2.43 and 2.51 times the Limit 1 for OpenFOAM, CORJET and VISJET models, respectively. Next, Limit 3 was reached at 1.2 times the Limit 2 for both experimental data and the OpenFOAM numerical model. In addition, both the physical tests and the results obtained with OpenFOAM found that the relationship between the Limit 3 and Limit 1 was 2.7. All these values were the average data of the four tests.
As shown in Table 4, the OpenFOAM numerical model best represented the evolution of the positively buoyant jet in all the proposed zones until it reached the lateral dispersion dominant zone (Zone 4). In the Limits 1 and 2, the OpenFOAM model had an RMSE approximately five times lower than the semiempirical models CORJET and VISJET. In the case of Limit 3, the OpenFOAM values cannot be compared with the semiempirical models since they cannot give results when the buoyant jet advances to the lateral dispersion dominant zone.
When comparing the zonation Limits calculated by the experimental data, the NRMSE shows that the OpenFOAM error was less than 7% in both Limits 1 and 2. For the case of CORJET, the differences were 39% and 22.6% in Limits 1 and 2, respectively, and 45.6% and 18.2% for the case of VISJET in those limits. As can be seen, the models are capable of giving results with less standardized error with the jet advance up to the moment of impact with the water-free surface.

Physycal Modelling
Physical tests have provided similar results to previous studies along the Zones 1 and 2 of positively buoyant jets [7,8,11,58].
A fundamental aspect of buoyant jets is the relationship between the dispersion of mass and energy (velocity) related to its centerline. This relationship (λ) is defined throughout the ratio between the maximum dimensionless concentrations C Maxi and velocities U Maxi in a selected profile (see Equation (9)). According to [59], the value of λ follows a non-linear pattern in the centerline advance along the longitudinal axis.
For instance, λ was obtained for the profiles X/D = 20, X/D = 60, X/D = 80 and X/D = 110, which are representative of the four different zones in the evolution of positively buoyant jets, showing an average value for the four physical tests of 0.73, 0.57, 0.51 and 0.60, respectively. The highest values of λ were located inside the momentum dominant zone (Zone 1) because this zone includes the most energetic vortices. λ values started to decrease in the transition zone (Zone 2) due to the reduction of the turbulent energy. In the Zone 3, λ continued decreasing because the turbulent energy is weaker. In the Zone 4, λ values increased due to the vortices destruction and the advance by dispersion. Finally, λ reached a constant value when the turbulent wake was fully developed in the lateral dispersion dominant zone (Zone 4).
From the Zone 3 until Zone 4, the positively buoyant jet was not axisymmetric because the water-free surface behaved like a wall. Therefore, the buoyant jet only presents dispersion in the half bottom of its volume. This fact reduces the differences between energy and mass dispersion (half dispersion of mass and energy) so the net values of both were more similar, considering that λ is less than one in all profiles that cut the advance of the buoyant jet. These results agree with [59] as the higher stretching vortex processes were in the first part of the buoyant jet evolution (increasing mixing processes), then they diminished in the impact zone with the water-free surface (dominated by advection) and then lightly increased with the turbulent wake.

Numerical Modelling
As shown in Figure 8, semiempirical models do not allow calculating the jet trajectory and dilution in the Zone 4. Due to their non-interaction with the water-free surface, semiempirical models do not change their trajectory when approaching it. As a consequence, semiempirical models presented a shorter distance of the impact point with the water-free surface than the experimental data or the CFD model. On the contrary, CFD models provide results in a continuous mode from the Zone 1 until the Zone 4.
Numerical models enable accurate results with less standardized error for the evolution of positively buoyant jets until the Limit 2. In the case of the OpenFOAM model, the only one capable of obtaining Limit 3, it can be seen that Limit 3 has a standard error that doubles the value obtained in the previous limits. The increase of the error in this area is due to the buoyant jet mechanics is governed by diffusion and, therefore, it is more sensitive to numerical diffusion phenomena, which were less appreciable in the areas where the buoyant jet moved either by momentum or by buoyancy gradients.
Overall, the obtained results suggest that semiempirical models can provide results in a range error of ±50% and CFD models obtain much better results for both centerline position and dilution in the four zones of a positively buoyant jet. Taking into account the accuracy of numerical models in every zone, model results suggest the use of semiempirical models for only calculations in Zones 1 and 2 while CFD models are applicable for anyone. However, depending on the case study, the high computational costs of CFD models should be always considered for the model selection due to the trade-off between computational time and model accuracy. Accordingly, CFD models could be used in further work to analyze the zonation of other buoyant jet configurations without laboratory testing, meanwhile, semiempirical models could not be used as additional information generators due to the displayed errors.

Zonation of Positively Buoyant Jets
The buoyant jets have been delimited by length scales based on physical modelling [11,[35][36][37]. However, in many cases, these length scales cannot be calculated with numerical modelling. This work tries to overcome these constraints by proposing a simple indicator, α, able to delimitate the evolution of positively buoyant jets from physical and numerical model results. Thus, classical methods based on the dimensional analysis of variables may prescribe uncertainty about the location of the different zones because the classification use length scales build through linear fits with an error higher than 30% [38].
The Limit 1 could correspond to the momentum length scale (L M ) proposed by [11] because both mark the limit of momentum dominant zone-based characteristics (single-port diameter, momentum and buoyancy). The obtained values were very similar to those obtained with the proposed methodology. Regarding the Limit 2, the classical length scales do not include a specific methodology to obtain it. For instance, the calculation of the buoyancy length scale (L B ), proposed by [11], is diffused and constrained to the existence of stratification in the receiving medium. Finally, according to the Roberts experiments [35][36][37], the Limit 3 is obtained when the vertical concentration profile in the layer where the buoyant jet is progressing presents longitudinal variations lower than 5%. This calculation was in line with the lengths obtained by the proposed method for the Limit 3. Table 5 shows the values of these scales for the laboratory tests.
As can be seen in Table 5, the values of Limit 1 obtained by L M are closer to the semiempirical models results than the laboratory test data, displaying an error of 44% compared with physical test data and displaying more error than modelling errors obtained by means of OpenFOAM (error value of 6.4% compared with test data). In the case of Limit 2, L B could not be calculated because the environment did not present stratification. Finally, the point marked by [35][36][37], as Limit 3, shows a variation of 25% compared with the physical tests; this variation is greater than the variation for OpenFOAM (a variation for the Limit 3 location of 11.9%).

Conclusions
A new method to delimitate the four characteristic zones of positively buoyant jets interacting with the water-free surface was proposed using the angle (α) shaped by the tangent of the centerline trajectory and the longitudinal axis. These zones are (1) momentum dominant zone (jet-like), (2) momentum to buoyancy transition zone (jet to plume), (3) buoyancy dominant zone (plume-like), and (4) lateral dispersion dominant zone. The validation of the proposed method was tested by the results obtained with physical and numerical modelling. The proposed method could be imposed as a more general strategy when analyzing the different zones in which the evolution of positively buoyant jets can be separated because it can be always calculated by physical and numerical model results.
The proposed parameter α can be seen as an indicator of the relation between momentum and buoyancy. The Limit 1 (Zone 1 to 2) is defined when α is equal or greater than 10 degrees, i.e., buoyancy begins to be significant. Next, the Limit 2 (Zone 2 to 3) is established when α is equal or greater than 45 degrees because this is the turning point where buoyancy becomes stronger than momentum. Finally, the Limit 3 (Zone 3 to 4) is specified when α is again equal or lower than 10 degrees, i.e., the buoyant jet evolution is again parallel to the longitudinal axis (X-axis) and the water-free surface.
Furthermore, the evolution of positively buoyant jets was studied with non-intrusive techniques (PIV and LIF) by analyzing four physical tests in the four characteristic zones. The coupled analysis of PIV and LIF shows that the dispersion of mass (concentrations) is always greater than the dispersion of energy (velocity) during the evolution of positively buoyant jets.
Regarding numerical modelling, on the one hand, the semiempirical models (CORJET and VISJET) for the study of positively buoyant jets underestimate the advance of the jets up to their impact with the water-free surface, where they end the calculation. In turn, these numerical models overestimate the actual dilution reached by positively buoyant jets, especially close to the impact zone. On the other hand, the CFD model (OpenFOAM) for the study of positively buoyant jets was able to faithfully reproduce the behavior of the buoyant jet for all the proposed zones according to the physical results. CFD can also avoid the simplifications made by semiempirical models for the study of jet dynamics (trajectory and dilution). In addition, CFD models are able to consider all such complexity in relation to the geometry of the diffusers and the receiving water body. The use of this technique allows, in a single calculation, results to be obtained in all the zones.
The application of the proposed method by using numerical modelling has shown that the semiempirical models underestimate the distance of the limits for the different zones proposed in the methodology. The comparison of the laboratory test data with the OpenFOAM results showed an error of 6.4% and 5.4% in the definition of Limits 1 and 2, respectively. Moreover, a comparison with semiempirical models showed higher errors, with error values of 39% and 22.6% in the case of CORJET and 45.6% and 18.2% in the case of VISJET for the location of the Limits 1 and 2, respectively. It should be noted that the semiempirical models cannot offer information on Limit 3. In contrast, the OpenFOAM model displayed an error value of 11.9% for Limit 3 with respect to the laboratory tests.
Therefore, CFD models increases the understanding about the evolution of positively buoyant jets in the aquatic medium and will facilitate the improvement of the discharge designing process. This fact could be a starting point for their use as information/data providers, i.e., a virtual laboratory, in order to generate better semiempirical models or advanced statistical models based on artificial intelligence techniques, enabling the analysis of new discharge system configurations without large computational costs or large physical modelling costs. Finally, CFD model results might be used to develop an abacus to gauge the influence of key parameters on the initial dilution, location, width and thickness of positively buoyant jets in order to predict the behavior of any discharge in the near-field and couple this with far-field models.