Mixing efficiency in large-eddy simulations of ... - Princeton University

1 downloads 0 Views 378KB Size Report
buoyancy scale Lb on mixing efficiency in LES of stratified turbulence. In addition ...... large-eddy simulation: application to a neutral atmospheric boundary layer.
Under consideration for publication in J. Fluid Mech.

1

Mixing efficiency in large-eddy simulations of stratified turbulence Sina Khani† Program in Atmospheric and Oceanic Sciences, Princeton University, Princeton, NJ 08544, USA (Received ?; revised ?; accepted ?.)

The irreversible mixing efficiency is studied using large-eddy simulations (LES) of stratified turbulence, where three different subgrid-scale (SGS) parameterizations are employed. For comparison, direct numerical simulations (DNS) and hyperviscosity simulations are also performed. In the regime of stratified turbulence where F rv ∼ 1, the irreversible mixing efficiency γi in LES scales like 1/ (1 + 2P rt ), where F rv and P rt are the vertical Froude number and turbulent Prandtl number, respectively. Assuming a unit scaling coefficient and P rt = 1, γi goes to a constant value 1/3, in agreement with DNS results. In addition, our results show that the irreversible mixing efficiency in LES, while consistent with this prediction, depends on SGS parameterizations and the grid spacing ∆. Overall, LES approach can reproduce mixing efficiency results similar to those from DNS approach if ∆ . Lo , where Lo is the Ozmidov scale. In this situation, the computational costs of numerical simulations are significantly reduced because LES runs require much smaller computational resources in comparison with expensive DNS runs. Key words:

1. Introduction Large-eddy simulation (LES) is an approach that has been considered in numerical simulations of stratified turbulence (Kaltenback, Gerz & Schumann 1994; Armenio & Sarkar 2002; Remmler & Hickel 2012; Khani & Waite 2014, 2015; Schaefer-Rolffs 2016). The main goal of LES is to decrease the computational costs of direct numerical simulation (DNS), which seems to be an unnecessarily expensive numerical method because most of the computational resources are dedicated to resolving the non-energetic and viscous-affected eddies around the Kolmogorov scale Ld (see e.g. Pope 2000). In LES, we aim to directly resolve just the large energy-containing eddies, and parameterize the effects of the smaller-scale motions on resolved large eddies using the subgrid-scale (SGS) models. It has been known that, unlike isotropic turbulence, grid spacing (i.e. the smallest resolved scale) and SGS parameterizations can significantly influence on the dynamics of anisotropic turbulent flows (see for example, Port´e-Agel et al. 2000; Khani & Waite 2014, 2015, for LES studies in atmospheric boundary layer and stratified turbulence). Khani & Waite (2014, 2015) have shown the importance of resolving the buoyancy scale Lb ∼ urms /N in LES in order to capture the fundamental characteristics of stratified turbulence using isotropic SGS models. Here, urms and N are the root-mean-square (r.m.s.) velocity and buoyancy frequency, respectively (see below for the explicit definition of † Email address for correspondence: [email protected]

2

S. Khani

urms ). Also, they have shown that LES will fail if the buoyancy scale is not explicitly resolved (Khani & Waite 2014, 2015). Mixing efficiency is defined as the ratio of the buoyancy flux to the mean flow turbulence production. The vertical mixing is a key parameter for researchers in the atmospheric sciences and physical oceanography. For example, it has been shown that eddy mixing has an important influence on breaking gravity waves and vertical transport of heat and tracers in geophysical turbulence (see e.g. Lilly et al. 1974; Weinstock 1978; Osborn 1980). Mixing efficiency has been studied numerically in stratified- turbulence and shear layers (e.g. Caulfield & Peltier 2000; Shih et al. 2005; Salehipour & Peltier 2015; Maffioli et al. 2016), in which turbulent diffusivity, mixing coefficient and flux Richardson number Rif , are modelled or estimated using DNS data. Using high resolution DNS results, Maffioli et al. (2016) have argued that mixing efficiency approaches to a constant value when the Froude number is small. This finding is in line with DNS study of Venayagamoorthy & Stretch (2010), where they have proposed that mixing efficiency becomes constant when the gradient Richardson number Ri is large. On the other hand, Shih et al. (2005) have shown that in weakly stratified turbulence, mixing efficiency depends on the buoyancy Reynolds number Reb when the Froude number is & 1. As a result, there is a controversy on parameters that mixing efficiency depends on, in the DNS approach. The ability of LES to reproduce the correct mixing efficiency is a promising approach and active research avenue in numerical simulations of stratified turbulence, where in LES, eddy- viscosity and diffusivity are parameterized in terms of resolved variables using SGS models. For example, Armenio & Sarkar (2002) have investigated the influences of flux Richardson number and mixing efficiency in LES of stratified channel flow. They have shown that increased stratification reduces the vertical buoyancy exchanges and turbulence levels, leading to a decrease in mixing efficiency. In this paper, we shall study mixing efficiency in LES of forced stratified turbulence, where three different SGS parameterizations are considered: the Kraichnan (1976), Smagorinsky (1963) and dynamic Smagorinsky (Germano et al. 1991) models. We compare LES results with those from DNS and hyperviscosity simulations. First, we intend to present a mathematical framework to define SGS mixing efficiency in LES. Also, we intend to estimate mixing efficiency using resolved variables by performing an scale analysis on corresponding equations in LES. In particular, we are interested in studying the effects of capturing (or not) the buoyancy scale Lb on mixing efficiency in LES of stratified turbulence. In addition, we aim to assess the performance of different SGS models in capturing the correct mixing efficiency by comparing to DNS results. The review of literature on stratified turbulence and mixing efficiency, along with mathematical formulations are given in §2. Section 3 shows simulations set-up and methodology. Results and discussion are represented in §4, and conclusion is given in §5.

2. Background At sufficiently small scales in the atmosphere and ocean where the Rossby number is large (i.e. horizontal scales below O(100) km in the atmosphere and O(1) km in the ocean), fluid motions are strongly affected by the stable stratification, but only weakly affected by the Earth’s rotation (Nastrom & Gage 1985; Riley & Lelong 2000; Lindborg & Cho 2001; Riley & Lindborg 2008; Waite 2014). In this range of scales, turbulent vortices are elongated horizontally and squeezed vertically (so-called pancake like vortices in stratified turbulence, e.g. Waite & Bartello 2004; Lindborg 2006; Brethouwer et al. 2007). Therefore, we might use different lengthscales in the horizontal and vertical directions in

Mixing efficiency in large-eddy simulations of stratified turbulence

3

stratified turbulence as are Lh ∼ u3rms /ǫ and Lv ∼ urms /N (e.g. Lindborg 2006; Khani & Waite 2013), where ǫ is the kinetic energy dissipation rate, and Lh and Lv are the horizontal and vertical length scales, respectively. The Navier-Stokes equations under the Boussinesq approximation for resolved variables in LES of stratified turbulence, can be written as (e.g. Khani 2015) g 1 ∂u ¯ p − ρ¯ez + ∇ · τ + f¯, (2.1) +u ¯ · ∇¯ u = − ∇¯ ∂t ρ0 ρ0 ∇·u ¯ = 0, (2.2) ∂ ρ¯ dΦ +u ¯ · ∇¯ ρ+w ¯ = ∇ · h, (2.3) ∂t dz where u ¯ = (¯ u, v¯, w) ¯ is the resolved velocity vector, ρ¯ and p¯ are the resolved density and pressure perturbations, respectively, and f¯ is the forcing term; the molecular- viscosity ν and diffusivity κ are negligible in our LES work. It is assumed that the background density field Φ varies linearly with the vertical displacement z (i.e. Φ = −ρ0 N 2 z/g); g and ρ0 are gravity and reference density, respectively; and ez is the unique vector in the vertical direction. Noteworthy, unresolved stress tensor τ and density flux h are unknown and required to be modelled using SGS parameterizations. The transport equations for the kinetic and potential energy can be obtained by multiplying u ¯ and β 2 ρ¯, where β = g/ (ρ0 N ), to equations (2.1) and (2.3), respectively, as follows   g 1 ∂E ∂ ∂E − p¯u¯j + u¯j τij − ρ¯w = ¯ − τij s¯ij + f¯j u ¯j , (2.4) +u ¯j ∂t ∂xj ∂xj ρ0 ρ0  ∂ ∂ g ∂ ∂ ρ¯ (P) − ρ¯w ¯= , (2.5) (P) + u¯j β 2 hj ρ¯ − β 2 hj ∂t ∂xj ρ0 ∂xj ∂xj where the kinetic energy E = (1/2)¯ ui u¯i , potential energy P = (β 2 /2)¯ ρρ¯, and rateof-strain tensor s¯ij = (1/2) (∂ u ¯i /∂xj + ∂ u ¯j /∂xi ). We can follow a similar procedure to obtain the transport equations for the kinetic and potential energy in wavenumber (Fourier) space if we write the equations of motions (2.1-2.3) in wavenumber space (see e.g. Khani & Waite 2014; Khani 2015, for more information). Osborn (1980) has proposed an equation for mixing efficiency based on the mixing coefficient Γ = B/ǫ, given as γr = B/ (B + ǫ), where B = gh¯ ρwi/ρ ¯ 0 is the buoyancy flux (angle bracket denotes averaging), and γr is also called the flux Richardson number Rif . By assuming a critical mixing coefficient Γc = 0.22, we can get γr = 0.18, a value which Osborn (1980) obtained using ocean measurements. Similar definition also used lately by Shih et al. (2005) in DNS of sheared stratified turbulence, where they concluded that mixing efficiency depends on the buoyancy Reynolds number Reb . They have shown that the flux Richardson number Rif is small when Reb ≪ 1, it increases to an approximately constant value Rif ≈ 0.2 when 7 < Reb < 100, and Rif decreases as Reb gets larger (Shih et al. 2005). Beyond this classical definition, Caulfield & Peltier (2000) has proposed an equation for irreversible mixing efficiency based on the kinetic and potential energy dissipation rates, as follows ǫp , (2.6) γi = ǫ + ǫp where ǫ and ǫp are the (molecular) kinetic and potential dissipation rates, respectively. This latter definition is more desirable in the study of turbulent mixing because unlike the former one, equation (2.6) does not consider the reversible part of mixing that is included in buoyancy flux and the flux Richardson number Rif . This new definition has also recently been used in DNS studies of stratified turbulence (Venayagamoorthy &

4

S. Khani

Stretch 2010; Salehipour & Peltier 2015; Maffioli et al. 2016). For example, Maffioli et al. (2016) have shown that γi is constant when the Froude number is small. In addition, using some experimental and observational evidences, it has been suggested that mixing efficiency could reach up to values about 0.4 to 0.5 (see e.g. Strang & Fernando 2001; Pardyjak et al. 2002; Peltier & Caulfield 2003). For example, Mashayek et al. (2017) have shown that mixing efficiency could increase to be ≈ 1/3 in stratified regions of the bottom ocean near topographies. In this paper, we use equation (2.6) to calculate the irreversible mixing efficiency in LES of stratified turbulence. Also, we should emphasize that the kinetic and potential dissipation rates in LES are given by SGS parameterizations as ǫ = τij s¯ij and ǫp = β 2 hj (∂ ρ¯/∂xj ), where τij and hj are parameterized in terms of the resolved variables using SGS models. We use the notation γi for irreversible mixing efficiency in LES, DNS and hyperviscosity simulations. However, we should note that γi is indeed the ‘SGS’ mixing efficiency in the LES approach because the molecular viscosity (or hyperviscosity) is neglected. Nevertheless, on average and for stationary turbulence, SGS dissipation rates should be equal to molecular dissipation rates (see tables 1 & 2). 2.1. SGS parameterization Common SGS parameterizations in physical space are Smagorinsky (1963) and dynamic Smagorinsky (Germano et al. 1991) models, in which the SGS momentum stress τij and density flux hj are related to the resolved rate-of-strain tensor s¯ij and density gradient ∂ ρ¯/∂xj using the eddy- viscosity νr and diffusivity νκ coefficients, respectively, given by τij = 2νr s¯ij ,

h j = νκ

∂ ρ¯ , ∂xj

νr = cs ∆2 (2¯ sij s¯ij )1/2 ,

νκ =

νr , P rt

(2.7)

where cs , ∆ and P rt are the Smagorinsky coefficient, grid spacing and turbulent Prandtl number, respectively. In classical Smagorinsky (1963) model, cs = (0.17)2 , as is proposed by Lilly (1967). In the dynamic Smagorinsky model, however, cs is not constant and is a time- and space-dependent variable, which is calculated dynamically by employing ˜ > ∆ (see e.g. Khani & Waite 2015, for more information about a coarse-filter scale ∆ cs calculations). Similarly, we can propose a turbulent-viscosity model in wavenumber (Fourier) space, for which the spectral SGS kinetic and potential energy transfers are related to the kinetic and potential energy spectra at the cutoff wavenumber kc using spectral eddy- viscosity and diffusivity coefficients, respectively (see e.g. Kraichnan 1976, for more information about the Kraichnan eddy-viscosity model).

3. Methodology Idealized LES runs of forced stratified turbulence in a cubic domain of side length L = 2π are considered. Random AR(1) forcing is applied to the rotational part of the horizontal velocity (vortical modes). The forcing is uncorrelated in space, but correlated in time with a correlation timescale of 10 timesteps ∆t (as in e.g. Waite & Bartello 2004; Khani & Waite 2014, 2015). The grid spacing ∆ is related to the cutoff wavenumber kc , given by ∆ = π/kc . The spectral transform method is considered for discretization of spatial derivatives, where the two-thirds rule (Orszag 1971) is employed in each direction to eliminate aliasing errors. Hence, the effective grid spacing is given by ∆ = 3L/(2n), where n is the number of grid points in the x, y and z directions. Also, the third-order Adams-Bashforth scheme is employed for time advancements. For comparison, we have also performed DNS and hyperviscosity simulations, in which the molecular viscosity ν and hyperviscosity ν4 are effective in the corresponding dissipation mechanisms (see Khani & Waite 2014, for more information about hyperviscosity simulations).

Mixing efficiency in large-eddy simulations of stratified turbulence

5

Simulations are initialized with random noise in low resolution (n = 256) hyperviscosity runs from time t = 0 to t = 300, and then spun up to run in LES, DNS and hyperviscosity simulations from time t = 300 to t = 450 (the time by which simulations results have already reached to statistically stationary state; similar approaches have been followed in Khani & Waite 2014, 2015). We consider simulations resolutions from n = 256 to n = 768, and the buoyancy frequency varies from N = 0.5 to 6 to ensure sufficiently small horizontal Froude number F rh = urms /N Lh . Results are averaged ˜ = 2∆, and turbulent over times 375 6 t 6 450, denoting by the angle bracket h·i, ∆ Prandtl number P rt = 1, which is consistent with the a priorip study in DNS of stratified turbulence (see e.g. Khani & Waite 2016). We use urms = hE(t)i since the vertical kinetic energy is much smaller than the horizontal; and Lb = 2πurms /N . Table 1 shows a list of parameters and averaged variables in LES runs, where the Smagorinsky, dynamic Smagorinsky and Kraichnan SGS models are considered. Also, simulation parameters and averaged variables for DNS and hyperviscosity simulations are shown in table 2, where kmax and kd are the maximum and Kolmogorov wavenumbers, respectively, and the ratio kmax /kd shows that the resolution of the Kolmogorov scale is reasonable for all DNS and hyperviscosity simulations (i.e. in all simulations, kmax /kd > 0.8, which is consistent with the criterion that is given by Moin & Mahesh 1998, for the resolution of 1/(6m−2) kd ). Here, kd = ǫ/ν 3 , where m = 1 and 4 in DNS and hyperviscosity simu1/2 lations, respectively. The Ozmidov scale Lo = 2π ǫ/N 3 is also shown in tables 1 & 2 for LES, DNS and hyperviscosity simulations. Moreover, we perform extra LES runs with different turbulent Prandtl numbers, ranging from 0.2 to 2, to provide a complete analysis on our scaling theory for mixing efficiency (see below).

4. Results and discussion 4.1. Kinetic energy and SGS dissipation spectra Figure 1 shows the time-averaged compensated horizontal wavenumber kinetic energy spectra for hyperviscosity simulations, DNS and LES, where the horizontal energy spec5/3 tra are multiplied by kh hǫi−2/3 . The horizontal axes are normalized by the buoyancy wavenumber kb , for which in each panel kh /kb = 1 shows the location of the buoyancy wavenumber. Constant horizontal compensated spectra in the intermediate horizontal wavenumbers imply a −5/3 power-law in the inertial subrange of the horizontal wavenumber kinetic energy spectra that is in line with the Lindborg (2006) hypothesis for stratified turbulence. For DNS, LES and hyperviscosity simulations, the −5/3 spectral slope is seen over the horizontal wavenumbers 0.5 . kh /kb . 1 when N = 0.5, 0.3 . kh /kb . 0.8 when N = 1, 0.2 . kh /kb . 0.9 when N = 2, and 0.1 . kh /kb . 0.8 when N = 4 (see panels a-d in figure 1). Also, these constant compensated spectra are followed by a bump (or shallowness) around kh ∼ kb that indicates a non-local energy transfer from larger scales to the buoyancy scale Lb = 2π/kb (consistent with hyperviscosity simulations, DNS and LES of Waite 2011; Khani & Waite 2013, 2014, 2015). These bumps are seen around kb ≈ 8, 16, 30 and 60, when N = 0.5, 1, 2 and 4 (see the vertical arrows in figure 1, indicating the location of the buoyancy wavenumber kb ). It is clear that the bump (or shallowness) occurs at kh /kb ≈ 1 in all panels of figure 1. For example, for the dynamic Smagorinsky simulation with n = 512 and N = 2, the bump (or shallowness) happens around kh ≈ 30 that is very close to kb = 29.9 (see table 1 and the dash double-dotted red line in figure 1c). Overall, these results show that LES

6

S. Khani

Table 1. List of numerical simulations with LES. Dynamic Smagorinsky N d5N0.5 d5N1 d5N2 d5N4 d5N6 d3N0.5 d2N0.5 d2N1 d2N2 d2N4 d2N6

0.5 1 2 4 6 0.5 0.5 1 2 4 6

Smagorinsky N

n

Lb

Lo



hǫi

hǫp i

512 512 512 512 512 384 256 256 256 256 256

0.79 0.39 0.21 0.11 0.08 0.85 0.82 0.41 0.21 0.12 0.09

0.162 0.060 0.023 0.008 0.005 0.186 0.177 0.066 0.023 0.008 0.005

0.019 0.019 0.019 0.019 0.019 0.025 0.038 0.038 0.038 0.038 0.038

8.3 × 10−5 9.1 × 10−5 1.1 × 10−4 1.1 × 10−4 1.2 × 10−4 1.1 × 10−4 9.9 × 10−5 1.1 × 10−4 1.1 × 10−4 1.1 × 10−4 1.2 × 10−4

3.9 × 10−5 3.6 × 10−5 3.7 × 10−5 3.3 × 10−5 3.1 × 10−5 5.0 × 10−5 4.4 × 10−5 4.1 × 10−5 3.2 × 10−5 2.7 × 10−5 2.7 × 10−5

n

Lb

Lo



hǫi

hǫp i

S7N2 S7N4 S7N6 S5N2 S5N4 S2N0.5 S2N1 S2N2

2 4 6 2 4 0.5 1 2

768 768 768 512 512 256 256 256

0.20 0.10 0.07 0.20 0.11 0.81 0.40 0.20

0.022 0.008 0.005 0.024 0.009 0.178 0.066 0.022

0.012 0.012 0.012 0.019 0.019 0.038 0.038 0.038

9.8 × 10−5 1.1 × 10−4 1.1 × 10−4 1.2 × 10−4 1.3 × 10−4 1.0 × 10−4 1.1 × 10−4 1.0 × 10−4

2.9 × 10−5 2.2 × 10−5 1.8 × 10−5 3.0 × 10−5 2.4 × 10−5 4.0 × 10−5 3.4 × 10−5 2.0 × 10−5

Kraichnan

N

n

Lb

Lo



hǫi

hǫp i

K7N2 K7N6 K2N0.5 K2N1 K2N2 K2N4 K2N6

2 6 0.5 1 2 4 6

768 768 256 256 256 256 256

0.21 0.08 0.82 0.40 0.20 0.11 0.08

0.023 0.005 0.174 0.063 0.023 0.009 0.005

0.012 0.012 0.038 0.038 0.038 0.038 0.038

1.1 × 10−4 1.2 × 10−4 9.6 × 10−5 1.0 × 10−4 1.1 × 10−4 1.2 × 10−4 1.3 × 10−4

4.0 × 10−5 3.6 × 10−5 4.6 × 10−5 4.1 × 10−5 3.7 × 10−5 2.9 × 10−5 2.7 × 10−5

hE(t)i Lh 0.0040 0.0039 0.0043 0.0049 0.0056 0.0046 0.0042 0.0043 0.0044 0.0055 0.0070

0.81 0.95 1.55 2.32 2.80 1.02 1.23 1.51 2.15 2.65 2.87

hE(t)i Lh 0.0039 0.0043 0.0048 0.0042 0.0050 0.0042 0.0040 0.0040

1.61 2.87 3.00 2.17 2.76 1.31 1.82 2.65

hE(t)i Lh 0.0045 0.0053 0.0043 0.0041 0.0042 0.0050 0.0060

0.98 2.03 1.03 1.22 1.72 2.60 2.90

Lv 0.61 0.41 0.40 0.31 0.28 0.93 1.02 0.72 0.55 0.53 0.52 Lv 0.36 0.30 0.32 0.41 0.42 1.06 0.73 0.68 Lv 0.33 0.20 0.94 0.65 0.46 0.39 0.38

F rh F rv 0.98 0.41 0.14 0.05 0.03 0.83 0.67 0.27 0.10 0.05 0.03

1.30 0.95 0.53 0.36 0.29 0.91 0.80 0.57 0.38 0.23 0.17

F rh F rv 0.12 0.03 0.02 0.09 0.04 0.62 0.22 0.08

0.56 0.33 0.22 0.49 0.26 0.76 0.55 0.30

F rh F rv 0.21 0.04 0.80 0.33 0.12 0.04 0.03

0.64 0.40 0.87 0.62 0.44 0.28 0.21

runs capture the fundamental characteristics of stratified turbulence that are in line with DNS and hyperviscosity simulations. The horizontal and vertical SGS energy transfer spectra (eddy dissipation spectra) are shown in figure 2 (see Appendix A, and also Khani & Waite 2014, 2015, for calculations of the SGS eddy dissipation spectra). Interestingly, eddy dissipation spectra peak at large horizontal and small vertical scales that is consistent with anisotropic dissipation in stratified turbulence (in line with effective eddy viscosity results of Khani & Waite 2013; Remmler & Hickel 2014; Khani & Waite 2016). The horizontal wavenumber eddy dissipation spectrum for the Kraichnan LES with n = 768 and N = 2 has a cusp at kh ≈ kc because the Kraichnan eddy viscosity model increases near kc (see e.g. Kraichnan 1976; M´etais & Lesieur 1992; Khani & Waite 2013). Increased stratification or decreased resolution increases the eddy dissipation by squeezing the buoyancy

7

Mixing efficiency in large-eddy simulations of stratified turbulence

Table 2. List of numerical simulations with DNS and hyperviscosity. DNS

N

n

Lb

Lo

D7N0.5 0.5 768 0.81 0.167 D7N0.75 0.75 768 0.53 0.092 1 768 0.39 0.060 D7N1 2 768 0.20 0.022 D7N2 D2N2 2 256 0.22 0.023

kmax /kd

hǫi

hǫp i

hE(t)i Lh

0.84 0.83 0.83 0.81 0.83

8.9 × 10−5 9.1 × 10−5 9.1 × 10−5 9.7 × 10−5 1.1 × 10−4

4.3 × 10−5 4.0 × 10−5 3.8 × 10−5 3.3 × 10−5 2.5 × 10−5

0.0042 0.0040 0.0039 0.0040 0.0047

Hyperviscosity

N

n

Lb

Lo

kmax /kd

hǫi

hǫp i

h7N2 h7N4 h7N6 h2N0.75 h2N1 h2N1.5 h2N2 h2N4 h2N6

2 4 6 0.75 1 1.5 2 4 6

768 768 768 256 256 256 256 256 256

0.21 0.11 0.08 0.55 0.41 0.28 0.21 0.10 0.05

0.016 0.006 0.004 0.063 0.044 0.027 0.019 0.008 0.005

1.38 1.37 1.34 1.43 1.39 1.38 1.36 1.34 1.34

5.0 × 10−5 5.9 × 10−5 8.8 × 10−5 4.3 × 10−5 4.8 × 10−5 6.2 × 10−5 7.6 × 10−5 1.1 × 10−4 1.1 × 10−4

2.3 × 10−5 2.8 × 10−5 3.8 × 10−5 2.5 × 10−5 2.7 × 10−5 3.0 × 10−5 3.0 × 10−5 2.8 × 10−5 2.4 × 10−5

0.80 0.86 0.97 1.52 2.6

Lv

F rh F rv Reb

0.61 0.48 0.41 0.31 0.93

1.01 0.62 0.40 0.13 0.09

hE(t)i Lh

Lv

0.0045 0.0048 0.0060 0.0043 0.0042 0.0044 0.0047 0.0044 0.0027

0.86 1.40 1.93 1.10 1.17 1.50 1.95 3.34 3.64

0.32 0.20 0.16 0.78 0.65 0.46 0.39 0.31 0.31

1.33 1.04 0.95 0.65 0.24

16.6 7.6 4.2 1.1 0.3

F rh F rv 0.24 0.08 0.04 0.50 0.35 0.19 0.11 0.03 0.01

scale thickness (vertical layer thickness) towards the dissipation scale or by increasing the dissipation scale (as seen in hyperviscosity simulations and DNS of Waite & Bartello 2004; Bartello & Tobias 2013). In addition, this anisotropic picture of eddy dissipation spectra is in agreement with the effective eddy dissipation spectra when a test-filter scale ∆test = π/kc is employed around the buoyancy scale Lb in DNS results (see below in figure 8). 4.2. Length scales We can calculate the vertical and horizontal length scales using the corresponding averaged vertical and horizontal energy spectra in DNS, LES and hyperviscosity simulations of stratified turbulence as follows, R∞ R∞ 2π 0 hE(kv )idkv 2π 0 hE(kh )idkh Lv = R ∞ , Lh = R ∞ , (4.1) kv hE(kv )idkv kh hE(kh )idkh 0 0

where hEi, kh and kv are the averaged directional resolved energy spectrum, and the horizontal and vertical wavenumbers, respectively. Figure 3 shows the ratio Lv /Lb versus Lb /∆, which is an indication to show how sufficiently the buoyancy scale Lb is resolved in LES of stratified turbulence. For comparison, the results from hyperviscosity simulations and DNS are also shown. Khani & Waite (2014, 2015) have shown that there is a threshold on the grid spacing ∆, for which LES of stratified turbulence captures the fundamental characteristics of stratified turbulence. This threshold depends on Lb and varies with different SGS models: Lb /∆ > 4.17 for the dynamic Smagorinsky model, Lb /∆ > 5.89 for Smagorinsky model, and Lb /∆ > 2.13 for the Kraichnan eddy viscosity model (see Khani & Waite 2014, 2015, these thresholds are also shown with vertical dashed lines in figure 3). It is clear that the ratio Lv /Lb approaches to 1 when  Lb is well resolved in LES, or when the buoyancy Reynolds number Reb = hǫi/ νN 2 is sufficiently large in

0.66 0.55 0.50 0.70 0.63 0.61 0.54 0.32 0.16

8

S. Khani (a) hE(kh )i × kh /hǫi2/3

10

0

5/3

10

D7N0.5 d5N0.5 d3N0.5 d2N0.5 S2N0.5 K2N0.5

10−2

0.1

1

10−2

D7N1 d5N1 d2N1 S2N1 K2N1 0.1

10

1

kh/kb

10

kh/kb (c)

(d)

kb hE(kh )i × kh /hǫi2/3

kb

10

0

5/3

100

5/3

hE(kh )i × kh /hǫi2/3

(b)

kb

5/3

hE(kh )i × kh /hǫi2/3

kb 0

10−2

D7N2 h7N2 d5N2 d2N2 S7N2 K7N2 0.1

1

10−2

h7N4 d5N4 d2N4 S7N4 S5N4 K2N4

10−4 0.01

10

0.1

kh/kb

1 kh/kb

Figure 1. The averaged compensated horizontal wavenumber kinetic energy spectra for DNS, LES and hyperviscosity simulations: (a) N = 0.5, (b) N = 1, (c) N = 2 and (d) N = 4. The horizontal axises are normalized by the buoyancy wavenumber kb , and the vertical arrows show the location of kb , where the bump (or spectral shallowness) starts. The spectra are averaged over 375 6 t 6 450. 0

0 (b) -2 kv Tr (kv ) ×105

kh Tr (kh ) ×105

-2

-4 d5N4 d5N2 d5N1 S7N2 K7N2 K2N2

-6

-8 1

(a)

10

100 kh

-4

-6 d5N4 d5N2 -8 d5N1 S7N2 K7N2 K2N2 -10 1

10

100 kv

Figure 2. The (a) horizontal and (b) vertical wavenumber SGS eddy dissipation spectra for LES. These spectra are calculated at t = 450, and are multiplied by wavenumber to preserve area on the log-linear axes.

DNS (note that ∆ ∝ Ld in DNS, so the ratio Lb /∆ is a measure of the size of the subbuoyancy-scale inertial range; see figure 3, and also figure 6 in below). This trend is in line with an important feature of stratified turbulence, in which the vertical length scale is scaled by the buoyancy scale Lb (i.e., Lv ∼ Lb , see e.g. Waite 2011; Khani & Waite 2013).

Mixing efficiency in large-eddy simulations of stratified turbulence

9

6 DSM SMG KRN DNS HYP

5

Lv /Lb

4 3 2 1 0 0 10

101

102

Lb /∆

Figure 3. The ratio Lv /Lb versus Lb /∆ for DNS, LES and hyperviscosity simulations in stratified turbulence. The dynamic Smagorinsky, regular Smagorinsky and Kraichnan SGS models are considered in LES runs. Vertical dashed lines (from left to right for the Kraichnan, dynamic Smagorinsky and Smagorinsky parameterizations, respectively) represent the threshold on Lb /∆, above which LES could capture fundamental characteristics of stratified turbulence (see e.g. Khani & Waite 2014, 2015). In DNS runs, ∆ ∝ Ld , so the ratio Lb /∆ is a measure of the size of the sub-buoyancy-scale inertial range.

Therefore, we can simplify the equation for vertical Froude number F rv = urms /(N Lv ) as follows urms Lb F rv ∼ ∼ ∼ 1, (4.2) N Lb Lb when Lb /∆ ≫ 1 in LES of stratified turbulence. On the other hand, when Lb /∆ . 1, Lv ≫ Lb , and so F rv =

urms Lb ∼ ≪ 1, N Lv Lv

(4.3)

in LES of stratified turbulence, as seen in figure 3 (consistent with the behavior of viscosity-affected stratified turbulence when Reb . 1, as seen in DNS run with n = 256 and N = 2; also see Brethouwer et al. 2007). In the next section, we will show that these two different vertical Frounde number regimes in equations (4.2 & 4.3) have critical influence on scale analysis of the SGS mixing efficiency in LES of stratified turbulence. 4.3. Mixing efficiency According to equations (2.4-2.7), we can write the irreversible SGS mixing efficiency in LES as follows  h β 2 /P rt (∂ ρ¯/∂xj ) (∂ ρ¯/∂xj )i . (4.4) γi = h2¯ sij s¯ij i + h(β 2 /P rt ) (∂ ρ¯/∂xj ) (∂ ρ¯/∂xj )i Next, we are interested in performing scale analysis on γi by using the scale of resolved variables in the right-hand side of equation (4.4). If we assume Lv ≪ Lh , which is a reasonable assumption in stratified turbulence (see table 1), we can write (following to Khani & Port´e-Agel 2017) g ∂ ρ¯ ∼ N 2, ρ0 ∂xj

s¯ij ∼

urms , Lv

(4.5)

and so, γi ∼

1 N 2 /P rt ∼ . 2 (u2rms /L2v ) + N 2 /P rt 2P rt F rv2 + 1

(4.6)

10

S. Khani

Moreover, in the stratified turbulence regime F rv ∼ 1 (as is given by equation 4.2, also see e.g. Billant & Chomaz 2001; Waite & Bartello 2004; Lindborg 2006; Waite 2011; Khani & Waite 2013; Khani 2015), which results in γi ∼

1 . 2P rt + 1

(4.7)

This scaling theory is supported by simulations at different turbulent Prandtl numbers P rt , which are discussed below in section 4.6. We can write an equality relation for the scaling theory (4.7) as follows c , (4.8) γi = 2P rt + 1 where c is a proportionality constant, and we may estimate it using LES data. We have employed the least-squares method to measure c using our LES runs, where the turbulent Prandtl number P rt varies from 0.2 to 2. It is shown that the estimated value for the coefficient c depends on the SGS scheme and grid spacing, but is fairly close to 1 (see Appendix B). Therefore, we set c ≈ 1 in the equality relation (4.8). Assuming P rt ≈ 1 (see below in figure 9 and also in Khani & Waite 2016), we can obtain a theoretical estimate for the irreversible SGS mixing efficiency γi in the stratified turbulence regime as follows 1 γi ≈ , (4.9) 3 which is slightly larger than the value that is observed in recent DNS studies (see e.g. Venayagamoorthy & Stretch 2010; Maffioli et al. 2016). It is important to mention that all parameters and variables are scaled except the turbulent Prandtl number, which is approximated by 1, to obtain equation (4.9); therefore, this equation gives an estimate for the scale of the irreversible mixing efficiency γi . In the rest of this section, we will show that the LES results are in good agreement with those from DNS and hyperviscosity simulations and the scale analysis in equation (4.9). Figure 4 shows the irreversible SGS mixing efficiency γi versus the vertical Froude number F rv in panel (a) and the horizontal Froude number F rh in panel (b) for dynamic Smagorinsky, regular Smagorinsky and Kraichnan SGS models. For comparison, the results from hyperviscosity simulations and DNS are also shown. Here, we have considered a 2π factor in definitions of F rv and F rh to remove the effects of 2π factor in definitions of Lv and Lh , respectively, as are given in equation (4.1). Interestingly, the curves for the irreversible SGS mixing efficiency in LES are fairly close to those of the irreversible mixing efficiency in DNS. The mixing efficiency γi is small when F rv is small, and it goes to the value 1/3 as F rv increases towards 1 (figure 4a). The violet cross sign in figure 4(a) shows the theoretical value of γi (as given by equation 4.9). Importantly, equation (4.9) is only valid in LES approach when F rv ∼ 1 because this is the situation that Lv ∼ Lb (see equation 4.2). As F rv → 1, LES results approach to the theoretical value γi ≈ 1/3 (see figure 4a). It is interesting that DNS and the dynamic Smagorinsky LES perform similarly (see black and blue lines in figure 4). The maximum irreversible mixing efficiency γi is given by 0.33 in DNS run with n = 768 and N = 0.5, and it is given by 0.32 (0.31) in the dynamic Smagorinsky simulation with n = 512 (256) and N = 0.5. This behavior shows that the dynamic Smagorinsky model performs reasonably well even at lower resolution simulations in comparison with high resolution DNS runs. The performance of the regular Smagorinsky model is less efficient in comparison with the Kraichnan model and the maximum of γi is 0.29 and 0.32, respectively, in the Smagorinsky and Kraichnan simulations with n = 256 and N = 0.5. As F rv (or F rh )

Mixing efficiency in large-eddy simulations of stratified turbulence

11

0.4 0.35

γi

0.3

DSM SMG KRN DNS HYP

(F rv = 1, γi = 1/3)

0.25 0.2 0.15 0.1 0

(a) 0.2

0.4

0.6

0.8

F rv =

1

2πurms NLv

1.2

1.4

0.4 (b) 0.35

γi

0.3 0.25 0.2

DSM SMG KRN DNS HYP

0.15 0.1 0

0.2

0.4

0.6

F rh =

2πurms N Lh

0.8

1

1.2

Figure 4. The irreversible mixing efficiency γi versus (a) the vertical Froude number F rv and (b) the horizontal Froude number F rh , for DNS, LES and hyperviscosity simulations. We have considered 2π factors in definitions of F rv and F rh to remove the effects of 2π factors, which are used in definitions of Lv and Lh , respectively, as seen in equation (4.1). The violet cross sign in panel (a) shows the theoretical scale analysis of γi for LES approach as given in equations (4.9).

decreases, the irreversible SGS mixing efficiency γi reduces (figure 4). Also, we should mention that simulations are mostly in the range of small horizontal Froude number, and F rh increases as N decreases (moving from the strongly stratified to weakly stratified regime, see section 4.4). The hyperviscosity simulations show slightly higher γi even at F rv ≈ 0.7 (figure 4a). This trend could be due to the effects of hyperviscosity that constrains the dissipation range to much smaller scales in comparison with regular viscosity, and so slightly decreases the averaged kinetic dissipation rates (see table 2). We have suggested different regimes for F rv in LES of stratified turbulence based on the resolution of the buoyancy scale Lb (see equations 4.2 & 4.3). Hence, it will be important to investigate if γi does depend on the resolution of Lb or not? We will do a similar investigation to study the dependence of γi on the resolution of Lo in the next section. Figure 5 shows γi versus Lb /∆ for different SGS parameterizations in LES, and also for hyperviscosity simulations and DNS (vertical dashed lines in figure 5 represent the thresholds on Lb /∆, respectively from left to right, for the Kraichnan, dynamic Smagorinsky and Smagorinsky models). It is clear that for all SGS parameterizations, as the ratio Lb /∆ increases, the irreversible SGS mixing efficiency γi rises and approaches to the theoretical value 1/3 (see equation 4.9 and figure 5). This trend suggests that if the buoyancy scale Lb is well resolved in LES, the mixing efficiency γi is in line with the results of the theoretical scale analysis in LES of stratified turbulence, and if the threshold Lb /∆ is violated, LES of stratified turbulence underestimates the irreversible SGS mixing efficiency γi (figure 5). This behavior is in agreement with DNS runs, in which the irreversible mixing efficiency γi increases monotonically with Lb /∆ (see black line in figure 5; note that ∆ ∝ Ld in DNS). Figure 6 shows the irreversible mixing efficiency

12

S. Khani 0.4 0.35

γi

0.3 0.25 DSM SMG KRN DNS HYP

0.2 0.15 0.1 0 10

101

102

Lb /∆

Figure 5. The irreversible mixing efficiency γi versus the ratio Lb /∆ for DNS, LES and hyperviscosity simulations. See the caption in figure 3 for more information about the vertical dashed lines.

0.34 0.32 0.3

γi

0.28 0.26 0.24 0.22 DNS

0.2 0.18 −1 10

100

Reb

101

102

Figure 6. The irreversible mixing efficiency γi versus the buoyancy Reynolds number Reb = hǫi/ νN 2 for DNS runs.

versus the buoyancy Reynolds number Reb in DNS runs. Clearly, γi increases respectively from 0.18 to 0.33 when the buoyancy Reynolds number increases from Reb ≈ 0.3, where the layered stratified flows are mainly affected by viscous sublayers, to Reb ≈ 17, where overturning and Kelvin-Helmholtz instabilities are ubiquitous in stratified turbulence (in line with DNS results of Shih et al. 2005). Overall, depending on the resolution of the buoyancy scale Lb , LES of stratified turbulence will underestimate γi if Lb in not adequately resolved (similar to DNS runs with Reb . 1), or LES will correctly resemble γi if Lb is well resolved, i.e. Lb /∆ & 10 as seen in figures 3 and 5 (similar to DNS runs with Reb & 10). In addition, the regular Smagorinsky parameterization still underestimates γi when N is large even at the highest resolution simulations, but both the Kraichnan and dynamic Smagorinsky models represents reasonable results when Lb is more likely well resolved (see figures 4 and 5). This behavior is due to the highly dissipative nature of the regular Smagorinsky model (as is also reported in Khani & Waite 2014, 2015). These results suggest that reproducing correct mixing efficiency in LES approach requires much better resolution of Lb in comparison with, for example, reproducing correct energy spectra in LES. In the next section, we study the effects of resolving Lo in mixing efficiency for LES of stratified turbulence.

Mixing efficiency in large-eddy simulations of stratified turbulence

13

4.4. Ozmidov scale effects In stratified turbulence, the ratio of the Ozmidov scale Lo to the buoyancy scale Lb can be scaled as Lo 1/2 ∼ F rh , (4.10) Lb where we have used Taylor’s hypothesis Lh ∼ u3rms /ǫ. When F rh ≪ 1, which is the case in strongly stratified turbulent flows, Lo ≪ Lb (see equation 4.10). Therefore, if Lo is resolved (i.e. Lo & ∆), the resolution of Lb is ensured (i.e. Lb ≫ ∆). In this situation Lv ∼ urms /N , and we can scale F rh as follows F rh =

Lv urms ∼ , N Lh Lh

(4.11)

where Lv ≪ Lh (anisotropic motions). Using equations (4.10 & 4.11), we can conclude that  1/2 Lo Lv ∼ , (4.12) Lb Lh in the strongly stratified turbulence regime. If the turbulence regime moves from strongly stratified to weakly stratified (and so forth to the isotropic Kolmogorov turbulence), Lv → Lh (isotropic motions), and thus Lo → Lb (see equation 4.12). As a result, weakly stratified means F rh ∼ O(1), which is in agreement with equation (4.10). Figure 7 shows the ratio Lv /Lb and γi versus Lo /∆ in panels (a) and (b), respectively, for DNS, LES and hyperviscosity simulations. It is clear that Lv ≈ Lb when Lo /∆ & 1 (figure 7a). This trend is expected because if Lo is more likely resolved in LES of stratified turbulence where F rh ≪ 1, it is guaranteed that Lb is well-resolved (see equation 4.10). In addition, the irreversible mixing efficiency γi moves towards the theoretical value 1/3 if Lo is resolved in LES of stratified turbulence (figure 7b). These results show that if ∆ . Lo in LES of stratified turbulence, mixing efficiency is correctly reproduced. Overall, these results are in line with those from figures 3 and 5. Consequently, if ∆ . Lo , irreversible mixing efficiency goes to the theoretical value 1/3 (see figures 5 and 7). In addition, the fundamental characteristics of stratified turbulence, such as anisotropic features, are more likely captured because the buoyancy scale Lb is well resolved (see e.g. Khani & Waite 2014, 2015; Khani & Port´e-Agel 2017) when ∆ . Lo . 4.5. An a priori study using DNS results In this section, we calculate the spectra of effective SGS energy transfer (eddy dissipation) and turbulent Prandtl number by applying a test cutoff wavenumber kc , which is around the buoyancy wavenumber kb , to DNS data (see Appendix A, and also Khani & Waite 2016, for calculations of the effective eddy dissipation spectra). Figure 8 shows the horizontal and vertical wavenumber spectra for effective eddy dissipation. These plots should be compared with the SGS eddy dissipation spectra that are shown in figure 2. We consider kc = 20 for DNS cases with N = 0.5, 0.75 and 1, and kc = 40 for the DNS case with N = 2 (these values are around the buoyancy wavenumber kb = 2π/Lb ; see table 2). The effective eddy dissipation spectra also peak at large horizontal and small vertical scales (in line with the anisotropic dissipation picture, and consistent with the corresponding SGS eddy dissipation that are shown in figure 2). These trends strongly suggest that the dissipation mechanism in stratified turbulence includes anisotropic features around the buoyancy scale Lb , and that is different than the dissipation mechanism in isotropic turbulence. When the test-filter wavenumber kc is fixed, increased buoyancy

14

S. Khani

Figure 7. (a) The ratio Lv /Lb and (b) irreversible mixing efficiency γi versus the ratio Lo /∆ for DNS, LES and hyperviscosity simulations. In DNS, ∆ ∝ Ld , and therefore Lo /∆ is a measure of Reb (i.e. large Lo /∆ in DNS refers to stratified turbulence regime with Reb ≫ 1). 0

0

-2 kv Tr (kv ) ×105

kh Tr (kh ) ×105

-2

-4

-6

D7N2 D7N1 D7N0.75 D7N0.5

(a) -8 0.1

k = kh /kc

-6 D7N2 D7N1 D7N0.75 D7N0.5

-8 (b) -10 1

+

-4

0.1

1 k+ = kv /kc

Figure 8. The (a) horizontal and (b) vertical wavenumber effective eddy dissipation spectra in DNS runs. The spectra are calculated at t = 450, and are multiplied by wavenumber to preserve area on the log-linear axes. We employ the test cutoff wavenumber kc = 40 in the DNS run with N = 2, and kc = 20 in other DNS runs. These plots should be compared to the corresponding SGS eddy dissipation spectra that are shown in figure 2.

frequency N results in increasing the nonlocal horizontal (i.e. kh ≪ kc ) and local vertical (i.e. kv ∼ kc ) effective eddy dissipation in simulations with N = 0.5, 0.75 and 1, because the buoyancy scale Lb becomes smaller than the filter width ∆test = π/kc , and thus the anisotropic motions are mainly filtered (figure 8 and table 2). Overall, resolving the buoyancy scale in LES approach is a critical threshold in order to capture the fundamental characteristics of stratified turbulence (in line with LES and DNS results of Khani & Waite 2014, 2015, 2016; Khani & Port´e-Agel 2017). The spectra of effective turbulent Prandtl number are shown in figure 9. We define

Mixing efficiency in large-eddy simulations of stratified turbulence

15

2 1.5

P rt (k|kc )

1 0.5 0 D7N2 D7N1 D7N0.75 D7N0.5

-0.5 -1 0.1

1 k/kc

Figure 9. The effective turbulent Prandtl number spectra for DNS runs. The spectra are calculated at t = 450. We employ the test cutoff wavenumber kc = 40 in the DNS run with N = 2, and kc = 20 in other DNS runs. The solid and dashed magenta lines indicate values of zero and 1, respectively.

P rt (k|kc ) as the ratio of the effective eddy viscosity νe (k|kc ) over the effective eddy diffusivity De (k|kc ) (as are given in Appendix A, and also in Khani & Waite 2016). The turbulent Prandtl number spectra show approximately constant trends when k/kc & 0.2 (consistent with the results of Khani & Waite 2016, and this also demonstrates that the assumption of P rt ≈ 1 in this current work is a reasonable approximation). 4.6. Turbulent Prandtl number effects In stratified turbulence, P rt ≈ 1 (as is shown in figure 9; and also see Salehipour & Peltier 2015; Khani & Waite 2016), which leads to a scale of around 1/3 for the irreversible mixing efficiency in LES of stratified turbulence (see equation 4.9) that is in line with DNS results. Nevertheless, it is still very interesting to test the scaling theory (4.7) with different turbulent Prandtl numbers. We have performed extra LES runs when P rt varies from 0.2 to 2. Tables 3, 4 and 5 show simulation parameters and averaged variables for these LES runs. Using these new simulations, we can evaluate the dependency of the irreversible mixing efficiency γi to the turbulent Prandtl number P rt as is given by scaling theory (4.7), and also we can estimate the proportionality coefficient c in equation (4.8), as shown in Appendix B. Figure 10 shows γi versus P rt for the dynamic Smagorinsky, regular Smagorinsky and Kraichnan SGS models. For comparison, the scaling theory (4.7) is also shown with black dashed line. It is clear that γi increases (decreases) when P rt decreases (increases), and that LES results are in line with the scaling theory (figure 10). When P rt & 1, the scaling theory (4.7) agrees very well with LES results, however, it slightly overestimates γi when P rt < 1. For example, the dynamic Smagorinsky model with n = 384, N = 1 and P rt = 2 gives γi = 0.20, which equals to the scaling theory value in equation (4.7) when P rt = 2. Similar results are given by the Smagorinsky model with n = 512 and N = 1, and the Kraichnan model with n = 384 and N = 1. When P rt = 0.5, the dynamic Smagorinsky model with n = 384 and N = 0.5 gives γi = 0.38,

16

S. Khani

Table 3. List of numerical simulations with LES when P rt = 0.2. Dynamic Smagorinsky N d5N0.5 d3N0.5 d3N1

Lb

Lo



hǫi

hǫp i

hE(t)i Lh

Lv

F rh F rv

0.5 512 0.79 0.149 0.019 7.0 × 10−5 5.5 × 10−5 0.0040 0.86 0.63 0.92 1.25 0.5 384 0.85 0.165 0.025 8.6 × 10−5 7.5 × 10−5 0.0046 1.11 0.97 0.77 0.88 1 384 0.41 0.059 0.025 8.7 × 10−5 9.5 × 10−5 0.0042 1.35 0.63 0.30 0.65

Smagorinsky N S5N1 S2N0.5

n

n

Lb

Lo



hǫi

hǫp i

hE(t)i Lh

Lv

F rh F rv

1 512 0.40 0.057 0.019 8.4 × 10−5 6.4 × 10−5 0.0041 1.31 0.60 0.30 0.67 0.5 256 0.81 0.160 0.038 8.1 × 10−5 6.7 × 10−5 0.0042 1.44 1.10 0.56 0.73

Kraichnan

N

K5N2 K3N1

2 1

n

Lb

Lo



hǫi

hǫp i

hE(t)i Lh

Lv

F rh F rv

512 0.20 0.019 0.019 7.7 × 10−5 6.1 × 10−5 0.0041 1.39 0.35 0.14 0.57 384 0.41 0.058 0.025 8.4 × 10−5 7.0 × 10−5 0.0042 1.05 0.59 0.39 0.69

Table 4. List of numerical simulations with LES when P rt = 0.5. Dynamic Smagorinsky N d5N1 d3N0.5 d3N1

Kraichnan K5N1 K2N0.5 K2N1

Lb

Lo



hǫi

hǫp i

hE(t)i Lh

Lv

F rh F rv

1 512 0.39 0.057 0.019 8.1 × 10−5 4.5 × 10−5 0.0039 1.00 0.42 0.39 0.93 0.5 384 0.85 0.176 0.025 9.8 × 10−5 6.0 × 10−5 0.0046 1.05 0.95 0.81 0.89 1 384 0.41 0.063 0.025 1.0 × 10−4 5.7 × 10−5 0.0043 1.26 0.64 0.33 0.64

Smagorinsky N S5N0.5 S2N0.5

n

n

Lb

Lo



hǫi

hǫp i

hE(t)i Lh

Lv

F rh F rv

0.5 512 0.83 0.173 0.019 9.5 × 10−5 5.4 × 10−5 0.0044 0.96 0.88 0.86 0.94 0.5 256 0.81 0.168 0.038 9.0 × 10−5 5.1 × 10−5 0.0042 1.36 1.10 0.60 0.74 N

n

Lb

Lo



hǫi

hǫp i

hE(t)i Lh

Lv

F rh F rv

1 512 0.40 0.061 0.019 9.6 × 10−5 5.3 × 10−5 0.0041 0.87 0.53 0.46 0.75 0.5 256 0.82 0.168 0.038 9.0 × 10−5 5.6 × 10−5 0.0043 1.05 0.94 0.78 0.87 1 256 0.40 0.060 0.038 9.2 × 10−5 5.4 × 10−5 0.0041 1.26 0.65 0.32 0.61

which is below the given value by the theoretical scale analysis curve with γi ∼ 0.5. Similar trends are seen in the Smagorinsky and Kraichnan models. Also when P rt = 0.2, the dynamic Smagorinsky model with n = 384 and N = 1 gives γi = 0.52, which is below the given estimate with scaling theory γi ∼ 0.71. It is worthwhile mentioning that in stratified turbulence regime in which P rt ≈ 1 (see e.g. Khani & Waite 2016), LES results agrees very well with the theoretical scale analysis value with γi ∼ 0.33.

17

Mixing efficiency in large-eddy simulations of stratified turbulence

Table 5. List of numerical simulations with LES when P rt = 2. Dynamic Smagorinsky N d5N0.5 d3N0.5 d3N1 d3N2

0.5 0.5 1 2

n

Lb

Lo



hǫi

hǫp i

512 384 384 384

0.79 0.85 0.41 0.21

0.170 0.195 0.069 0.025

0.019 0.025 0.025 0.025

9.2 × 10−5 1.2 × 10−4 1.2 × 10−4 1.3 × 10−4

3.2 × 10−5 3.9 × 10−5 3.1 × 10−5 2.7 × 10−5

n

Lb

Lo



hǫi

hǫp i

Smagorinsky N S5N1 S2N0.5 Kraichnan K3N1 K2N0.5

hE(t)i Lh 0.0040 0.0046 0.0043 0.0045

0.78 1.00 1.15 1.66

hE(t)i Lh

Lv 0.60 0.91 0.63 0.46 Lv

F rh F rv 1.01 0.85 0.36 0.13

1.32 0.93 0.65 0.46

F rh F rv

1 512 0.40 0.069 0.019 1.2 × 10−4 2.8 × 10−5 0.0041 1.11 0.59 0.36 0.68 0.5 256 0.82 0.186 0.038 1.1 × 10−4 3.1 × 10−5 0.0043 1.29 1.05 0.64 0.78 N

n

Lb

Lo



hǫi

hǫp i

hE(t)i Lh

Lv

F rh F rv

1 384 0.41 0.069 0.025 1.2 × 10−4 3.3 × 10−5 0.0042 0.95 0.58 0.43 0.71 0.5 256 0.82 0.186 0.038 1.1 × 10−4 3.8 × 10−5 0.0043 1.02 0.93 0.80 0.88

DSM

0.7

SMG

γi

0.6

KRN

0.5

1 2P rt +1

0.4 0.3 0.2 0.1

0.5

1

P rt

1.5

2

Figure 10. The irreversible mixing efficiency γi versus the turbulent Prandtl number P rt in LES. The black dashed line shows the scaling theory (4.7).

Overall, these results suggest that the scaling theory (4.7) agrees very well with LES results when P rt & 1. When P rt < 1, LES results and theoretical scale analysis still show similar behavior (i.e. decreased P rt leads to an increase in γi in LES runs and scaling theory). Nevertheless, γi is overestimated when P rt < 1, and this trend might be due to the dynamics of turbulent flows with P rt ≪ 1 compared to those with P rt & 1, including stratified turbulence, and that we might require to use different SGS parameterizations for eddy diffusivity rather than those are used in equation (2.7) when P rt ≪ 1. This topic is an interesting research plan for future work, and is out of the scope of this current paper about stratified turbulence.

5. Conclusions Mixing efficiency is studied in LES of stratified turbulence using the Smagorinsky, dynamic Smagorinsky and Kraichnan SGS models, and the results are compared with

18

S. Khani

those from DNS and hyperviscosity simulations. If Lv ∼ Lb , we can scale the irreversible SGS mixing efficiency as a function of the vertical Froude number F rv and the turbulent Prandtl number P rt (see equation 4.6). In the stratified turbulence regime where F rv ∼ 1, and assuming P rt ≈ 1 and a unit scaling coefficient, we can obtain a theoretical value 1/3 for γi . This value in LES of stratified turbulence is in line with DNS and hyperviscosity results (and is slightly higher than those are reported in DNS studies of Shih et al. 2005; Venayagamoorthy & Stretch 2010; Maffioli et al. 2016). In addition, our results show that Lv ∼ Lb when ∆ . Lo in LES of stratified turbulence (figures 3 and 7). Clearly, γi → 1/3 when Lo is sufficiently resolved in all three SGS parameterizations in LES, or when Reb ≫ 1 in DNS (figures 5, 6 and 7). In other words, as F rv approaches to 1, γi moves towards the theoretical value 1/3. The irreversible SGS mixing efficiency γi in lower-resolution LES runs with ∆ . Lo , are very close to the irreversible mixing efficiency in high resolution DNS runs (figures 4 and 5). The SGS energy transfer spectra show anisotropic features around the buoyancy scale Lb : the SGS eddy dissipation spectra peak at large horizontal and small vertical scales (see figures 2a,b). This trend is also seen in effective eddy dissipation spectra through an a priori study, in which we employ a test filter scale ∆test ∼ Lb to DNS data (see figures 8a,b). In addition, by increasing (decreasing) the turbulent Prandtl number P rt , the irreversible mixing efficiency γi decreases (increases), in agreement with the theoretical scaling (4.7), see also figure 10. Overall, this work suggests that LES approach can obtain mixing efficiencies that are in agreement with high resolution DNS data if ∆ . Lo in LES. This finding is indeed highly promising in numerical approach of stratified turbulence because LES requires much smaller computational resources in comparison with DNS. For example, the irreversible mixing efficiency γi in the DNS run with n = 768 and N = 0.5 is almost the same as the irreversible SGS mixing efficiency γi in LES approach with n = 256 and N = 0.5, which is 34 = 81 times faster than the DNS run (∆ increases by a factor of 3 in x, y and z directions, and ∆t increases by the same factor; see figures 4, 5, and tables 1, 2). We can write Reb ∼ ReF rh2 using Taylor’s hypothesis (see e.g. Shih et al. 2005; Brethouwer et al. 2007; Khani & Waite 2013), and so for a DNS run with Reb ≫ 1, increased buoyancy frequency by a factor of 2 requires increasing the Reynolds number Re by a factor of 4 4/3 to keep Reb to be still ≫ 1. Since Re ∝ kd , the Kolmogorov wavenumber kd should increase 43/4 ≈ 2.8 times if Re becomes fourfold. As a result, if we double N in a DNS run with n = 768, we require to increase n to be around 2100 to have the same Reb , which will be an expensive simulation. However, LES is an alternative approach that can yield similar results in mixing efficiency through much cheaper simulations because, unlike DNS, LES do not require to resolve the Kolmogorov scale Ld and it just require to adequately resolve the Ozmidov scale Lo (see figure 7). For future work, we are interested in performing LES of stratified turbulence in more practical cases with non-periodic boundary conditions, and with larger turbulent Prandtl (Schmidt) number, which might be relevant to salinity stratification that is common in realistic ocean models (see e.g. Elliott & Venayagamoorthy 2011).

6. Acknowledgement This work benefited from comments by Michael Waite and three anonymous reviewers. Computing resources from Adroit cluster of Research Computing at Princeton University are gratefully appreciated. Financial support from the Atmospheric and Oceanic Sciences (AOS) Program in Princeton University is gratefully acknowledged.

Mixing efficiency in large-eddy simulations of stratified turbulence

19

Appendix A. Eddy dissipation spectrum Let us assume that m(k, ˆ t) is the Fourier transform of the eddy dissipation term −∇ · τ (x, t) in equation (2.1), written as X −∇ · τ (x, t) = m(k, ˆ t)eik·x . (A 1) k Note that τ in LES is parameterized using SGS models, but in DNS, τ can be directly measure by employing a test filter function such as the cutoff or Gaussian filter to DNS data (see e.g. Khani 2015; Khani & Waite 2016; Khani & Port´e-Agel 2017). If we multiply m(k, ˆ t) with the complex conjugate of the Fourier transform of velocity field u ˆ(k, t), and add the result to the complex conjugate of m(k, ˆ t) with u ˆ(k, t), we can compute the kinetic energy transfer spectrum in the following form Tr (k, t) = m ˆ i (k, t)ˆ u∗i (k, t) + m ˆ ∗i (k, t)ˆ ui (k, t),

(A 2)

where Tr (k, t) is the SGS (effective) kinetic eddy dissipation spectrum in the LES (DNS) approach, and ∗ sign denotes the complex conjugate. The eddy dissipation spectrum Tr (k, t) in equation (A 2) can be binned over constant kh and kv to calculate the onedimension horizontal and vertical SGS (effective) eddy dissipation spectrum Tr (kh , t) and Tr (kv , t), respectively, in LES (DNS), as follows X Tr (kh , t) = Tr (k, t), (A 3) |k|=kh X Tr (k, t). (A 4) Tr (kv , t) = |k |=kv In addition, we can define the effective eddy- viscosity and diffusivity coefficients in filtered DNS data as follows Tr (k, t) , 2k 2 E(k, t) Tp (k, t) , De (k|kc , t) = − 2 2k P(k, t) νe (k|kc , t) = −

(A 5) (A 6)

p where k = kh2 + kv2 ; Tp and P are the effective potential- eddy dissipation and energy spectra, respectively (see e.g. Khani & Waite 2016).

Appendix B. A regression analysis for the coefficient c Using the least-squares method, we can estimate the proportionality coefficient c in equation (4.8) as follows PN n n n=1 γi Λ , (B 1) c = PN n n n=1 Λ Λ

where Λ = 1/ (2P rt + 1), and N is the number of LES runs we use to fit a value for c using the least-squares method. We have used those LES runs that meet the grid spacing requirement for the irreversible mixing efficiency γi (i.e. ∆ . Lo ) in this regression analysis, with different P rt ranging from 0.2 to 2. If we use the results of the dynamic Smagorinsky model in our analysis, the fitted value for c is approximately 0.75, and if we use the results of all SGS models,

20

S. Khani

c ≈ 0.72. If the analysis is restricted to P rt > 1, then a larger c ≈ 0.87 is found. In all cases, c is found to be close to 1.

REFERENCES Armenio, V. & Sarkar, S. 2002 An investigation of stably stratified turbulent channel flow using large-eddy simulation. J. Fluid Mech. 459, 1-42. Bartello, P. & Tobias, S. M. 2013 Sensitivity of stratified turbulence to the buoyancy Reynolds number. J. Fluid Mech. 725, 1-22. Billant P. & Chomaz J.-M. 2001 Self-similarity of strongly stratified inviscid flows. Phys. Fluids 13(6), 1645-1651. Brethouwer G., Billant P., Lindborg E. & Chomaz J.-M. 2007 Scaling analysis and simulation of strongly stratified turbulent flows. J. Fluid Mech. 585, 343-368. Caulfield, C. P. & Peltier, W. R. 2000 The anatomy of the mixing transition in homogeneous and stratified free shear layers. J. Fluid Mech. 413, 1-47. Elliott Z. A. & Venayagamoorthy S. K. 2011 Evaluation of turbulent Prandtl (Schmidt) number parameterizations for stably stratified environmental flows. Dyn. Atmos. Oceans 51, 137150. Germano M., Piomelli U., Moin P. & Cabot W. H. 1991 A dynamic subgrid-scale eddy viscosity model. Phys. Fluids A 3(7), 1760-1765. Kaltenback H.-J., Gerz T. & Schumann U. 1994 Large-eddy simulation of homogeneous turbulence and diffusion in stably stratified shear flow. J. Fluid Mech. 280, 1-40. Khani S. 2015 Large eddy simulations and subgrid scale motions in stratified turbulence. UWSpace. http://hdl.handle.net/10012/9209. Khani S. & Port´e-Agel F. 2017 Evaluation of non-eddy viscosity subgrid-scale models in stratified turbulence using direct numerical simulations. Eur. J. Mech. B/Fluids 65, 168-178. Khani S. & Waite M. L. 2013 Effective eddy viscosity in stratified turbulence. J. of Turbl. 14(7), 49-70. Khani S. & Waite M. L. 2014 Buoyancy scale effects in large-eddy simulations of stratified turbulence. J. Fluid Mech. 754, 75-97. Khani S. & Waite M. L. 2015 Large eddy simulations of stratified turbulence: the dynamic Smagorinsky model. J. Fluid Mech. 773, 327-344. Khani S. & Waite M. L. 2016 Backscatter in stratified turbulence. Eur. J. Mech. B/Fluids 60, 1-12. Kraichnan R. H. 1976 Eddy viscosity in two and three dimensions. J. Atmos. Sci. 33, 1521-1536. Lilly D. K. 1967 The representation of small-scale turbulence in numerical simulation experiments. In In NCAR Manuscript 281, National Center for Atmospheric Research, Boulder, CO, USA, 99?164. Lilly D. K., Waco D. E. & Adelfang S. I. 1974 Stratospheric mixing estimated from high-altitude turbulence measurements. J. Appl. Meteor. 13, 488-493. Lindborg E. 2006 The energy cascade in strongly stratified fluid. J. Fluid Mech. 550, 207-242. Lindborg E. & Cho J. Y. N. 2001 Horizontal velocity structure functions in the upper troposphere and lower stratosphere. 2. Theoretical considerations. J. Geophys. Res. 106(10), 1023310241. Maffioli A., Brethouwer G. & Lindborg E. 2016 Mixing efficiency in stratified turbulence. J. Fluid Mech. 794, R3, doi:10.1017/jfm.2016.206. Mashayek A., Salehipour H., Bouffard D., Caulfield C. P., Ferrari R., Nikurashin M., Pletier W. R. & Smyth W. D. 2017 Efficiency of turbulent mixing in the abyssal ocean circulation. Geophys. Res. Lett., doi:10.1002/2016GL072452. M´etais O. & Lesieur M. S. 1992 Spectral large-eddy simulation of isotropic and stably stratified turbulence. J. Fluid Mech. 239, 3157-194. Moin P. & Mahesh K. 1998 Direct numerical simulation: a tool in turbulence research. Annu. Rev. Fluid Mech. 30, 539-578. Nastrom G. D. & Gage K. S. 1985 A climatology of atmospheric wavenumber spectra observed by commercial aircraft. J. Atmos. Sci. 42, 950-960.

Mixing efficiency in large-eddy simulations of stratified turbulence

21

Orszag S. A. 1971 On the elimination of aliasing in finite-difference schemes by filtering highwavenumber components. J. Atmos. Sci. 28, 1074-1074. Osborn T. R. 1980 Estimates of the local rate of vertical diffusion from dissipation measurements. J. Phys. Oceanogr. 10, 83-89. Pardyjak E. R., Monti P. & Fernando H. J. S. 2002 Flux Richardson number measurements in stable atmospheric shear flows. J. Fluid Mech. 459, 307-316. Peltier W. R. & Caulfield C. P. 2003 Mixing efficiency is stratified shear flows. Annu. Rev. Fluid Mech. 35, 135-167. Pope S. B. 2000 Turbulent Flows. Cambridge University Press, Cambridge. Port´e-Agel F., Meneveau C. & Parlange, M. B. 2000 A scale-dependent dynamic model for large-eddy simulation: application to a neutral atmospheric boundary layer. J. Fluid Mech. 415, 261-284. Remmler S. & Hickel S. 2012 Direct and large eddy simulation of stratified turbulence. Int. J. Heat Fluid Flow 35, 13-24. Remmler S. & Hickel S. 2014 Spectral eddy viscosity of stratified turbulence. J. Fluid Mech. 755, R6. http://dx.doi.org/10.1017/jfm.2014.423. Riley J. J. & Lelong M.-P. 2000 Fluid motions in the presence of strong stable stratification. Annu. Rev. Fluid Mech. 32, 613-657. Riley J. J. & Lindborg E. 2008 Stratified turbulence: a possible interpretation of some geophysical turbulence measurements. J. Atmos. Sci. 65, 2416-2424. Salehipour, H. & Peltier, W. R. 2015 Diapycnal diffusivity, turbulent Prandtl number and mixing efficiency in Boussinesq stratified turbulence. J. Fluid Mech. 775, 464-500. Schaefer-Rolffs U. 2016 A generalized formulation of the dynamic Smagorinsky model. Meteorol. Zeitschrift, PrePub, doi:10.1127/metz/2016/0801. Shih L. H., Koseff J. R., Ivey G. N. & Ferziger J. H. 2005 Parameterizations of turbulent fluxes and scales using homogeneous sheared stably stratified turbulence simulations. J. Fluid Mech. 525, 193-214. Smagorinsky J. 1963 General circulation experiments with the primitive equations. I. The basic experiment. Mon. Weather Rev. 91(3), 99-164. Strang E. J. & Fernando H. J. S. 2001 Entrainment and mixing in stratified shear flows. J. Fluid Mech. 428, 349-386. Venayagamoorthy S. K. & Stretch D. D. 2010 On the turbulent Prandtl number in homogeneous stably stratified turbulence. J. Fluid Mech. 644, 359-369. Waite M. L. 2014 Direct numerical simulations of laboratory-scale stratified turbulence, in: T. von Larcher, P. Williams (Eds.), Modeling Atmospheric and Oceanic Flows: Insights from Laboratory Experiments. American Geophys. Union, Washington, DC, 159-175. Waite M. L. 2011 Stratified turbulence at the buoyancy scale. Phys. Fluids 23, 066602. Waite M. L. & Bartello P. 2004 Stratified turbulence dominated by vortical motion. J. Fluid Mech. 517, 281-303. Weinstock J. 1978 Vertical turbulent diffusion in a stably stratified fluid. J. Atmos. Sci. 35, 1022-1027.