Modeling Nonconfined Density Currents 3 Using 3 D Hydrodynamic Models

5 Abstract: Density currents generated by marine brine discharges, e.g., from desalination plants, can have a negative impact on marine 6 ecosystems. It is therefore important to accurately predict their behavior. Predictions are often made using computational hydrodynamic 7 models, which should be validated using field or laboratory measurements. This paper focuses on the setup and validation of three8 dimensional (3D) models for estimating the transport and mixing processes that occur in these types of flows. Through a comprehensive 9 sensitivity analysis based on the reproduction of several laboratory-generated density currents, a set of recommendations are made regarding 10 the modeling aspects, including the domain discretization, the treatment of momentum at the density current source, the hydrostatic 11 hypothesis and the selection of turbulence closure models. Finally, the proposed numerical model setup is validated using different exper12 imental data showing good agreement in terms of the main variables considered: errors of less than 1.3% for dilution and of 6% for velocity. 13 This study serves as a first step toward the full validation of these 3D hydrodynamic models for the simulation of field-scale density currents. 14 DOI: 10.1061/(ASCE)HY.1943-7900.0001563. © 2018 American Society of Civil Engineers.

the TCM (e.g., Eidsvik and Brørs 1989;Brørs and Eidsvik 1992;Choi and Garcia 2002).The κ-ε model has also been applied to density currents plunging into reservoirs by Farrell and Stefan (1988) and Bournet et al. (1999).In recent years, a number of direct numerical simulations (DNSs) of density currents have been reported in the literature (e.g., Härtel et al. 2000;Lowe et al. 2005; Birman et al. 2005;Cantero et al. 2006Cantero et al. , 2007)).These more sophisticated simulations are capable of capturing interfacial vortex dynamics such as Kelvin-Helmholtz instabilities and the formation of lobe-cleft structures at the current head.Other authors such as Patterson et al. (2005Patterson et al. ( , 2006) ) have conducted simulations of axisymmetric density currents using implicit large eddy simulations (LESs) (Almgren et al. 1996) relying on the use of subgrid scale modeling (SGS).Nowadays, there are several remarkable studies focused on LESs of different kinds of density currents (e.g., Ottolenghi et al. 2016aOttolenghi et al. , b, 2017aOttolenghi et al. , 2018)).However, DNS and LES are still prohibitively expensive in terms of computational time and especially when considering field-scale simulations.
An alternative to these models based on the resolution of the hydrodynamic equations is given by the lattice Boltzmann method (LBM), defined in the framework of the kinetic theory, which describes the flow in terms of probability density functions.The Focusing on the most developed hydrodynamic equationsbased models mentioned and their application, integral models can be considered adequate for field-scale practical studies of water resources management where a coarse approximation of the characteristic flow quantities at an equilibrium state may be sufficient.Conversely, the most complex numerical approximations (e.g., DNS and LES) are highly time-demanding computationally and are typically applied at the laboratory scale under controlled conditions.However, intermediate complexity 3D hydrodynamic models can be used to solve field-scale applications; as an example, Bombardelli and Garca (2002) where production terms denoted by the shear P and buoyancy G 219 values are estimated by (in our case the salinity), the present study considers a second-order central-upwind scheme originally based on the Kurganov and Petrova (2007) scheme and adapted by Bourban (2013).For κ-ε variables, the default advection scheme established by TELEMAC-3D is the method of characteristics.However, a sensitivity test based on the advection scheme for these variables was carried out in this study.Finally, note that relative accuracies, number of iterations, and preconditionings for each variable required for the recommended iterative solver (conjugate gradient) are established in accordance with the recommendations made in the TELEMAC-3D user manual (LNHE 2007).

Experimental Databases
Two experimental databases were used to establish the validated shows a schematic diagram of the experimental setup outlined.The main initial characteristics of this set of laboratory-generated density currents are listed in Table 1 and are outlined in Fig. 2. Note that Case C1 was selected as the reference case on which all other cases were modified by one of the initial parameters (thickness h 0 , flow rate Q 0 , slope α, or density difference ρ a − ρ 0 ).This way, a comprehensive set of laboratory-generated density currents was carried out to experimentally characterize these kinds of flows under different flow-expected conditions.
To obtain high quality velocity and concentration measure-   The specific methodology of the sensitivity analysis involved several simulations of the base application Case C1 (Table 1) by varying by one numerical aspect while keeping the remainder unchanged.Table 2 summarizes the numerical aspects and the variations considered in this study.The significance of the numerical aspect considered for the prediction of density current behavior was analyzed by comparing the numerical results of characteristic magnitudes such as horizontal density current spreading (b), and velocity (U) and dilution (S) evolution within the density current.
Dilution was calculated using the following expression: Table 1.Main characteristics of laboratory configurations      1).Furthermore, Figs.2001).This method is based on equations for the conservation of volume (Eq.8) and solute mass (Eq.9) where U, β, and h are the mean values for velocity, fractional density [β ¼ ðρ − ρ a Þ=ρ a ] generated from saline concentrations, and the density current thickness for each location, respectively.
The value of the current thickness h is defined as the distance from the base where the salinity concentration is less than 10% of the maximum concentration value at the corresponding location.
Combining Eq. ( 8) with Eq. ( 9 Fig. 4 shows the value of the average normal entrainment coefficient E N relative to the number of vertical layers within the density current N Zh ¼ h 0 =Δz min for different horizontal discretizations and for both kind of tests, namely with and without the application of the TCMv.The variable E N represents a single average value calculated from X=L P ¼ 15 to eliminate the effects of the density current's initial adjustment on the normal flow state. Fig. 4 reveals that with the exception of those for the most coarsely resolved tests, the global modeled entrainment rates (i.e., entrainment rates of cases with the TCMv set to ML) converge on the order of 10 −2 while numerical entrainment rates (i.e., entrainment rates of cases with the TCMv turned off, NT) decrease nearly exponentially as the vertical resolution increases.For tests involving values of h 0 =Δz min higher than 8, the difference between the two entrainment rates reaches close to or higher than an order of magnitude of 10 −2 (specifically, for h 0 =Δz min ¼ 8, global modeled entrainment rates are less than 1.2 × 10 −2 and numerical rates are higher 4 × 10 −3 ).Thus, for cases involving more than eight layers within the density current (N zh ¼ h 0 =Δz min > 8), numerical mixing ceases to dominate global modeled mixing.We note that when Q an injection velocity V.This section presents our analysis of the 551 effects of these numerical aspects on the density current behaviour.

552
The numerical parameters for these simulations are presented in 553   For the turbulence closure model of the horizontal direction (TCMh), it is common for hydrodynamic models to use a constant turbulence model.Thus, the user must calibrate the horizontal eddy viscosity value ν th depending on the particular flow being studied and depending on domain discretization (Madsen et al. 1988).However, in taking advantage of options programmed into TELEMAC-3D, simulations varying the TCMh between the constant and Smagorinsky models were compared.For the studied case, while using the Prandtl mixing length TCMv, no appreciable differences were observed in the results.For these tests, ν th was defined according to the variable grid resolution (applying scaleddown absolute values of 4.7 × 10 −5 to 1.6 10 −4 m=s 2 ), and the calibration parameter of the Smagorinsky model (C s ) was defined as 0.1, a common value for anisotropic flows (e.g., flow in a canal).
The κ-ε model was also set as the TCMh model, but in this case the κ-ε model was mandatory for the TCMv, and so the influence of the TCMh could not been extracted.In addition, the simulations became so time-consuming and unstable (due to numerical prob-  1993).In addition to our numerical results, measured data obtained via PIV-PLIF techniques are also displayed.Fig. 6 shows cross profiles derived from different simulations with constant TCMv based on varying ν tv values of ν th =1000 to ν th =10.In these simulations the eddy diffusivity value Γ is defined by the Schmidt number (σ c ¼ ν tv =Γ v ), which is considered to take values of close to one for these types of flows.Both the velocity and salinity concentration cross profiles show that the simulation ν tv ¼ ν th =10 better fits the experimental data.However, the simulation does not represent the approximate shape of the first cross sections, i.e., those with higher momentum and concentration values.For the remainder of the cross sections, the velocity agrees with the experimental data while the dilution is underestimated (i.e., higher concentrations than were expected).As is shown by the similarity cross profiles, the normalized mean horizontal velocity and mean concentration cross profiles collapse well into single profiles, which agree with the experimental similarity profiles.While the velocity cross profiles match fairly well with experimental data for the whole density current body, the concentration profiles significantly underestimate dilution levels found in the region close to the bottom for areas located closer areas to the discharge.Taking the similarity cross profiles into account, i.e., the profile shape profile, the constant, and the κ-ε TCMv are the models that flows (e.g., far-field jets and plumes for which the velocity difference across the flow represents only a small fraction of the convection velocity), P was found to be significantly different from ε, and c μ was found to take different values (Rodi 1975).Rodi (1972) correlated experimental data and proposed a function of c μ ¼ fðP=εÞ that is only valid for thin shear layers (similar to the density currents studied).
To study effects of these empirical constants on density currents modeling, the results of several simulations varying the values of these constants are analyzed.The range of values for c 3ε and c μ were chosen based on state-of-the-art findings, namely c 3ε values of 1-0.6 and c μ values of 0.09-2.5.Another important parameter is the Schmidt number σ c , which contributes to the eddy diffusivity definition and which is valued at between 0.7 and 1.
However, due to its lesser impact on the results compared to the c 3ε and c μ constants' impact, it was assumed to be equal to 0.7, a common value for heat and salinity transport.Based on the experimental results, comparisons are drawn with results obtained from different simulations using the root mean-square absolute error where  4.
Table 4 shows that the results of the simulation with c 3ε equal to 0.7 and c μ equal to 0.2 offer the smallest errors, namely 0.135 psu for salinity and 0.005 m=s for velocity (2.8% and 4.8%, respectively).Fig. 9 shows the downstream variation and similarity cross profiles for the simulation.As is shown, the numerical and experimental data agree fairly well, both in terms of shapes and absolute values.
The optimal value of c 3ε obtained from this study (∼½0.7 − 0.6)    5. Apart from the previously calculated relative errors regarding the evolution of main variables in the cross profiles, in applying these recommendations to Case 1, a relative error of the front position (NRMSE X f exp ) of less than 6% is obtained.
For the evolution of the density current in the streamwise direction, wherein the dilution of the current is the lowest (line 0°),   spreading and hinder lateral spreading (i.e., higher front velocities and lower half-widths), and a higher density difference favors longitudinal and lateral spreading (i.e., higher front velocities and higher half-widths).• Full momentum source specification is advisable, i.e., both flow rate and velocity information (Q and V) should be detailed.
• The hydrostatic hypothesis is considered to be appropriately assumed, as it significantly (1.5-2 times) reduces computation time and because not applying this hypothesis does not generate significantly improved results.

1 Environmental Hydraulics 2 Institute
are generally referred to 17 as density or gravity currents, are continuous underflows that 18 travel downslope due to their negatively buoyant characteristics, 19 i.e., because they are heavier than the surrounding fluid.This 20 phenomenon occurs widely in natural environments and is caused 21 by either human activities or natural processes (Simpson 1997; 22 Huppert 2006).Currently, in coastal and marine environments, 23 some of the most common density currents are those generated by 24 brine discharge from desalination plants.Hodges et al. (2011) make 25 an analogy between the behavior of a natural salt wedge and such 26 brine discharges into shallow waters, both of which are governed 27 by the density difference, by the hydrodynamics of the surrounding 28 area (Shao et al. 2008), and by the bottom slope.Due to the po-29 tentially negative impacts of these human-induced currents on the 30 environment (Lattemann and Höpner 2008; Sánchez-Lizaso et al. 31 2008; Laspidou et al. 2010; Dawoud and Al Mulla 2012) there is a 32 growing interest in obtaining accurate predictions of their behavior.33 Dense underflows have been widely investigated in labora-34 tory experiments (Alavian 1986; Garcia 1993; Gerber et al. 35 2011; Ottolenghi et al. 2017b) and field studies (Hebbert et al. 36 1979; Dallimore et al. 2001; Fernandez and Imberger 2006; Hodges et al. 2011).Major efforts have also been made to predict the behavior of these currents from different modeling techniques and through comparisons with previous laboratory experiments (Choi 1999; La Rocca et al. 2008; Lombardi et al. 2015; Sciortino et al. 2018).As a broad classification, two modeling techniques are available for studying these flows numerically using the hydrodynamic equations (i.e., continuity, momentum, and transport equations): integral models and those that determine the vertical structure of the flow.The integral model for density currents was first introduced by Ellison and Turner (1959) and was further developed by Alavian (1986) and Parker et al. (1986) among others, primarily focused on turbidity currents (Akiyama and Stefan 1985; Parker et al. 1986; Garcia 1993; Choi and Garcia 1995; Bradford et al. 1997; Imran et al. 1998; Choi 1999).In general, these integral models assume a hydrostatic pressure distribution within the density current and use vertical depth-integrated equations.They have been designated single-or double-layer shallow water models depending on whether they consider only the heavier layer (e.g., Ungarish 2007a, b; La Rocca et al. 2008; Lombardi et al. 2018) or divide the entire depth into two layers (e.g., Ungarish 2008; La Rocca and Pinzon 2010; La Rocca et al. 2012).Note that although it has been demonstrated that these integral models are capable of providing good results under laboratory controlled conditions, they are not capable of taking into account the complexity of real environmental conditions that may occur in nature and may affect the evolution of the density current flow.Conversely, there are numerical studies that have used models based on the resolution of three-dimensional (3D) Navier-Stokes (N-S) equations with differing degrees of simplification to solve the vertical structure of density current flows.Specifically, the vertical distribution of the main variables of density currents has been numerically analyzed from hydrodynamic models that solve N-S equations taking into account the Reynolds approximation (Reynolds-averaged Navier-Stokes equations or RANS).In these applications, a turbulence closure model (TCM) estimates the Reynolds stress in conjunction with wall functions.Stacey and Bowen (1988a, b) solved the vertical distribution of onedimensional turbidity currents using a mixing length model as the TCM.Other authors have employed the κ-ε model (Rodi 1984) as "IH Cantabria," Universidad de Cantabria, 39011 Santander, Spain (corresponding author).ORCID: https:// orcid.org/0000-0002-1987-2605.Email: perezdb@unican.es simplicity and versatility of the LBM has encouraged its development in the computations fluids dynamics within the last decade (e.g., Aidun and Clausen 2010).Regarding the reproduction of complex flows such as density currents, Rocca et al. (2012) developed a LBM for two-layered shallow-water flows by considering two separate sets of LBM equations, one for each layer, obtaining good agreement between the LBM numerical results and the experimental results when the evolution of the flow does not depend on the viscosity.Recently, Ottolenghi et al. (2018) revels that this alternative can be also applied for three-dimensional numerical simulations of density currents for different Reynolds number by implementing an equivalent large eddy simuation model in the LBM framework.Nevertheless, although LBM has been successfully applied to simulate density currents generated by laboratory experiments, both considering regular and complex geometries (e.g., Prestininzi et al. 2016), significant research still needs to be done to strengthen the LBM for simulating real field-scale density currents.
assessed the potential development of density currents in the Chicago River while capturing their spatial variability.Kulis and Hodges (2006) carried out a layernumber sensitivity test to numerically simulate density currents of the Corpus Christi Bay in Texas using a sigma-coordinate 3D hydrodynamic model based on RANS equations and while taking into account the Boussinesq approximation and the hydrostatic hypothesis.Applying similar 3D hydrodynamic models, Firoozabadi et al. (2009) and Mahgoub et al. (2015) simulated density currents and validated their results against some laboratory measurements.Nevertheless, to the authors' knowledge, none of the reviewed studies have been fully validated, i.e., considering both horizontal spreading and the vertical structure of the main flow variables (velocity and concentration).In addition, none of these studies provide 140 recommendations to accurately reproduce these kind of flows.This 141 will be very useful in practical purposes such as the design of brine 142 discharges into seawater, which have to meet strict water quality 143 criteria regarding the salinity concentration far from the discharge 144 point (e.g., Sánchez-Lizaso et al. 2008).145 The current paper focuses on establishing a suitable setup of a 146 3D hydrodynamic model based on RANS equations, while taking 147 into account the Boussinesq approximation and the hydrostatic or 148 nonhydrostatic pressure hypothesis, for simulating nonconfined 149 density currents.The study numerically reproduces a set of labo-150 ratory experiments carried out in the Environmental Hydraulics 151 Institute (IHCantabria) by applying advanced optical techniques 152 (Pérez-Díaz et al. 2016, 2018), as well as reproducing other ex-153 periments under different flow conditions presented by Choi and 154 Garcia (2001).The numerical simulations are fully calibrated and 155 validated against laboratory measurements by comparing the main 156 flow and mixing characteristics.Therefore, the present paper out-157 lines an optimum modeling setup that predicts the behavior of these 158 types of flows using 3D hydrodynamic models.In addition, taking 159 into account that these models are also capable of simulating real 160 environmental conditions that may affect field-scale density cur-161 rents (Shao et al. 2008), the findings of this study are also presented 162 as a starting point for future field-scale studies.163 The paper is presented as follows.First, the methods used are 164 introduced, then the calibration results obtained through a compre-165 hensive sensitivity analysis are presented and discussed.Third, the 166 validation results are presented, and finally conclusions are drawn.diffusivities according to the grid resolution 201 and characteristic velocity of the type of flow motion studied 202 (Madsen et al. 1988).Mixing length and Smagorinsky TCMs are 203 based on the mixing-length concept proposed by Prandtl.While 204 mixing-length TCMs such as the standard Prandtl model (Rodi 205 1984) and the Nezu and Nakagawa model (Nezu and Nakagawa 206 1993) are only applied as vertical TCMs, the Smagorinsky model 207 (Smagorinsky 1963) is a subgrid turbulence model wherein mixing 208 length is dependent on the grid and on a dimensionless coefficient 209 according to the type of flow involved (anisotropic or isotropic 210 flow).Finally, the most complex TCM used in this study is the 211 two-equation κ-ε model.In this paper, eddy viscosities are evalu-212 ated by applying the following Kolmogorov-Prandtl relationship 213 per direction: setup of the TELEMAC-3D model for predicting the behavior of nonconfined density currents.The main and largest database was generated from a set of laboratory-generated density currents tested at IHCantabria's facilities.These laboratory experiments, presented in Pérez-Díaz et al. (2018), consisted of saline density currents that evolved over a gentle-slope (α) plastic-material base within a 3 × 3 × 1 m 3 test tank filled with freshwater to simulate the receiving body (both fluids at the same temperature but with different saline concentration).The constant-flux saline effluent was discharged through a rectangular height-adjustable slot (b o × h 0 ) at the base [Fig.1(b)], simulating the start of far field region of mixing that commonly forms when brine is discharged by submerged jets (Papakonstantis and Christodoulou 2010; Palomar 2014).Fig. 1(a) ments in the longitudinal profile of the density currents tested, the mentioned set of experiments used nonintrusive laser optical techniques, namely particle image velocimetry (PIV), and planar laser induced fluorescense (PLIF).The PIV technique consists of capturing the movements of small seeding particles within the flow between consecutive laser pulses that illuminate the flow, while PLIF consists of an indirect way to measure concentration based on capturing the light re-emision of a fluorescent dye once it is illuminated by the laser at determined wavelength.In the mentioned laboratory study, to capture the maximum covering area (1,400 mm), two LaVision Imager ProX 4 (CCD) cameras with a resolution of 2048 × 2048 pixels were located adjacently [Fig.1(a)] with an overlapping zone.This configuration together with the appropriate selection of PIV and PLIF calibration parameters enabled the accurate measuring of concentrations and velocities.In addition, to study the lateral and front spreading [Fig.2(b)], a set of photos were taken with a camera located at the upper part of the test tank that was able to capture the whole plan view.Further information on the experimental setup and on the advanced PIV-PLIF measurement techniques used to obtain this experimental database can be found in Pérez-Daíz et al. (2018).The other experimental database used in this study was developed by Choi and Garcia (2001), who published the results of a set of density current experiments carried out in a 2.44 × 3.66 × 1.22 m 3 tank.These experiments were chosen for this study because they generated nonconfined saline density currents with different initial conditions and because the authors provided all reproduce them numerically.Choi and 304 Garcia (2001) studied the spreading rates of nonconfined density 305 currents on sloping beds while varying the density difference and 306 slope angle.The sloping bed used was composed of fiberglass, 307 i.e., with a roughness equivalent to that of glass or Plexiglass.The 308 initial flow parameters of these experiments are clearly detailed in 309 317predictions.The simulations of the calibration stage were used to 318 create a numerical setup for predicting the behavior of these types 319 of flows.This calibration stage consisted of a sensitivity analysis of 320 key numerical aspects that may influence the numerical prediction of density currents such as the domain discretization, the treatment of momentum at the density current source, the hydrostatic hypothesis, and the selection of the TCM.Once the sensitivity analysis was conducted, the proposed setup was validated from a final set of simulations that reproduced the complete set of laboratorygenerated density currents.

Fig. 1 . 7 341
(a) Schematic diagram of the PIV-LIF experimental setup; and (b) photographs of the discharge device.= initial salinity concentration of the source; C a = 338 surrounding fluid salinity concentration; and C = salinity concen-339 tration at the study point within the density current.340 Initial Model Setup The initial and boundary conditions of the TELEMAC-3D appli-342 cations were defined with the aim of numerically simulating the real conditions of the experimental setup.Accordingly, nodes across the domain were initialized with constant corresponding elevation values and with zero velocities (stagnant receiving water).For the boundary conditions, the free surface, open boundaries (i.e., liquid boundaries), and base and rigid walls were taken into account.At the free-surface boundary, a rigid-lid approximation that considers zero gradients and zero fluxes perpendicular to the boundary was applied (i.e., ∂U i =∂z ¼ ∂κ=∂z ¼ ∂ε=∂z ¼ ∂T=∂z ¼ 0).At the open boundaries, streamwise gradients of all of the variables (i.e., velocities, tracer, and fluxes) were set to zero and a prescribed elevation was applied.Strictly speaking, boundary conditions for velocities subjected to rigid wall reflect a no-slip condition (i.e., Dirichlet conditions U i ¼ 0).However, due to the presence of turbulence and of a boundary layer, the velocity close to the base quickly becomes non-zero, and the no-slip condition is replaced with tangential stress [i.e., τ ¼ μð∂ Ũ=∂nÞ ] due to friction subjected to the base.This tangential stress is replicated by a turbulence model for the bottom using the friction or shear velocity τ ¼ −ρðU Ã Þ 2 and the distance to the bottom z.Assuming that the flow is hydraulically rough [i.e., the characteristic roughness size of the base is greater than the thickness of the viscous sublayer (Hervouet 2007)], the velocity profile close to the base was defined by a logarithmic law function of the Nikuradse coefficient k s representing the roughness size.As the base material used for the experiments was plastic, a Nikuradse coefficient of 10 −5 m was set.For the rigid vertical walls, slip conditions were assumed (i.e., without friction).Tracer concentration gradients were also set to zero for the rigid walls (base and vertical rigid walls).For the turbulent kinetic energy κ and its dissipation rate ε, the boundary conditions defined by Rodi (1984) for rigid walls were applied.Note that while in other numerical experiments (e.g., Firoozabadi et al. 2009) dense fluid enters the domain through an open-liquid boundary with the slot dimensions, in this study the saline flow rate was determined based on a series of discrete source terms.The number of source terms was determined via the slot dimensions and the domain discretization.This way, the findings of this study can be applied to future field applications (e.g., brine discharges from desalination plants) wherein the saline outflows are typically located within the study domain rather than along boundaries.As these types of hydrodynamic models and corresponding grid tools are designed and generally configured (e.g., accuracy levels and number of iterations) to model coastal and ocean processes, laboratory tests should be scaled up to prevent numerical problems from emerging.Froude similarity [i.e., the relevant forces are the inertial and gravity forces(Heller 2011)] and mechanical similarity (i.e., geometric, kinematic, and dynamic similarity) are expected to be achieved and thus fewer scale effects are expected to be obtained.However, as the Reynolds number, R, of the case studied was not sufficiently high (i.e., R ≫ 2000) to directly neglect the viscous force, a previous sensitivity analysis varying 10 the scale Sc F ¼ L NUM =L LAB (considering the Froude similarity) was performed.Scale factors studied include the following: Sc F ¼ 1; 10; 20; 50; 80; 100; 200; 1000.This analysis showed that after converting all of the results to the same scale, the relative difference of main quantities (concentration and velocity) for the scalesensitivity cases (at geometrically equivalent locations) was always lower than 2%.This negligible difference attributable to numerical and scale effects shows that Froude similarity can safely be assumed.Thereafter, a Sc F ¼ 100 scale factor was used so that the modeled density currents have similar characteristics to the far field region of brine discharges.
computational grid design is critical to simulate real physical processes while avoiding artificial numerical effects.As a general rule, a fine grid domain discretization (i.e., high resolution) is needed when high spatial-temporal gradients of the variables modeled are anticipated.For negatively buoyant density current flows, high variability areas are the zone closest to the bottom (high vertical gradients) and the surroundings of the discharge location (high vertical and horizontal gradients).This section presents the results of our domain discretization sensitivity tests, in which the horizontal and vertical discretization were varied, leaving the rest of the numerical aspects constant according to Table3.We note that sigma-layer coordinates are necessary to resolve the vertical structure of the density current.The number of sigma planes and their spacing was therefore investigated as part of the sensitivity tests.First, sensitivity to the horizontal discretization was studied.For this purpose, vertical discretization close to the bottom was set as Δz ¼ ho=20, following the recommendations ofKulis and Hodges (2006).As shown in Table3, two horizontal discretization parameters were defined, namely Δx 1 , the highest resolution close to the discharge location, and Δx 2 , the lowest resolution in the area furthest away from the discharge location.To compare corresponding patterns of horizontal spreading, Fig.3(a) shows the front positions versus time in the plane of symmetry and in the plane at 45°relative to the symmetry plane (see lines 0°and 45°of Fig. 2).These spatial and temporal quantities are normalized by the characteristic length and time scales of plume-like behavior flows (Chu and Jirka 1987; Choi and Garcia 2001),

Fig. 3 (
Fig. 3(a) shows that for values of Δx 1 smaller than half of the horizontal slot dimension, i.e., bo=2, the spreading converges to similar values for all of the cases.From the longitudinal profile of the normalized maximum velocity shown in Fig. 3(b), higher differences can be observed in the region close to the discharge point.Using the test with the smallest value of Δx 1 as a reference (a value close to the one), it is evident that a discretization of at least Δx 1 < b 0 =4 is needed to capture flow motion in the discharge surroundings.Additional tests varying Δx 2 from b 0 =4 to 3b o show that in areas positioned far from the discharge location, the horizontal discretization can be coarser and the results do not show significant differences.In this way, the computational time can be minimized, but due to the high levels of horizontal discretization variability, special attention must be paid to the TCMh when it is set as the constant model.In such cases, it is recommended that different eddy viscosity and diffusivity coefficient values are applied along

Fig. 3 .
Horizontal discretization sensitivity tests: (a) dimensionless F3:2 front position versus dimensionless time for lines 0°(continuous line) F3:3 and in line 45°(dashed line); (b) longitudinal profile of normalized F3:4 maximum velocity; and (c) longitudinal profile of minimum dilution.from that depth to the surface, spacing was gradually increased by a factor of 1.5.The results of these tests (with and without TCMv) were used to evaluate the relationship between global modeled and numerical mixing (i.e., mixing only due to numerical effects not related to the real physics of the process) at different vertical discretizations.For cases in which the TCMv was turned off (no turbulence NT tests), all mixing can be attributed to molecular Brownian motion.As in these NT tests the kinematic viscosity (i.e., molecular) coefficient was set to 10 −9 m 2 =s, vertical mixing should effectively be zero.Therefore, any mixing observed in the NT tests can be attributed to numerical effects or so-called numerical mixing (or numerical diffusion).Conversely, for cases in which the TCMv was turned on, mixing can be attributed to the combination of turbulent and numerical mixing, henceforth referred to as global modeled mixing.The aim of this specific analysis is to detect the number of layers from which numerical mixing can be considered negligible as compared with the global modeled mixing.Calculating the vertical entrainment coefficient (E) in the symmetrical longitudinal profile (line 0°) for both types of tests generates a quantitative measure for mixing.The entrainment coefficient is calculated from the salinity longitudinal profile using the method developed by Dallimore et al. ( ) generates equation Eq. (10) where dC=dx is the variation in concentrations when the current travels with the mean velocity.By assuming similarity of the concentration profiles (Parker et al. 1987; Pérez-Díaz et al. 2018), C can be set as 522 considering the horizontal resolution sensitivity tests for the each 523 of the corresponding number of layers, the previous pattern is 524 maintained and always generates higher entrainment rates for tests 525 involving coarser horizontal resolutions.Fig. 4 shows that both 526 horizontal and vertical resolutions affect the global modeled mix-527 ing; the higher the resolution, the lower the degree of numerical 528 diffusion.To ensure that the effects of numerical diffusion are neg-529 ligible, our testing shows that h 0 =Δz min values of 16-20 are 530 needed.As shown in Fig. 4, this ensures rates of numerical diffu-531 sion that are approximately an order of magnitude below the global 532 modeled mixing.533 Source Input and Hydrostatic Hypothesis 534 Having defined a suitable computational grid, numerical aspects 535 that can also affect the initial region of the density current are stud-536 ied.As cited in the numerical model description, the TELEMAC-537 3D model allows both hydrostatic and nonhydrostatic pressure 538 fields to be assumed.The hydrostatic assumption is mainly valid 539 for anisotropic flows for which scales of motion are substantially 540 larger in the horizontal than in the vertical.As density currents are 541 primarily horizontal flows, the hydrostatic hypothesis should be 542 appropriate.However, as Mahgoub et al. (2015) show, vertical 543 accelerations may not be negligible relative to gravitational ac-544 celeration for the starting region of the density current due to 545 local effects found in the near field region.Numerically, with a 546 nonhydrostatic pressure field, the vertical momentum equation 547 is fully solved without simplification, and the pressure equation 548 is split up into hydrostatic pressure and a dynamic pressure terms.549 TELEMAC-3D also allows the specification of a liquid flow rate 550

559Fig. 5 1 Fig. 4 . 15 Q
Fig.5shows the effects of the source specification type and lems found at open liquid boundaries based on Neumann boundary conditions for the κ and ε equations) that this option for the TCMh was not considered.Due to the strong stratification associated with density currents, the turbulence closure model of the vertical direction (TCMv) is a key numerical aspect for accurately predicting their behavior.To analyze the TCMv's influence on the vertical structures of the flow, Figs.6-8 show downstream variations in the velocity (the main velocity component, U x ) and salinity (C) cross profiles found along the symmetry plane (see line 0°of Fig. 2) for each TCMv case in the constant, mixing length, and κ-ε models.For all of the simulations, the TCMh is held as constant according to the results presented previously.Furthermore, these figures present normalized cross profiles that collapse into a single profile, called a similarity profile (Ellison and Turner 1959; Parker et al. 1987; Garcia

Fig. 7
Fig. 7 shows cross profiles for the TCMv mixing length models (both the Prandtl forumation and the Nezu and Nakagawa formulation).For these cases, the damping function addressed by Munk and Anderson (1948) is used to govern vertical mass and momentum exchanges.The results of the two simulations are almost identical, showing logarithmic velocity and concentration profiles for the most parts of the density current body.In these cases, the profiles are the result of agreement between the ML model, the damping function and the turbulence model for the bottom based on the previously defined roughness size (i.e., the Nikuradse coefficient).

F5: 1 Fig. 5 .
Source-input and hydrostatic-hypothesis sensitivity tests: F5:2 (a) dimensionless front position versus dimensionless time for lines F5:3 0°(continuous line) and 45°(dashed line); and (b) longitudinal profile F5:4 of normalized maximum velocity.addition, the normalized mean horizontal velocity and mean 654 concentration cross profiles collapse into single profiles but they 655 disagree with the experimental similarity profiles.This shows that 656 the presented TCMv option does not capture the vertical shape of 657 the analyzed density current.

Fig. 8 1 Fig. 6 .
Fig.8shows the cross profiles derived from two simulations with the κ-ε TCMv applied.The two simulations differ in terms of advection scheme (AdSch κε ) applied.Whereas in one case the method of characteristics was used, for the other case the most conservative scheme programmed into the TELEMAC-3D was

F7: 1 Fig. 7 .
Downstream variation cross profiles for cases with TCMv ¼ ML: (a) horizontal velocity; (b) salinity concentration; and the corresponding F7:2 similarity cross profiles for (c) horizontal velocity; and (d) salinity concentrations.Note that both lines overlap in all figures.locations studied, although the shape of the profile obtained 672 is noticeably different from one case to the other.Regarding the 673 similarity profiles, from which a proper shape comparison can 674 be undertaken, they show that the simulation applying the most conservative advection scheme presents better agreement with the experimental similarity curves.

F8: 1 Fig. 8 .
Downstream variation cross profiles for cases with TCMv ¼ κ-ε: (a) horizontal velocity; (b) salinity concentration; and the corresponding F8:2 similarity cross profiles for (c) horizontal velocity; and (d) salinity concentrations.vertical structure of the density current.As the constant model only uses ν tv as a calibration parameter and as the results presented in Fig. 6 for ν tv ¼ ν th =10 are its best results, we focus on the κ-ε TCMv.As described in the previous section, the κ-ε model includes several empirical constants obtained via data fitting for a broad range of flows.Of these empirical constants, the c 3ε and c μ affect the modeling of density currents, and their values have been the subject of debate.Hossain and Rodi (1982), Rodi (1987), and Choi and Garcia (2002) have suggested that c 3ε values of 1-0.6 show a good agreement with experimental results for density currents.Conversely, the standard value of the other controversial empirical constant (c μ ¼ 0.09) was used on the basis of experiments on flows for which production P and dissipation ε of the turbulence energy were in approximate balance.For weak shear x exp = experimental values of the variable studied (in this case the salinity concentration and horizontal velocity); x num represents the corresponding numerical value; and N = number of pairs of comparable values.In this case, due to the large amount of data obtained from the PIV-PLIF experiments, N corresponds to the number of cross profiles multiplied by the number of data points measured within the density current thickness (which varies depending on the variable).To estimate the relative error, the normalized RMSE (NRMSE) was obtained by dividing the RMSE by the initial value of the evaluated magnitude (C 0 or U 0 ) and by multiplying this value by 100 to obtain the percentage.As a calibration assessment, the error obtained for each simulation while varying the c 3ε and c μ constants is shown in Table Fig. 10(d) presents the production term (P), which shows very the previous sensitivity analysis, a proposed modeling 775 setup to simulate the behaviour of density currents is shown in 776 Table

F9: 1 Fig. 9 .
Downstream variation cross profiles for cases with calibrated TCMv ¼ κ-ε: (a) horizontal velocity; (b) salinity concentration; and the F9:2 corresponding similarity cross profiles for (c) horizontal velocity; and (d) salinity concentrations.roughly 4% for the main velocity component 783 (NRMSE U 0 ) and is less than 1.3% for the dilution (NRMSE S min ) 784 value.Fig. 11 illustrates the evolution of the main variables de-785 scribed (i.e., the front position, the dilution, and the main velocity 786 component) and their graphical comparison with the experimental 787 data.Moreover, Fig. 12(a) presents the good agreement between 788 the numerical and experimental spreading results for different times 789 using a plan-view comparison of concentration values.The longi-790 tudinal symmetry profile views of the numerical and experimental 791 concentration results are shown in Figs.12(b and c), for consis-792 tency with the previous results presentation.793 To ensure that the proposed modeling setup is valid for den-794 sity currents generated under different flow conditions, consider-795 ing turbulent flow and supercritical regime (i.e., R>1000 and 796 F>1), it was applied to the complete set of experimental density 797 currents (Table1).Comparisons were made between the experi-798 mental and numerical results from the minimum dilution obtained 799 from the last section of density currents (S minF ), which represents the result of the transport and mixing processes and the parameter on which the environmental regulations are commonly based.

Fig. 13 (
a) shows that the predicted evolution of the dimensionless half-width (b 1=2 ) agrees well with the experimental data.In addition, Fig. 13(b) compares longitudinal spreading values based on front velocity evolution.In this case, velocity is fairly well reproduced when the initial high momentum region is developed.Beyond this comparison, conclusions regarding effects of the different densities and slopes on the described behavior can be drawn from the graphs.In agreement with previous studies (Choi and Garcia 2001; Alavian 1986), steeper slopes favor rapid longitudinal

•
sensitivity and validation analysis based on the reproduction of several laboratory-generated density currents, this study shows that suitably configured 3D hydrodynamic models can simulate the behavior of saline density currents with a high level of accuracy.Moreover, as these laboratory-generated currents were scaled up to prevent numerical effects and to have similar characteristics to the far field region of brine discharges, the sensitivity and validation analysis as well as the recommendations resulted are extended to real field-scale saline current flows.The main numerical guidelines for solving such flows using these models are as follows: Variable spacing horizontal discretization (e.g., unstructured grids) is recommended as a means to obtain high resolutions close to the source.In this way, a lower resolution can be set to the farthest region of the source, rendering application more computationally efficient.Assuming a slot-shaped source (b 0 × h 0 ), a common approximation for the beginning of the far field region of brine discharges, minimal spacing of at least equal to the width of the slot b 0 must be applied to ensure the momentum and mass conservation as much as possible.Specifically, for the cases analyzed in this study, a spacing of greater than b 0 =4 is recommended.•In the vertical direction, a high resolution must be applied to minimize vertical numerical diffusion.For the type of density currents studied in this study and considering the a slot-shaped source, vertical spacing should be at least h 0 =16.As this fine resolution cannot be maintained along the entire water column (i.e., there would be too many numbers of layers), gradual vertical spacing with the highest resolution close to the bottom (within the density current body) is recommended.In such cases, a sigma layer coordinate (i.e., terrain following) for vertical domain discretization should be applied to keep the finest layers at the bottom.

F11: 1 Fig. 11 .
Comparison between the numerical and the experimental F11:2 results: (a) dimensionless front position versus dimensionless time for F11:3 lines 0°(continuous line) and 45°(dashed line); (b) longitudinal profile F11:4 of normalized maximum velocity; and (c) longitudinal profile of mini-F11:5 mum dilution.F12:1 Fig. 12.Comparison between the numerical and the experimental spreading results for different times: (a) plan view (experimental results are shown F12:2 as the dashed line); (b) profile view of the modelled results; and (c) profile view of the LIF experimental results.

•• 1 Fig. 13 .
The constant model is recommended as a horizontal turbulence closure model (TCMh) varying eddy coefficient values according to the grid resolution.Several vertical turbulence closure models (TCMv) such as constant, mixing length, or κ-ε models can be successfully applied.However, given the demonstrated influential role a TCMv has on the numerical simulation of such flows, calibration should be applied to each case study.Specifically, for the cases analyzed in this study, the calibrated κ-ε model (empirical constants c 3ε and c μ are equal to ∼0.7 and ∼0.2, respectively) in conjunction with the most conservative advection scheme for κ-ε equations F13:Comparison between the numerical and experimental results of Choi and Garcia (2001): (a) dimensionless maximum half-width versus time; F13:2 and (b) velocity front.

Table 1 of
Choi and Garcia (2001).Specifically, experiments from 310DEN1 to DEN9 (see Table1of Choi and Garcia 2001) were used 311 in the present study.312 Methodology 313 The present study involved of a series of more than 90 numeri-314 cal simulations that reproduced the aforementioned laboratory-315 generated density currents.Simultaneously, data collected from 316 the experiments were used to calibrate and validate the numerical

Table 2 .
Numerical aspects and their alternatives a V = injection velocity.b Constant model.c Smagorinsky model.d Mixing-length models.e Method of characteristics.f Second-order Kurganov and Petrova scheme.

Table 3 .
Numerical aspects in each of the sensitivity tests considered

Table 3 .
554Note that while physically there is only one source of h 0 × b 0 555 dimensions, due to domain discretization, it is numerically trans-556 formed in several discrete source terms.For instance, 160 discrete 557 sources of a liquid flow rate of Q 0 =160 are needed for the speci-558 fications shown in Table3(Δx 1 ¼ bo=8 and Δz min ¼ ho=20).

Table 4 .
Calibration of empirical coefficients c μ and c 3ε of the κ-ε model

Table 6
minF N maintain the correlation revealed by the experiments S minF E .That is, establishing the value of C1 as the base case S minF Nb (corresponding to the case analyzed in the sensitivity analysis), the rates S minF N =S minF Nb agree with the corresponding experimental ratio S minF E =S minF Eb , showing that steeper slopes (C4 shows the comparison between the experimental and numerical results, with subindices E and N, respectively.As is shown, a good agreement between the experimental and numerical values was obtained for all cases.Furthermore, the modeled values S

Table 1 of
Choi and Garcia 2001).

Table 5 .
Optimal numerical aspects for predicting density currents

Table 6 .
Minimum dilution comparison for the last section of the studied density currents T6:1 Cases Modified parameter S minF E S minF N S minF E =S minF Eb S minF N =S minF Nb results.We recognize that applying the κ-ε 882 turbulence model, which solves two more equations, is demand-883 ing in terms of computation time.Both the mixing length model 884 and the constant model can generate good approximations 885 within a more reasonable timeframe for field applications.886 The results obtained through this study show that by applying 887 the previous guidelines, 3D hydrodynamic models can reproduce 888 density current flows in stagnant receiving waters with errors of 889 less than 1.3% for a minimum dilution line and of 6% for a maxi-