CFD Model of the Density-Driven Bidirectional Flows through the West Crack Breach in the Great Salt Lake Causeway

Abstract

Stratified flows and the resulting density-driven currents occur in the natural environment and commonly in saline lakes. In the Great Salt Lake, Utah, USA, the northern and southern portions of the lake are divided by an east-to-west railroad causeway that disrupts natural lake currents and significantly increases salt concentrations in the northern section. To support management efforts focused on addressing rising environmental and economic concerns associated with varied saltwater densities throughout the lake, the causeway was recently modified to include a new breach. The purpose of this new breach is to enhance salt exchange between the northern and southern sections of the lake. Since construction, it typically exhibits a strong density-driven bidirectional flow pattern, but estimating flows and salt exchange has proven to be difficult. To obtain much needed insights into the ability of this hydraulic structure to exchange water and salt between the two sections of the lake, a field campaign coupled with CFD modeling was undertaken. Results from this study indicate that the vertical velocity profile in the breach is sensitive to density differences between flow layers along with breach geometry and water surface elevations. The CFD model was able to accurately represent the bidirectional flows through the breach and provides for improved estimates of water and salt exchanges between the north and south sections of the lake.

1. Introduction

Saline lakes represent 23 percent by area, and 44 percent by volume, of all lakes on Earth [1]. Saline lakes can be found in diverse environments, though most of them are terminal lakes and located in places with an arid climate. Some of the major saline lakes are the Caspian Sea (coastline shared by Azerbaijan, Iran, Kazakhstan, Russia, and Turkmenistan), the Great Salt Lake (GSL) (Unites States), and Lake Urmia (Iran). The importance of these lakes is indicated by their proximity to urban centers (e.g., GSL is near Salt Lake City) and their thriving shipping, fishing, mineral, and tourism industries. In addition to direct economic benefits, saline lakes and the surrounding wetlands provide unique habitats for birds and other fauna [2,3]. Climate-change propelled droughts and increases in human consumption have resulted in lowering water levels and volumes of these saline lakes across the world [4]. One of the first recorded incidences of human consumption leading to saline lake desiccation occurred in the Tarim Basin, leading to the downfall of the Loulon Kingdom in 645 CE [5]. Some examples of saline lake desiccation due to over consumption in the last 100 years include Lake Urmia in Iran [6] and Owens Lake in the United States that was entirely dried up in 1940. Desiccation of saline lakes and the resulting alarming decline in related ecosystems have been documented in the past [4,7,8]; this has led to a renewed interest in preserving these systems as a subset of these lakes have seen calamitous decline in lake levels in the last 20 years [9]. The GSL is facing similar issues with the lowest levels being recorded in summer 2021 due to over consumption throughout the watershed. Potential desiccation and related ecosystem issues have been recently highlighted and are a major concern [4].

As the saline lake levels drop due to climate change and over consumption, there is a clear need to understand changes in water and salt transport throughout these lakes. This is of even greater importance where impoundments or structures alter or restrict lake water and salt movement. For example, in the GSL a railroad causeway divides the lake into northern and southern sections and limited connectivity is provided at causeway breaches. In these settings, the water and salt exchanges through these hydraulic structures are difficult to quantify, but critical for management efforts related to long-term ecosystem health and related industries via the evolution of the GSL’s density profile [10,11].

To provide insight regarding the water and salt exchanges needed for GSL management efforts, a field campaign coupled with a computational fluid dynamics (CFD) modeling of the density-driven flows through the recently commissioned West Crack (WC) Breach was undertaken. The modeling effort provides a better understanding of the bidirectional buoyancy-driven flow through the WC breach that provides GSL management and stakeholders with insights into how the openings in the causeway could be used for effective brine management in the future [12,13].

2. Great Salt Lake and the West Crack Breach

The GSL is a significant resource to Utah, contributing more than $1.3 billion annually to Utah’s economy through mineral extraction and its world-leading brine shrimp production [14,15]. Recreational attractions connected to the GSL also play a significant part in the local economy, primarily due to the GSL’s unique ecosystem that provides habitat for almost five million migratory birds every year [16]. As a result, significant research and monitoring efforts have been undertaken for decades [17], especially since the construction of the Union Pacific Railroad Company (UPRR) causeway in 1959. This railroad causeway divides approximately 1/3 of the lake from the southern section (see Figure 1) and has altered its characteristics including natural circulation and the water and salt balances [18]. For example, since 1959 the northern part of the GSL has become more saline while the salt content of the southern part with significant freshwater inflows has reduced [17,19]. The density, ρ, of North saltwaterhas on-an-average been ~100 to 150 kg/m³ higher than the South from 1966 to 2019, except during high lake levels [20]. Two culverts were built in the causeway in 1959 to facilitate water connectivity between the northern and southern sections. Subsequently, a breach was created in 1984 near the western end of the causeway. The original culverts were decommissioned in December 2013 with a second WC Breach (see Figure 2) opened in December 2016. Over the past 3 years, the water surface elevation (WSE) of the northern section of the lake has dropped significantly, illustrating the importance of the flow though the breach on both water and salt balances of the different sections of the GSL.

The difference in density, primarily due to varied salt concentrations, between the northern and southern sections of the GSL results in the formation of density-difference driven flow through the WC breach. In general, the WSE of the southern section is always higher than the north section. If there was no difference in density between the two sides, the flow through the breach would consistently be from South to North. However, the higher density on the North results in south water floating over the north water and the north water plunging in order to flow south. The result is bidirectional flow through the breach, where in the upper portion of the water column typically flows northward and the lower portion, with higher density, flows south. Consequently, this flow results in formation of a stable and stratified, but hydraulically complex, deep brine layer (DBL) in the southern part of the GSL. The connection between the density gradient-driven gravity current flowing through the breach and the formation of DBL is succinctly illustrated by the disappearance of the DBL during a three-year period when the last culvert was closed in December 2013 and a new breach was opened in December 2016 [22]. Although the DBL may be influenced by other factors, analyses had shown strong sensitivity of deeper salinity concentrations in the southern section of the GSL to the geometry and adjacent submerged berms of the new breach [12,13]. This essentially points to the fact that the new WC Breach is the location where most of the exchange of water and salt occurs between the northern and southern section of the lake. Thus, to correctly estimate the water and salt mass balance between two parts of the lake it is imperative to model the flow and salt transport through the WC Breach accurately.

Prior to the WC Breach several modeling efforts were undertaken. One of the first dedicated hydraulics models built to estimate the flow going through the culverts in the causeway was developed by Holley and Waddell [23], and later advanced by Loving et al. [10] to estimate the GSL saltwater densities and bidirectional flow through the causeway. While this model is capable of calculating the bidirectional flow through the causeway openings, the calculated flows were not in good agreement with measured flows, especially with the creation of the new WC breach. This is not completely unexpected due to the fact that Holley and Waddell’s model was developed for bidirectional flow through culverts, which are significantly different geometrically from the bridge built over the breach (see Figure 2). The newer breach has components (e.g., piers) that are expected to make the flow highly three-dimensional (3D). Thus, the complex bidirectional flow though the WC breach warranted the use of a CFD modeling approach to accurately capture the dynamics. In spite of the flow in the GSL being highly three-dimensional and complex, there has been very few studies that have used 3D hydrodynamic modeling to study different aspects of the lake. Spall [24,25] used 3D hydrodynamic models that incorporated baroclinic and barotropic responses, tidal forcing, wind stresses, Coriolis effects, surface thermal forcing, inflows, outflows, and transport of salt, heat, and passive scalars. Spall investigated the effects of surface heat flux and wind forcing on the spatio-temporal variations in the flow patterns within the southern section of the GSL. Neither of the studies modeled the spatio-temporal dynamics of the DBL in the GSL or the bidirectional flow passing through the causeway openings. To address these modeling gaps, a CFD model is required that captures the complex interaction between the bidirectional flow and the hydraulic structure that uses a multiphase CFD solver that can track the interface (fluid–fluid or fluid–air) between different phases. These solvers have been used successfully to study density flows and various complex multiphase flow problems (e.g., multiphase sewage flows in water reclamation plants [26] or estimation of flow discharge through hydraulic structures [27]).

3. Materials and Methods

3.1. Field Campaign

Field measurements were conducted at the WC Breach of the Great Salt Lake causeway, located approximately 53 km northwest of Salt Lake City and about 27 km west of Promontory Point. At this location lake depths were on average less than about 6 m. Field data were collected from 1 October 2020, to June of 2021 with supplemental data obtained from the U.S. Geologic Survey (USGS). USGS field data have been collected from the time the WC Breach was completed (Dec. 2016). Field observations were made as part of this modelling effort via shore access and by kayak (see Figure 3) that focused on specific conductance (σ), total dissolved solids (TDS), and water temperature measurements at and adjacent to the breach.

To better understand the relationship between σ, TDS, and ρ, water samples were collected at multiple locations north and south of the breach at specific flow depths (e.g., 0.1 m to 4.6 m at approximately 0.3 m spacing). Locations were GPS located and depths were measured with a field tape accurate to ±30 mm. A simple tubing and syringe system was used to collect samples, with the tubing placed at the prescribed depth and then a volume secured in the large syringe. Flushing was performed between each sample to minimize contamination and ensure accurate sampling. Samples were refrigerated and analyzed at the Utah Water Research Laboratory for TDS using standard methods for the analysis of water and wastewater [28]. TDS was then used to compute ρ for CFD modeling. At these same sampling locations, σ and water temperature measurements were made using an In-Situ AquaTroll 600 multi-parameter water quality sonde with reported instrument accuracy of ±2% and ±0.1 °C, respectively, for readings from 100,000 μS/cm to 200,000 μS/cm [29]. These sensors were able to be calibrated within 5% at 20 °C. Additional point vertical profiles were made with the sonde at locations approximately within 1 km north and 1 km south of the breach and approximately 2 km east-to-west. These general measurements indicated generally uniform σ values adjacent to the density-driven currents through the breach. Therefore, it was appropriate to install longer term monitoring locations approximately 400 m from the breach on the south bank of the causeway to capture a representative concentration of south water (Figure 3a) and on a pier approximately 10 m from shore and 800-m west of the WC Breach to capture north water concentrations. Instruments were placed within a perforated PVC casing that allowed for natural water circulation. The southern instrumentation casing extended out into the lake and was approximately 1/3 the depth from the lake bottom. These longer-term north and south sonde measurements were initially collected at 30 min intervals, but after November 2021, sampling rates were increased to 10 min intervals. Data were downloaded from the sondes every 1–2-months with routine calibrations performed each visit. The temporal σ data were used to consider a range of ρ for sensitivity analyses in the CFD.

In addition to the field data collected as part of this study, the USGS and other government agencies have GSL data extending over more than 50 years. The USGS data includes various buoys measuring σ at monthly intervals located in the vicinity of the causeway (Figure 1 and Figure 4), but at a considerable distance from the breach. However, the USGS has collected samples at the WC breach at monthly intervals (locations 10010026 and 10010025) since 2017. They also have an acoustic doppler velocimeter meter (ADVM) uplooker at the WC breach that is adjacent to a row of piers, directly below the north face of the bridge decking (see Figure 4). The ADVM sensor is approximately 0.5 m above the lake bottom and begins sampling 1.3 m above the lake bottom. The USGS estimates velocity profiles at this location using the Index Velocity Method [30] with corresponding index rating and 9 vertical bins. During this study ADVM data were shared with USU, but the official review and online publication of the data by the USGS was underway and had not yet been completed. Due to uncertainty in the specific location of the ADVM relative to the pier in the transverse direction, the CFD portion of this study considered vertical profiles at 10 distances from the respective piers selected based upon information gathered from the USGS (Figure 4b). In addition to salt concentrations and velocity measurements, the USGS collects WSE data at an adjacent pier in the northern section (See Figure 3b), at a pier to the east of the WC breach in the southern section, and at the breach (Figure 2a).

The CFD model required density inputs for the boundary conditions on the north and south sections of the lake. Field instrumentation provided σ and water temperature data to estimate water densities, ρ, using linear interpolation of the USGS data [31], data provided to the USGS from monitoring during and after construction [32,33], and data collected as part of this modeling effort. In addition to σ, the TDS was also referenced for selecting density, with a conversion from TDS to density via a conversion equation employed by the USGS [17]. The final range of north and south densities selected for CFD testing was based upon the continuous σ measurements, TDS measurements taken for this study (USU UWRL), and the monthly USGS TDS measurements.

3.2. Numerical Model

3.2.1. Model Overview and Flow Equations

FLOW-3D^TM by Flow Science, Inc.© was used to numerically simulate the density-driven flow through the WC breach for the typical bidirectional case (surface flows to the north with subsurface flows to the south). This solver uses the Navier–Stokes equations with a finite v-volume method with the free-surface (air–water boundary) resolved spatially and temporally via the FLOW-3D Volume-of-Fluid (TruVOF) method, [34]. The CFD model development included an initial phase where different aspects of the model were tested in a sectional slice (2D setup) with the spanwise direction represented using a single cell thick layer of elements. This was followed by the main phase in which simulations were conducted with the computational domain covering the entire breach and adjacent portions of the lake and breach (3D setup). The computational element sizes used in the 2D and 3D setup simulations were the same. The 2D simulations were primarily used to ascertain the most appropriate turbulence model to be used. FLOW-3D has multiple turbulence closures, including one that uses a large eddy simulation (LES) like approach to calculate the turbulent eddy viscosity. The Reynolds-averaged Navier–Stokes equations under the Boussinesq approximation are solved along with the conservation of mass, and a scalar/density transport equation that models the transport of salt. The equations solved are:

$$\frac{\partial u_i}{\partial t}+\frac{1}{V_F}\left(u_j{A}_j\frac{\partial u_i}{\partial x_j}\right)=-\frac{1}{\rho }\frac{\partial p}{\partial x_i}+\frac{\left(\rho -{\rho }_S\right)}{{\rho }_S}{G}_i+f_i$$

(1)

$$\frac{\partial }{\partial x_i}\left(u_i{A}_i\right)=0$$

(2)

$$\frac{\partial \rho }{\partial t}+\frac{1}{V_F}\left(A_i\frac{\partial \rho u_i}{\partial x_i}\right)=\frac{1}{V_F}\frac{\partial }{\partial x_i}\left(K{A}_i\frac{\partial \rho }{\partial x_i}\right)$$

(3)

where V_F is the fractional volume of flow, x_i represents the selected cartesian coordinate, G_i and f_i are the body and viscous accelerations, respectively. A_i is the fractional area open to flow in a particular direction, u_i is the mean velocity, p is the pressure, $\rho$ is the density of the flow, ${\rho }_S$ is the reference density (which is the same as the density of south side of the lake), and K is summation of molecular and turbulent diffusivity of the transported salt. In Equation (1), f_i is defined as:

$$\rho V_F{f}_i={\tau }_{b,i}-\frac{\partial A_j{S}_{ij}}{\partial x_j}$$

(4)

$$S_{ii}=-2{\mu }_{tot}\left[\frac{\partial u_i}{\partial x_i}\right], S_{ij}=-2{\mu }_{tot}\left[\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right],$$

(5)

S_ij is the strain rate tensor, ${\tau }_{b,i}$ is the wall shear stress, ${\mu }_{tot}$ is the total dynamic viscosity, which encapsulates both the molecular and turbulent eddy viscosity (${\mu }_T$). The differences between all the turbulence models tested for the study is due to the different approaches adopted to caclulate the turbulent eddy viscosity (${\mu }_T=\rho {\nu }_T$) [35].

This study tested four turbulence schemes: the k-ε, the RNG k-ε, k-ω, and LES. Equations (6)–(8) are used to determine k and ε in the standard k-ε model.

$$\frac{\partial k}{\partial t}+u_i·\frac{\partial k}{\partial x_i}=\frac{\partial }{\partial x_i}\left(\frac{{\nu }_t}{{\sigma }_k}\frac{\partial k}{\partial x_i}\right)+{\nu }_t\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)\frac{\partial u_i}{\partial x_j}-\epsilon$$

(6)

$$\frac{\partial \epsilon }{\partial t}+u_i·\frac{\partial k}{\partial x_i}=\frac{\partial }{\partial x_i}\left(\frac{{\nu }_t}{{\sigma }_{\epsilon }}\frac{\partial \epsilon }{\partial x_i}\right)+C_1\frac{\epsilon }{k}{\nu }_t\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)\frac{\partial u_i}{\partial x_j}-C_2\frac{{\epsilon }^2}{k}$$

(7)

$${\nu }_t=C_v\frac{k^2}{\epsilon }$$

(8)

with C₁ = 1.44, C₂ = 1.92, C_ν = 0.09, and the Prandtl number of turbulence for ε and k are σ_ε = 1.3 and σ_k = 1.0, respectively. The RNG k-ε model (Equations (9) and (10)) with turbulence viscosity determined via Equation (11) and μ = νρ:

$$\frac{\partial k}{\partial t}\left(\rho k\right)+\frac{\partial }{\partial x_i}\left(\rho k{u}_i\right)=\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\frac{\partial k}{\partial x_j}\right]+P_k-\rho \epsilon$$

(9)

$$\frac{\partial }{\partial t}\left(\rho \epsilon \right)+\frac{\partial }{\partial x_j}\left(\rho \mathsf{\epsilon }{u}_i\right)=\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_{\epsilon }}\right)\frac{\partial \epsilon }{\partial x_j}\right]+C_{1\epsilon }\frac{\epsilon }{k}{P}_k-C_{2\epsilon }\rho \frac{{\epsilon }^2}{k}$$

(10)

$${\mu }_t=\rho C_{\mu }\frac{k^2}{\epsilon }$$

(11)

where k is the turbulent kinetic energy (TKE), ε the rate of dissipation of k, μ dynamic viscosity, μ_t turbulent dynamic viscosity, and P_k is the production of k, and C_1ε = 1.39, C_2ε = 1.42 and C_μ = 0.085. The k-ω model is Equations (12)–(14):

$$\frac{\partial k}{\partial t}+u_i·\frac{\partial k}{\partial x_i}={\tau }_{ij}\frac{\partial u_i}{\partial x_i}+\frac{\partial }{\partial x_i}\left[\left(\nu +\frac{{\nu }_t}{{\sigma }_k}\right)\frac{\partial k}{\partial x_j}\right]+\beta k\omega$$

(12)

$$\frac{\partial \omega }{\partial t}+u_i·\frac{\partial \omega }{\partial x_i}=\frac{\partial }{\partial x_i}\left[\left(\nu +\frac{{\nu }_t}{{\sigma }_{\omega }}\right)\frac{\partial \omega }{\partial x_j}\right]-{\beta }^{*}{\omega }^2+\alpha \frac{\omega }{k}{\nu }_t\left({\tau }_{ij}\frac{\partial u_i}{\partial x_j}\right)$$

(13)

$${\nu }_t=\frac{k}{\omega }$$

(14)

where σ_ω = 0.5 and σ_k = 0.5, α = 0.52, β = 0.072f_B and β^* = 0.09f^*_B. Finally, the LES solver used herein is the Smagorinsky model [36,37] that computes all turbulent flow structures that can be resolved by the computational grid. Sub-grid turbulence effects are modeled through the turbulent eddy viscosity (${\nu }_T$), calculated using Equations (15)–(17):

$${\nu }_t={\left(cl\right)}^{2}S$$

(15)

$$l={\left(\Delta x\Delta y\Delta z\right)}^{1/3}$$

(16)

$$S=\sqrt{2S_{ij}{S}_{ij}}$$

(17)

where c is the Smagorinsky coefficient defined as 0.1 [38], l is the sub-grid length scale and S is the strain rate tensor, and $S_{ij}$ is the strain rate tensor and is defined as $S_{ij}=\frac{1}{2}\left[\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right]$. $\Delta x$ is the size of the element in the x-direction.

The aforementioned turbulence closure schemes were tested during the 2D modeling phases, and the best performance (relative to field data) was given by the LES and the k-ω models. Interestingly, use of the LES model (compared to k-ω) resulted in a reduction in wall-clock time by about 30% for 2D simulations. Due to these factors, the turbulence closure used for the final 3D phase of modeling was LES. The turbulence closure models did not have a term explicitly to account for stratification driven turbulence dissipation. Though requirement of that component in the models has been found to be more important for modeling turbulence in a stably-stratified setup [39]. Previous studies have shown that accounting for stratification effects in turbulence closures have less significance for modeling accelerating gravity current dynamics [35,40].

The first-order momentum advection scheme and first-order volume-of-fluid advection were used with an implicit pressure solver with the generalized minimum residual (GMRES) method (no over- or under- relaxation) with the subspace size of 15, inter-block boundary coefficient of 0.25, and volume fraction cleanup of 0.05. The bathymetry, causeway, and piers were numerically represented as solids by the Fractional Area/Volume Obstacle Representation (FAVOR) method. In order to reduce numerical diffusion and increase accuracy in solving the density transport equation, the second-order monotonicity preserving scheme was used. Other lower-order options (e.g., first-order approximation) were also tested and were found to give erroneous results.

3.2.2. Model Domain and Boundary Conditions

The 3D numerical domain was 295 m by 275 m centered on the breach (Figure 5). Hexahedral cells of different dimensions were employed to resolve different parts of the domain (Figure 6a). Grid convergence [41] was studied for both 2D and 3D models, where the grid convergence for the 2D model included eight cell sizes refinements ranging from 0.3 m to 0.03 m, and the grid convergence for the 3D model included three cell size refinements ranging from 0.3 m to 0.08 m. Results from the grid-convergence indicated solution independence was achieved for cells of 0.15 m (see Figure 6b). For this study the final discretization of the main phase domain was comprised of two mesh blocks of 0.3 m cubic cells linked in series, with two nested blocks with 2:1 cell refinement at the boundary (i.e., the second nested block cell size was x = y = z = 0.08 m). The number of cells in the final 3D simulations were 23.4 million.

The northern boundary conditions were assigned a WSE_N of 1277.76 m and density of ρ_n = 1180.22 kg/m³, based upon field data (noted in red region as P in Figure 5). Similarly, the southern boundary conditions were assigned a WSE_S of 1277.82 m, or 0.06 m higher with density of ρ_s =1097.75 kg/m³ (noted as P in blue region in Figure 5). The boundaries where WSE is assigned are specified as pressure boundaries of the type that allows free inflow/outflow of water and corresponding flow properties, including variable water densities, are tracked and corresponding advective scalars are prescribed at the boundary along with a hydrostatic pressure distribution. The top boundary was defined as a Dirichlet boundary condition with atmospheric pressure set to 0 or the gage pressure. Note that diffusion at the fluid–no fluid interface is managed by a fluid-volume correction (<±0.01%). The wall boundary condition (bottom, all solids, etc.) were prescribed as a no-slip surface condition where the solver checks for laminar or turbulent flow conditions in a cell at the boundary and then prescribes an appropriate log-velocity profile and a zero tangential velocity at the wall. FLOW-3D also considers smooth vs. rough wall logarithmic expressions via a modified viscosity that reflects enhanced momentum exchange due to a roughness height:

$$u_{tan}=u_s\left[\frac{1}{\kappa }\ln \left(\frac{\rho u_{s}y}{\mu +\rho a{u}_s{k}_s}\right)+5\right]$$

(18)

where u_tan is the tangential velocity at the wall at a distance from the wall y, κ is the von Karman coefficient, u_s is shear velocity, a is a constant equal to 0.246 (rough wall boundary), and k_s is a user-defined roughness height. Due to the macro-roughness features through the breach, k_s is not significant and with a negligible roughness height, good agreement was obtained in the velocity profile discussed in Section 4. Note that local velocities are used as expected to compute wall shear stresses and are included as required as transport parameters in corresponding turbulence models. Based upon the final selected mesh size y+ values close to the boundary are about 5000, which is much higher than those used by boundary-resolved LES simulations where y+ < 10.

Initial boundary conditions were defined for each main mesh block with fixed water surface elevations and water densities corresponding to the respective northern and southern boundaries. All the simulations were monitored and run long enough to reach a steady state for the given boundary conditions.

4. Results and Discussion

4.1. Density and Water Surface Elevations

As previously noted, USGS has multiple buoys in the GSL where grab samples are analyzed for TDS along with corresponding σ. USGS calculates ρ as a function of TDS via Equation (19) from [17]:

$$TDS=\frac{\left(\rho -1\right)1000}{0.63}$$

(19)

where TDS is in g/mL and ρ is in g/L at 20 °C. Therefore, to directly compare ρ collected from this study (identified in Figure 7 as UWRL data) to USGS data the same technique was adopted. Naftz et al. [42] developed an equation of state for the hypersaline water from the southern part of GSL, where ρ can be computed as a function of temperature and salnity. The Naftz et al.’s equation of state was not useful for the current study due to lack of a relationship between specific conductance and salinity, and also the fact that the aforementioned equation was developed only using data from the south part of GSL. Note that since 2016, a private consultancy HDR was also monitoring ρ at the new breach on behalf of the Union Pacific Railroad and they also adopted the USGS approach; these three data sets are presented in Figure 7.

These data in addition to the σ measurements made by the authors served as the basis of choosing ρ_S and ρ_N used in the model simulations. Note that the CFD simulations require that ρ be in units of kg/m³, thus the density difference Δρ is provided accordingly.

In addition to Δρ, the bidirectional flow patterns in the WC Breach are also influenced by water surface elevations; the continuous measurements of USGS for north and south water surface elevations near the breach for the entire period since the WC Breach was completed are also shown (Figure 8). From this data, a specific day corresponding to ADVM measurements and TDS measurements was selected for modeling (

Table 1. Density difference Δρ and the difference in WSE imposed for the 3D simulations. A range of Δρ was tested due to field data and corresponding uncertainty in estimation of ρ_N.

Sim (#)	ρS (g/cm3)	ρN (g/cm3)	Δρ (kg/m3)	ΔWSE (m)
1	1.0978	1.1699	72.1	0.061
2	1.0978	1.1802	82.4	0.061
3	1.0998	1.1957	95.9	0.061
4	1.0978	1.2008	103.0	0.061
5	1.09782	1.2127	115.5	0.061

).

Some variability in field measurements is likely due to instrument accuracy (Figure 7 and Figure 8) and also due to lake dynamics and temporal variations in response to climate, weather patterns, etc. These factors were considered when selecting boundary conditions for the CFD simulations. A case that is representative of the density driven bidirectional flow patterns was selected. Specific velocity data and corresponding WSEs were used in the modeling and a sensitivity analysis of densities and velocity profile location was completed.

4.2. Model Comparison with Field Measurements

Once the mesh resolution study was completed, the accuracy of the 3D model was tested against the unpublished USGS velocity profile measured at the WC breach (Figure 6). Several unknown parameters needed to be considered while comparing the CFD model results with the field measurements. First, the exact density difference Δρ between the northern and southern sections of the lake is uncertain. As described in Section 3.1 and Figure 4, the field data, including the time-series of σ, indicated that there was a small range of possible ρ values for the northern part of the lake. Thus, for the 5 cases considered below, different densities between the northern and southern sections of the lake were simulated (Table 1).

The second source of uncertainty was the exact location of the ADVM uplooker with respect to the central column of piers at the WC breach. It is known that the uplooker is located slightly north of the northern-most pier in the center column, and within 1.5–2.5 m east of the that pier. However, to ascertain the potential location of the uplooker, the streamwise velocity was probed at multiple locations east of the pier (Figure 4) and the locations were investigated using different density differences (Figure 9). The results show that the correct density difference, Δρ, between the northern and southern sides is ~82.4 kg/m³, and the location of the uplooker is about 1.5 m east of the pier. Given the uncertainty in the field data and its location, the comparison between the measured and computed velocity (V) profiles is satisfactory (Figure 10). The only points in the profile that do not compare well, which have more than 50 percent error, are V at the interface of the north–south and south–north flows. This is not unexpected, as there is some uncertainty in the field velocity profile at this interface and this section is more challenging to simulate numerically as there is a sudden change in density and velocity.

4.3. Simulated Flow Field

Based on the field observations (Figure 11) and field data, the general dynamics of the flow is clear; hypersaline water from the north side of the lake plunges below the south side water that is flowing south to north (S-N) just before entering the breach. This creates a relatively thick and denser bottom current flowing north to south (N-S), and a relatively thinner and less dense S-N current flowing at the surface. As this top current exits the breach it spreads longitudinally and laterally, and thins as it travels north, with its thickness decreasing with distance away from the breach. To provide a general overview of the hydrodynamics at the breach, ρ and the S-N (y-direction) velocity V has been plotted for the vertical-plane orthogonal to the x-axis and going through the center of the breach (Figure 12 and Figure 13). CFD results clearly capture the dynamics observed in the field and the simulation results provide more insights about the formation of the density driven bidirectional flow pattern (N-S and S-N currents). In agreement with field observations, the thickness of the distinct layers of the flow are similar with the bottom layer being relatively thick with a non-negligible velocity. The hypersaline water from the north maintains its thickness some distance from the causeway, eventually forming the DBL. The hypersaline northern current also appears to go through a subsurface hydraulic jump just after it plunges and before entering the central part of the breach, which is a function of the submerged dike or berm specifically placed to manage flows through the WC breach. The hypersaline N-S current can be seen to accelerate just after plunging under the incoming S-N flow. The S-N flow also accelerates while going through the breach and can be seen to travel much further into the northern side than the density plots indicate. In order to get a clearer view of the flow evolution, vertical planes orthogonal to the S-N (y-direction) direction have been plotted at 12-m intervals starting south of the breach (Figure 14).

It can be observed that the flow has a clear bidirectional form, though the thickness of the N-S bottom current and the S-N top layer flow change spatially (Figure 14). It was observed that a distinct N-S flowing bottom layer, which is still relatively thick, is primarily confined to the center of the channel at location 1 (Figure 14 top). The flow has expanded laterally, but is still mostly flowing along the channel that was dredged during breach construction. Charting the evolution of the N-S current, starting from the northernmost slice at location 7 (Figure 14 bottom), it is observed that the current accelerates and decreases in thickness as it approaches the breach. The current accelerates from about 0.2 m/s to about 1 m/s in a span of about 24 m. The approximate densimetric Froude number of the bottom current at the location of the uplooker measurement is about 0.7. Thus, there is a chance that the flow might become supercritical while accelerating towards the breach, justifying that hydraulic jump is plausible. Additionally, the piers of the breach can be observed to have induced mixing and decrease the density stratification. This can be inferred from the sudden increase in bottom layer thickness going from location 5 to location 4 (Figure 14). The effect of the piers on the structure of the N-S current can be observed 12 m downstream (see Figure 14, location 3). Going further south, the bottom current further decreases in thickness (see Figure 12), but retains velocities around 0.8–0.9 m/s. This N-S current transforms to become the dense DBL.

The S-N flow, which is confined to the top layer of the water column, can be observed to be accelerating while approaching the breach. While flowing among the piers, the S-N flow covers more area in the cross-section than the N-S flow. Though by the time the flow passes through the breach (Figure 14, location 5), the water column has been divided into two distinct parts. Within 24 m, the S-N current gets confined within 1/10th of the water column, while accelerating to about 0.6 m/s. This S-N current now has a densimetric Froude number of about 1.2 and is supercritical.

In order to further understand what the flow looks like at different elevations, density and velocity magnitudes at four different layers have been plotted in Figure 15. The flow at the surface is primarily S-N, and the current can be seen to be expanding radially away from the WC breach. This matches with the field observations, where this current confined to the top-layer of the water column was observed beyond 1000 m north of the WC breach. The next layer at 1.1 m from the surface, the flow has relatively high velocity while going through the breach, especially once the S-N flow has crossed the piers and is pushed further up in the water-column. The next layer at 2.5 m below the surface, the flow is dominated by the N-S current, though there are parts where this layer is near/passing-through the interface between the two currents. The bottom current can be observed passing through the breach and continuing along the central-channel (Figure 15g), while a portion travels laterally. This plunging of the N-S current might lead to overturning of the flow, leading to enhanced mixing near the edges of the breach on the south side (Figure 14 location 3). The N-S current also seem to go through relatively high acceleration after entering the breach (Figure 14, location 6), and this is also the zone where the N-S current plunges towards the bottom. The layer at about 3.4 m from the surface (Figure 15g,h) is dominated by the bottom hugging N-S current. High velocities are observed within the dredged channel.

These modeling results also highlight the role of bottom topography in the bottom current behavior. The distinct berm (bottom sill) is observed North of the breach, which can be used to reduce the N-S flow, and consequently the DBL. Though how effective it will be, or how high it has to be made to impact the DBL is not well understood. It can also be observed that the bottom of the lake on the southern side increases in elevation just after the dredged region of the breach. This will have a significant effect on the mechanism of the DBL formation, and will be studied in detail in the future. This CFD model can also be modified and used to study the effect of the berm.

For the purpose of developing an accurate water and salt balance model of the lake, it is important to accurately estimate the north and south flows through the breach. The results of this study indicate that the CFD model can be used to simulate other lake levels and water density variability that could be combined to develop a more robust rating curve for the breach. Furthermore, this model can be a used as a forecasting tool to anticipate future conditions as lake levels continue decreasing to unprecedented conditions resulting from extended periods of drought.

5. Conclusions

The density-driven bidirectional current through the recently completed West Crack breach at the Great Salt Lake causeway was simulated with an LES CFD model to support management efforts that include environmental and economic concerns. Modeling was based upon field campaigns and datasets provided by the USGS. Results from this study indicate that the vertical velocity profile in the breach had good agreement with unpublished USGS velocity profile data with an R² = 0.9578 and generally within 20 to 50% at different depths throughout the water column in the breach. The model is sensitive to density differences between flow layers along with breach geometry and water surface elevations, which is also the case for lake flows through the WC breach. The model resulting from this work can be used to better approximate water and salt exchanges between the north and south sides of the GSL, provide insight into future management of the lake, and be applied as a forecasting tool to anticipate future lake conditions.