Veriﬁcation and Validation of URANS Simulations of the Round Buoyant Jet in Counterﬂow

: This paper presents a study on the veriﬁcation and validation (V&V) of numerical solutions for round buoyant jets in counterﬂow. The unsteady ﬂow was simulated using an unsteady Reynolds-averaged Navier–Stokes (URANS) solver with a two-phase mixture model. This work aimed to quantitatively investigate the reliability and applicability of various uncertainty estimators in the simulation of a buoyant jet in counterﬂow. Analysis of the discretization uncertainty estimation results revealed that the factor of safety (FS) and the modiﬁed FS (FS1) methods were the appropriate evaluation estimators in the simulation of a buoyant jet in counterﬂow. Validation by comparison with the experimental data indicated that the area without achieving the validation at the validation level was strongly related to the shear layer between the jet ﬂow and the ambient ﬂuid. Moreover, the predicted concentration contours, coherent structures, and centerline concentration were strongly affected by the grid resolution.


Introduction
The effluent discharge in a receiving water body through outfalls and diffusers is generally in the form of a turbulent jet or plume. The flow characteristics and mixing process of jets are affected by the receiving environment (e.g., still water, density stratification, and crossflow). Numerous experimental and numerical investigations have been conducted on the flow characteristics of jets, and some typical phenomena have been identified [1][2][3][4][5]. The situation in which a jet is injected horizontally to the ambient fluid in an opposite direction is called a jet in counterflow, which is found in wastewater outfalls, particularly a rosette jet group [6]. Some experimental studies investigate the behavior of a jet in counterflow using laser Doppler anemometry (LDA) [7,8], particle imaging velocimetry (PIV), and laser induced fluorescence (LIF) [9]. Traditional studies of counterflowing jets involve velocity and concentration fields [10,11]. From the literature, it can be seen that the studies of non-buoyant jets in counterflow have been studied widely. In ocean outfalls, wastewater is lighter than the ambient seawater, and the discharging wastewater is in the form of a buoyant jet. Therefore, analysis of the mixing process becomes complex if both buoyancy and a jet in counterflow are considered, but a limited number of works have done so [12,13].
The obtained reliable numerical results suggest that numerical solutions are more detailed than experimental data in reflecting the performance of jet characteristics. At present, Reynolds-averaged Navier-Stokes (RANS) and large eddy simulation (LES) are widely used for numerical simulations [14,15]. However, limited attention has been directed toward the verification and validation (V&V) of a buoyant jet in counterflow through RANS or LES methods, although the mesh influence has been explored by a few works [16,17]. V&V is the precondition for evaluating the reliability before numerical solutions are accepted and flow characteristics are analyzed. In the context of CFD (computational fluid dynamics)

Numerical Methods and Uncertainty Estimators
A two-phase mixture model [14,37,38] in the commercial code FLUNET (ANSYS, Canonsburg, PA, USA,) was adopted in this paper. The wastewater jet flow was composed of fluid f and particle s phases, which were separate, and these phases were allowed to interpenetrate where α f and α s are the volume fractions of the fluid and particle phases, respectively. Following the concept of slip velocities, the mixture model allowed the phases to move at different velocities. A Schiller-Naumann drag model was adopted to describe the interaction between the fluid and particle phases.

Governing Equations
The continuity equation for the mixture was: where ρ m was the mixture density, as in: where α i was the volume fraction of phase i, n was the number of phases (n = 2). u m was the mass-averaged velocity, as in: The momentum equation for the mixture was found by summing the individual momentum equations for all phases as: α i ρ i u dr,i u dr,i ) + ρ m g + F, (4) where n was the number of phases (n = 2), F was a body force, µ m was the viscosity of the mixture α i µ i ), and u dr,i was the drift velocity for secondary phase i (u mdr,i = u i − u m ). The drift velocity for the secondary phase u dr,i and the relative velocity u sf were linked by the expression: where α s was the volume fraction of the secondary phase, and u sf was the velocity of the secondary phase (s) relative to the velocity of the primary phase (f ) as given by: where τ s was the particle relaxation time given as τ s = ρ s d 2 s /(18µ f ), and d s was the diameter of the particles of the secondary phase s. The default drag function f drag was taken as: where Re was the relative Reynolds number. The relative Reynolds number for phases f and s was defined as [39]: For the continuity equation for the secondary phase s, the volume fraction equation could be given by: ∂ ∂t (α s ρ s ) + ∇ · (α s ρ s u m ) = −∇ · (α s ρ s u dr,s )

Turbulent Equations for the Mixture
The renormalization group k − ε (RNG k − ε) turbulence model was applied to close the governing equations. The turbulent kinetic energy k and its rate of dissipation ε could be expressed as follows: where µ t,m was the turbulent viscosity (µ t,m = ρ m C µ k 2 /ε), G k,m was the production term of the turbulence kinetic energy (G k,m = µ t,m [∇u m + (∇u m ) T ] : ∇u m ), and C * 1ε was defined as . In these equations, the model constants had the following values:

Geometry and Computational Conditions
Lam et al. conducted an experiment on a round buoyant jet in counterflow [16,17]. Figure 1 shows a computational domain with dimensions of 50 d × 20 d × 40 d in the abscissa (x), ordinate (y), and vertical (z) directions, respectively, where d is the diameter of the jet nozzle. The computational domain was based on the jet diameter, the area affected by the jet, and the computational resources. A round buoyant jet with an initial jet velocity U jet , jet density ρ jet , and jet nozzle diameter d (d = 0.5 cm) was injected into the ambient fluid at an ambient velocity U a and ambient density ρ a (Table 1). According to the referenced experimental data, salt water (ρ a = 1.0355 g/cm 3 ) was used as the ambient fluid, and water (ρ jet = 0.9984 g/cm 3 ) was used as the negatively buoyant jet effluent. The jet-to-current velocity ratio was R = U jet /U a = 7.5, the densimetric Froude number was F = U jet / (∆ρ/ρ a )gd = 5, where ∆ρ was the density difference between the densities of the discharging jet ρ jet and the ambient fluid ρ a [16]. The Reynolds number of the jet nozzle was Re jet = U jet d/υ = 1067.
The renormalization group k−ε (RNG k−ε) turbulence model was applied to close the governing equations. The turbulent kinetic energy k and its rate of dissipation ε could be expressed as follows: , , (11) where μt,m was the turbulent viscosity ( In these equations, the model constants had the following values: 0.0845

Geometry and Computational Conditions
Lam et al. conducted an experiment on a round buoyant jet in counterflow [16,17]. Figure 1 shows a computational domain with dimensions of 50 d × 20 d × 40 d in the abscissa (x), ordinate (y), and vertical (z) directions, respectively, where d is the diameter of the jet nozzle. The computational domain was based on the jet diameter, the area affected by the jet, and the computational resources. A round buoyant jet with an initial jet velocity Ujet, jet density jet  , and jet nozzle diameter d (d = 0.5 cm) was injected into the ambient fluid at an ambient velocity Ua and ambient density a  (Table 1) .
According to the referenced experimental data, salt water

Numerical Setup
In this paper, the simulations were performed using a finite volume method to discretize the governing equations by using a coupled solver. The SIMPLEC scheme was used for the pressure-velocity coupling. The second-order upwind scheme was adopted to discretize the diffusion and convection terms in the governing equations, and the QUICK scheme was utilized for the convection terms in the volume fraction equations.
The velocity inlet boundary condition was set at the jet nozzle (x = 0) and the inlet of the counterflow (x = 25 d). The pressure outlet boundary condition with a static (gauge) pressure (x = −25 d) was set as the outlet of the counterflow. No slip boundary condition was adopted for the left, right, and bottom walls, and the free surface was a frictionless rigid lid, at which a plane of symmetry was applied.
Grid studies were conducted by four grids, which provide two separate grid studies (grids 1-3 and 2-4) for V&V analysis. The grids were created with a uniform grid refinement ratio where h i+1 and h i were the grid spacing between two successively refined grids. The detailed information of the grid series is shown in Table 2. An O-H (O-shaped and H-shaped) grid topology was applied to the four grids. The grid details near the jet nozzle are presented in Figure 2. The near-wall spacing y+ was less than 1 to guarantee that the first node was in the viscous sublayer.

Numerical Setup
In this paper, the simulations were performed using a finite volume method to discretize the governing equations by using a coupled solver. The SIMPLEC scheme was used for the pressure-velocity coupling. The second-order upwind scheme was adopted to discretize the diffusion and convection terms in the governing equations, and the QUICK scheme was utilized for the convection terms in the volume fraction equations.
The velocity inlet boundary condition was set at the jet nozzle (x = 0) and the inlet of the counterflow (x = 25 d). The pressure outlet boundary condition with a static (gauge) pressure (x = −25 d) was set as the outlet of the counterflow. No slip boundary condition was adopted for the left, right, and bottom walls, and the free surface was a frictionless rigid lid, at which a plane of symmetry was applied.
Grid studies were conducted by four grids, which provide two separate grid studies (grids 1-3 and 2-4) for V&V analysis. The grids were created with a uniform grid refinement ratio , where hi+1 and hi were the grid spacing between two successively refined grids.
The detailed information of the grid series is shown in Table 2. An O-H (O-shaped and H-shaped) grid topology was applied to the four grids. The grid details near the jet nozzle are presented in Figure 2. The near-wall spacing y+ was less than 1 to guarantee that the first node was in the viscous sublayer.    In consideration of the computational time and the convergence, 20 iterations per time step with a 10 −4 residual criterion were implemented to balance the accuracy and the efficiency. The solution was iterated within a computational time step of the simulations for 0.0056 d/U a time units. The simulations were initially conducted under steady flow conditions. Then, the unsteady flow, with the steady flow results as the initial conditions, was simulated over approximately 336 d/U a time units.

Uncertainty Estimators
The current uncertainty estimators were based on the Richardson extrapolation method. This study aims to realize discretization uncertainty estimation with truncated power series expansions. The theoretical error δ RE was given by: where S h was the solution with the grid spacing h, S 0 was the exact solution of the discrete equations, α f was a constant, and p f was the formal order of accuracy. The constant grid refinement ratio r was given by: where ∆x 1 , ∆x 2 , and ∆x 3 were the grid spacing with fine, medium, and coarse grids, respectively. If the solutions for the fine, medium, and coarse grids were S 1 , S 2 , and S 3 , respectively, the solution change e, the convergence ratio R, and the observed order of accuracy p k were defined as: R = e 21 /e 32 (15) p k = ln(e 32 /e 21 ) ln(r) The basic estimated errors adopted in the seven estimators were defined as: When the 0 < R < 1 monotonic convergence was achieved, and the Richardson extrapolation method was used to estimate the numerical error, the reliability and applicability of the seven uncertainty estimators, namely, CF, FS, FS1, GCI, GCI−OR, GCI−LN, and GCI−R, were investigated quantitatively. Detailed descriptions for the seven uncertainty estimators are as follows.
(1) Correction factor method (CF) The correction factor method developed by Stern et al. [31] and modified by Wilson et al. [32] is given by: where CF is the correction factor, which is an indicator of the distance from the asymptotic range, as in: (2) Factor of safety method (FS) and its modified version (FS1) The factor of safety method developed by Xing and Stern [33] is given by: where P is an indicator similar to the correction factor that shows the distance from the asymptotic range: FS1 is the modified version of FS proposed by Xing and Stern [34] and is given by: (3) GCI and its modified versions (GCI−OR, GCI−LN, and GCI−R) The GCI developed by Roache [27] is given by: The modified version GCI−OR proposed by Oberkampf and Roy [28] is given by: The modified version GCI−LN proposed by Logan and Nitta [29] is given by: The modified version GCI−R proposed by Roache [30] is given by: where P GCI−R = min(p k , p f ). Among these seven uncertainty estimators, CF and GCI have been widely used in various fields in the literature [40,41]. Xing and Stern [33,34] pointed out the drawbacks of CF and GCI and proposed FS and FS1 to overcome them. However, the applicability and reliability of these seven uncertainty estimators for the round buoyant jet in counterflow have yet to be investigated. The analyses and results of the analysis are as follows.

Verification and Validation Procedures
Uncertainty assessment was performed using the V&V method following the quantitative procedures proposed by Stern et al. [42] with an improved factor of safety [33]. It is worth noting that uncertainties were defined as deficiencies in the stages or results of the modeling process that were due to the lack of knowledge. On the other hand, errors are recognizable deficiencies in any phase of modeling and simulation that are not due to lack of knowledge [31]. In this paper, simulation numerical uncertainty U SN was expressed in terms of graphical methods for the iterative uncertainty U I and a generalized Richardson extrapolation for the grid uncertainty U G as follows: The validation uncertainty U V considers the simulation numerical uncertainty U SN and the experimental uncertainty U D , as in: The comparison error |E| is defined by the absolute value of the difference between the experimental data D and the simulation values S as |E| = |D − S|. If the comparison error |E| is within ±U V , the solutions are validated at the level of U V , and vice versa.

Statistic Convergence Analysis
Two quantities, namely, stream-wise velocity and concentration, were extracted from the time history to study the statistical convergence. The monitoring points near the jet flow are shown in Figure 3. Points A1, B1, and C1 were adopted to investigate the statistical convergence. Figure 4 shows a portion of the iteration history on grid 1 for the normalized stream-wise velocity u/U a and the normalized concentration c/c max at points A1, B1, and C1. This portion reflected a calculation run from a previous solution and not the total iterative history. In Figure 4a,c, the variations of the stream-wise velocity and concentration were relatively small. In Figure 4b,c, the magnified views of the stream-wise velocity and concentration at point C1 showed that the variation in flow quantities was less than 0.01%S, where S was the simulation result of the quantities. The level of iterative uncertainty for grid 1 was at least two orders of magnitude lower than the corresponding grid uncertainty, suggesting that the iteration uncertainties on grid 1 could be ignored. In addition, the iterative uncertainties calculated for grids 2-4 from their histories were less than 0.01%S, implying that for all four grids, the iteration uncertainties were negligible, and the variation of the quantities exerted an insignificant effect on the estimated uncertainty.
within ±UV, the solutions are validated at the level of UV, and vice versa.

Statistic Convergence Analysis
Two quantities, namely, stream-wise velocity and concentration, were extracted from the time history to study the statistical convergence. The monitoring points near the jet flow are shown in Figure 3. Points A1, B1, and C1 were adopted to investigate the statistical convergence. Figure 4 shows a portion of the iteration history on grid 1 for the normalized stream-wise velocity u/Ua and the normalized concentration c/cmax at points A1, B1, and C1. This portion reflected a calculation run from a previous solution and not the total iterative history. In Figure 4a,c, the variations of the stream-wise velocity and concentration were relatively small. In Figure 4b,c, the magnified views of the stream-wise velocity and concentration at point C1 showed that the variation in flow quantities was less than 0.01%S, where S was the simulation result of the quantities. The level of iterative uncertainty for grid 1 was at least two orders of magnitude lower than the corresponding grid uncertainty, suggesting that the iteration uncertainties on grid 1 could be ignored. In addition, the iterative uncertainties calculated for grids 2-4 from their histories were less than 0.01%S, implying that for all four grids, the iteration uncertainties were negligible, and the variation of the quantities exerted an insignificant effect on the estimated uncertainty.    Figure 5 presents the sensitivity of the normalized stream-wise velocity and the normalized concentration results to the variation of the refinement ratios at points A1, B1, and C1. A fitted curve by least squares approach to the simulated results at the monitoring points and a curve fitted with  Figure 5 presents the sensitivity of the normalized stream-wise velocity and the normalized concentration results to the variation of the refinement ratios at points A1, B1, and C1. A fitted curve by least squares approach to the simulated results at the monitoring points and a curve fitted with the theoretical order of accuracy p = 2 are also plotted in Figure 5. As shown, the quantities with different refinement ratios were scattered around the smooth curves. As the refinement ratio approached 1, i.e., a finer grid, the quantities became more closely distributed. This phenomenon qualitatively illustrated that grid refinement played a role in obtaining more accurate solutions. The fitted curves at monitoring points B1 and C1 were closer to the theoretical curves with the theoretical order of accuracy p = 2, implying that the numerical solution converged as expected. By contrast, the fitted curve at monitoring point A1 did not fit the theoretical curve well. For instance, the fitted order of accuracy of concentration p = 3.37 was slightly larger than the theoretical order of accuracy. The large fitted order of accuracy of concentration indicated that the numerical calculation converged faster than expected but necessarily with a higher accuracy. The poor convergence at point A1 could be due to the flow pattern at point A1. Point A1 is located at the shear layer between the jet flow and the entrained counterflow. If the flow pattern corresponding to the entrainment and shear-layer vortices are complex, the accuracy of the numerical simulation results may be diminished. Table 3 lists the uncertainty estimation of the seven estimators for the concentrations at the monitoring points shown in Figure 3. The uncertainties between the two sets of grid analysis for grids 1-3 and 2-4 were first compared. The uncertainties of grids 2-4 at the monitoring points were mostly larger than those of grids 1-3. In particular, the uncertainties at points A3, B1, B3, C1, and C3 on grids 2-4 were relatively large due to the coarse grids. This finding indicated that as the grids became refined, the numerical solutions converged more. The hydraulic and turbulent performances varied depending on the location of the monitoring points. The uncertainties at the jet boundary, such as points A3, B1, B3, C1, and C3, were larger than that along the jet centerline (i.e., points A2, B2, and C2). This was due to strong shearing and entrainment existing in the jet boundary, which increased the difficulty of obtaining the variables at this location. Additionally, the appropriate uncertainty estimator among the seven uncertainty estimators, namely, CF, FS, FS1, GCI, GCI−OR, GCI−LN, and GCI−R, was selected. In general, the P (P = p k /p f ) values varied considerably, which is common in practical applications [33]. As the P values approached zero, the uncertainties measured by GCI−OR were approximately 86% smaller than that by FS and FS1. When the P values were between 0.5 and 0.9, the uncertainties using GCI and its three modified versions were either relatively large or relatively small. When the P values were close to 1, i.e., P = 1.00 at point A3 on grids 1-3, the uncertainties using all other methods were balanced except for the small uncertainties using the CF method. When the P values were much larger than 1, i.e., P = 1.76 at point B3 on grids 1-3, the uncertainties using the CF method are were small, followed by GCI, GCI−OR, GCI−LN, and GCI−R. Thus, comparison of the uncertainties under seven estimators at the monitoring points revealed that the CF method, along with GCI and its modified versions, were not hindered by the drawback of a relatively small estimated uncertainty pointed out by Xing and Stern [33,34]. Therefore, the uncertainties using FS and FS1 were more reliable, followed by CF and GCI−R.  Table 3 lists the uncertainty estimation of the seven estimators for the concentrations at the monitoring points shown in Figure 3. The uncertainties between the two sets of grid analysis for grids 1-3 and 2-4 were first compared. The uncertainties of grids 2-4 at the monitoring points were mostly larger than those of grids 1-3. In particular, the uncertainties at points A3, B1, B3, C1, and C3 on grids 2-4 were relatively large due to the coarse grids. This finding indicated that as the grids became refined, the numerical solutions converged more. The hydraulic and turbulent performances varied depending on the location of the monitoring points. The uncertainties at the jet boundary, such as points A3, B1, B3, C1, and C3, were larger than that along the jet centerline (i.e., points A2, B2, and C2). This was due to strong shearing and entrainment existing in the jet boundary, which increased the difficulty of obtaining the variables at this location. Additionally, the appropriate uncertainty estimator among the seven uncertainty estimators, namely, CF, FS, FS1, GCI, GCI−OR, GCI−LN, and GCI−R, was selected. In general, the P (P = pk/pf) values varied  Table 4 presents the uncertainty estimation of seven estimators for the stream-wise velocity at monitoring points. The trend of the uncertainty comparison for the stream-wise velocity between the two sets of grid analysis was similar to the trend of the concentration above. Although the P varied considerably, they did not approach zero, as shown in Table 4. When 0 < P < 0.5, the uncertainties using GCI and GCI−LN were 38% smaller than those using FS and FS1, which was slightly different from the regularity in Table 3. Although slight differences were noted in the uncertainty comparison of the seven estimators, FS and FS1 were still identified to be more suitable for evaluating the uncertainty for a buoyant jet in counterflow. Nonetheless, comparing the magnitude of uncertainties could not determine the optimal solution but only measure the influence of grid refinement. Therefore, the numerical solution should be validated, as will be discussed in the subsequent section.

Relationship between the Buoyant Jet in Counterflow and the Uncertainty Estimation
The validation results for the concentration at the monitoring points on grids 1-3 and 2-4 are given in Tables 5 and 6. In these validation procedures, U D is computed for 5% of the experimental data as the error analysis in accordance with the literature [16,17]. Table 5 lists the simulation uncertainty U SN , the validation uncertainty U V and the comparison error |E| at the monitoring points. Points A1, B3, C2, and C3 evaluated using all uncertainty estimators denoted the areas without achieving the validation at the U V level, indicating that modelling errors existed in the simulations and that the model must be improved by incorporating |E|. Point A2 using all uncertainty estimators was the area with validation at the U V level. Points A3, B1, B2, and C1 by CF, GCI, GCI−OR, GCI−LN and GCI−R were the areas without validation, whereas these points using FS and FS1 were the areas with validation. Therefore, the reliability of FS and FS1 for evaluating the uncertainty of the buoyant jet in counterflow can be proved further in the validation section.
In Table 6, points A2, A3, and B1 on grids 2-4 were located at the area with validation at the U V level using all uncertainty estimators, and validation was achieved at points C1 and C2 using FS and FS1. Compared with the validation results in Table 5, the area with validation became smaller for a lower grid resolution, indicating that grid refinement improved the accuracy of the numerical solutions. After the reliability of the seven uncertainty estimators was analyzed, the relationship between the area without validation and the buoyant jet in counterflow was examined. The results showed that the areas without validation at the U V level were mostly distributed in the shear layer between the jet flow and the entrained counterflow (Figure 6), such as for points A1, B3, and C3. The areas near the jet centerline were the regions with validation, such as points A2 and B2. This phenomenon may be related to the complex entrainment of the jet shear layer. Complex entrainment refers to the scenario when the jet in the near-field begins to bend, roll up, and shed, and the vortices that developed along the jet boundary (external and internal shear-layer vortices) are strongly sheared as they propagate along the jet flow. A possible reason for the large error at the shear layer was that the strong shearing prevented a relatively stable and accurate value from being obtained due to some experimental measurement error and modelling errors caused by mesh skewness, the shortcoming of the turbulent model, the simplified boundary condition, and so on.

Evaluation of Grid Resolution and Flow Analyses
After the quantitative analyses of the grid independence, the intuitive results are presented as follows. Figure 7 shows the comparison of the experimental and predicted normalized concentration c/cmax in the plane of y = 0 on grids 1-4. As shown in Figure 7, the predicted concentration contours agreed well with the experimental LIF image when the grid resolution was improved. It matched well in the near-field of the jet flow, and a considerable deviation existed in the boundary and far-field of the jet flow even at a fine grid resolution.

Evaluation of Grid Resolution and Flow Analyses
After the quantitative analyses of the grid independence, the intuitive results are presented as follows. Figure 7 shows the comparison of the experimental and predicted normalized concentration c/c max in the plane of y = 0 on grids 1-4. As shown in Figure 7, the predicted concentration contours agreed well with the experimental LIF image when the grid resolution was improved. It matched well in the near-field of the jet flow, and a considerable deviation existed in the boundary and far-field of the jet flow even at a fine grid resolution.

Evaluation of Grid Resolution and Flow Analyses
After the quantitative analyses of the grid independence, the intuitive results are presented as follows. Figure 7 shows the comparison of the experimental and predicted normalized concentration c/cmax in the plane of y = 0 on grids 1-4. As shown in Figure 7, the predicted concentration contours agreed well with the experimental LIF image when the grid resolution was improved. It matched well in the near-field of the jet flow, and a considerable deviation existed in the boundary and far-field of the jet flow even at a fine grid resolution.  In general, typical vortices were present in the jet flow, such as shear layer vortices, counter-rotating vortex pair, horseshoe vortices, and wake vortices. Figure 8 presents the coherent structures of the buoyant jet in counterflow for the four grids based on the Q-criterion [43]. The coherent structures became more complex as the grid resolution improved from coarse to fine. The results for grid 1 captured the shear-layer vortices to show the entrainment between the jet flow and the ambient fluid, whereas the results for grids 2-4 obtained the counter-rotating vortex pair. The counter-rotating vortex pair is a distinct vortex structure in the mean flow field but not in the instantaneous flow field [17,44], implying that grids 2-4 were too coarse to reflect the flow characteristics. This confirmed the point that the predicted coherent structures could reflect the strong shearing between the jet flow and the ambient fluid only at a high level of grid resolution. In general, typical vortices were present in the jet flow, such as shear layer vortices, counter-rotating vortex pair, horseshoe vortices, and wake vortices. Figure 8 presents the coherent structures of the buoyant jet in counterflow for the four grids based on the Q-criterion [43]. The coherent structures became more complex as the grid resolution improved from coarse to fine. The results for grid 1 captured the shear-layer vortices to show the entrainment between the jet flow and the ambient fluid, whereas the results for grids 2-4 obtained the counter-rotating vortex pair. The counter-rotating vortex pair is a distinct vortex structure in the mean flow field but not in the instantaneous flow field [17,44], implying that grids 2-4 were too coarse to reflect the flow characteristics. This confirmed the point that the predicted coherent structures could reflect the strong shearing between the jet flow and the ambient fluid only at a high level of grid resolution. The centerline concentration based on the longitudinal setting of the axial coordinate system is depicted in Figure 9. As shown, radial lines perpendicular to the jet centerline were cut along the jet flow. ζ and η represent the axial and radial directions in the jet-trajectory-based coordinate system. η was assumed positive inside the jet flow, and xp, zp, and ζp were the abscissa, vertical, and axial coordinates of the penetration point, respectively. In Figure 9b, the concentration dropped after the jet was injected from the jet nozzle, decayed at a rapid rate, and then slowed down. The predicted centerline concentrations along the jet trajectory agreed well when the grid resolution became finer. On the basis of the uncertainty analyses, concentration contours, coherent structures, and centerline concentration, as well as the relationship between the numerical accuracy and computational resources, the buoyant jet in counterflow was appropriate on grid 1 for F = 5 and R = 7.5. The radial concentration profile on grid 1 is explored in Figure 10 to explore the decay of the jet concentration in detail. In the near-field of the jet, such as ζ/d = 1.1 and 3.1, the radial concentration decayed rapidly, both inside and outside the jet flow. Along the jet trajectory, the radial concentration inside the jet decayed slower than that outside. The probable reason is that the jet flow outside was facing the counterflow and strongly sheared, whereas the jet flow inside was forced to roll up and shed. Therefore, the area of the jet flow outside boundary was more difficult to validate, which was consistent with the results of the uncertainty analyses. The centerline concentration based on the longitudinal setting of the axial coordinate system is depicted in Figure 9. As shown, radial lines perpendicular to the jet centerline were cut along the jet flow. ζ and η represent the axial and radial directions in the jet-trajectory-based coordinate system. η was assumed positive inside the jet flow, and x p , z p , and ζ p were the abscissa, vertical, and axial coordinates of the penetration point, respectively. In Figure 9b, the concentration dropped after the jet was injected from the jet nozzle, decayed at a rapid rate, and then slowed down. The predicted centerline concentrations along the jet trajectory agreed well when the grid resolution became finer. The centerline concentration based on the longitudinal setting of the axial coordinate system is depicted in Figure 9. As shown, radial lines perpendicular to the jet centerline were cut along the jet flow. ζ and η represent the axial and radial directions in the jet-trajectory-based coordinate system. η was assumed positive inside the jet flow, and xp, zp, and ζp were the abscissa, vertical, and axial coordinates of the penetration point, respectively. In Figure 9b, the concentration dropped after the jet was injected from the jet nozzle, decayed at a rapid rate, and then slowed down. The predicted centerline concentrations along the jet trajectory agreed well when the grid resolution became finer. On the basis of the uncertainty analyses, concentration contours, coherent structures, and centerline concentration, as well as the relationship between the numerical accuracy and computational resources, the buoyant jet in counterflow was appropriate on grid 1 for F = 5 and R = 7.5. The radial concentration profile on grid 1 is explored in Figure 10 to explore the decay of the jet concentration in detail. In the near-field of the jet, such as ζ/d = 1.1 and 3.1, the radial concentration decayed rapidly, both inside and outside the jet flow. Along the jet trajectory, the radial concentration inside the jet decayed slower than that outside. The probable reason is that the jet flow outside was facing the counterflow and strongly sheared, whereas the jet flow inside was forced to roll up and shed. Therefore, the area of the jet flow outside boundary was more difficult to validate, which was consistent with the results of the uncertainty analyses. On the basis of the uncertainty analyses, concentration contours, coherent structures, and centerline concentration, as well as the relationship between the numerical accuracy and computational resources, the buoyant jet in counterflow was appropriate on grid 1 for F = 5 and R = 7.5. The radial concentration profile on grid 1 is explored in Figure 10 to explore the decay of the jet concentration in detail. In the near-field of the jet, such as ζ/d = 1.1 and 3.1, the radial concentration decayed rapidly, both inside and outside the jet flow. Along the jet trajectory, the radial concentration inside the jet decayed slower than that outside. The probable reason is that the jet flow outside was facing the counterflow and strongly sheared, whereas the jet flow inside was forced to roll up and shed. Therefore, the area of the jet flow outside boundary was more difficult to validate, which was consistent with the results of the uncertainty analyses.

Conclusions
V&V procedures were performed using uncertainty estimate methods in the study of the buoyant jet in counterflow. The main conclusions are as follows.
(1) The series of verifications of the numerical solutions for four grids and the validation with the experimental data revealed that FS and FS1 showed better reliability and applicability in the round buoyant jet in counterflow among the seven uncertainty estimators.
(2) The area without validation at the UV level showed a strong relationship with the shear layer between the jet flow and the ambient fluid. The strongly sheared boundary of the jet led to a higher difficulty in reaching the validation level.
(3) The predicted concentration contours, coherent structures, and centerline concentration agreed well with the experimental data when the grid resolution was improved. In particular, the predicted coherent structures reflected the strong shearing between the jet flow and the ambient fluid only at a high level of grid resolution. Experimental, analytical, or numerical benchmark data for the buoyant jet in counterflow were difficult to obtain due to the limitations of the experimental measurement, the improvement in theory and the direct numerical simulation (DNS). In view of these difficulties, the validation was difficult to process, and the corrected numerical solutions were difficult to derive. Nonetheless, the conclusions drawn based on the V&V procedure are still applicable to the study of the buoyant jet in counterflow. As noted in the literature [33,45], the influence of the grid and time-step should be studied. However, in the current work, the influence of the grid resolution at a fixed time-step was investigated to overcome the difficulties of interpreting the uncertainty results by the grid and time-step simultaneously. Moreover, whether the parallel computing affected the calculation accuracy compared with single-core computing must be explored in future works.
The tendency to employ the LES method for V&V has recently emerged. Compared with RANS, LES can deliver more accurate numerical solutions. However, LES has been rarely reported in the literature, probably due to various factors such as turbulent model, boundary condition, subgrid-scale model, grid resolution, or computational resources. Xing [46] recently proposed a general framework for the V&V of LES. Dutta and Xing [47] adopted this framework in the application of channel flows. In the study of buoyant jets in counterflow, the V&V in LES will be applied with the guidelines proposed by Xing [46] and Dutta and Xing

Conclusions
V&V procedures were performed using uncertainty estimate methods in the study of the buoyant jet in counterflow. The main conclusions are as follows.
(1) The series of verifications of the numerical solutions for four grids and the validation with the experimental data revealed that FS and FS1 showed better reliability and applicability in the round buoyant jet in counterflow among the seven uncertainty estimators.
(2) The area without validation at the U V level showed a strong relationship with the shear layer between the jet flow and the ambient fluid. The strongly sheared boundary of the jet led to a higher difficulty in reaching the validation level.
(3) The predicted concentration contours, coherent structures, and centerline concentration agreed well with the experimental data when the grid resolution was improved. In particular, the predicted coherent structures reflected the strong shearing between the jet flow and the ambient fluid only at a high level of grid resolution. Experimental, analytical, or numerical benchmark data for the buoyant jet in counterflow were difficult to obtain due to the limitations of the experimental measurement, the improvement in theory and the direct numerical simulation (DNS). In view of these difficulties, the validation was difficult to process, and the corrected numerical solutions were difficult to derive. Nonetheless, the conclusions drawn based on the V&V procedure are still applicable to the study of the buoyant jet in counterflow. As noted in the literature [33,45], the influence of the grid and time-step should be studied. However, in the current work, the influence of the grid resolution at a fixed time-step was investigated to overcome the difficulties of interpreting the uncertainty results by the grid and time-step simultaneously. Moreover, whether the parallel computing affected the calculation accuracy compared with single-core computing must be explored in future works.
The tendency to employ the LES method for V&V has recently emerged. Compared with RANS, LES can deliver more accurate numerical solutions. However, LES has been rarely reported in the literature, probably due to various factors such as turbulent model, boundary condition, subgrid-scale model, grid resolution, or computational resources. Xing [46] recently proposed a general framework for the V&V of LES. Dutta and Xing [47] adopted this framework in the application of channel flows. In the study of buoyant jets in counterflow, the V&V in LES will be applied with the guidelines proposed by Xing [46] and Dutta and Xing [47] in our future work.