Large Eddy Simulations in Astrophysics

In this review, the methodology of large eddy simulations (LES) is introduced and applications in astrophysics are discussed. As theoretical framework, the scale decomposition of the dynamical equations for neutral fluids by means of spatial filtering is explained. For cosmological applications, the filtered equations in comoving coordinates are also presented. To obtain a closed set of equations that can be evolved in LES, several subgrid scale models for the interactions between numerically resolved and unresolved scales are discussed, in particular the subgrid scale turbulence energy equation model. It is then shown how model coefficients can be calculated, either by dynamical procedures or, a priori, from high-resolution data. For astrophysical applications, adaptive mesh refinement is often indispensable. It is shown that the subgrid scale turbulence energy model allows for a particularly elegant and physically well motivated way of preserving momentum and energy conservation in AMR simulations. Moreover, the notion of shear-improved models for inhomogeneous and non-stationary turbulence is introduced. Finally, applications of LES to turbulent combustion in thermonuclear supernovae, star formation and feedback in galaxies, and cosmological structure formation are reviewed.


Introduction
Turbulent flows with high Reynolds numbers are often encountered in computational astrophysics. Examples are the solar wind, stellar convection zones, star-forming clouds, and probably the gas in galaxy clusters. This review concentrates on computational methods that treat turbulence in the limit of high Reynolds numbers by explicitly solving the compressible Euler equations for the large-scale dynamics of the flow, while incorporating small-scale effects such as viscous dissipation into a subgrid-scale model. Since the non-linear turbulent interactions between different scales are at least partially resolved, this type of simulation is called large eddy simulation (LES).
The relative importance of non-linear interactions and viscous damping is specified by the Reynolds number. It is determined by the characteristic velocity V of the flow, its integral length scale L, and the microscopic viscosity ν: The flow becomes turbulent if the non-linear interactions are much stronger than viscous damping. Generally, this happens if Re reaches values greater than a few 10 3 , but Re can become much greater than that. For instance, an estimate for the turbulent convection zone of the Sun is Re ∼ 10 14 [19].
In principle, we can also define a scale-dependent Reynolds number Re( ) = v ( ) /ν, where v ( ) is the typical magnitude of velocity fluctuations on the length scale . The length sale of strong viscous damping is then given by Re( K ) ∼ 1. For incompressible turbulence, substitution of the Kolmogorov-Obukhov scaling law v ( ) ∼ ( ) 1/3 yields [36] 1/3 4/3 K ν ∼ 1 .
Since the mean dissipation rate ∼ V 3 /L, it follows that The problem of high Re is thus a problem of largely different length scales or, equivalently, a high number of degrees of freedom. In a numerical simulation of turbulence, the range of length scales is limited by the grid scale ∆, which is simply the linear size of the grid cells. Only if ∆ K , turbulence can be fully resolved by a so-called direct numerical simulation (DNS). However, DNS become infeasible for very large Re because the total amount of floating point operations (FLOPs) increases with (L/∆) 4 (L/ K ) 4 ∼ Re 3 . The scaling may differ for highly compressible turbulence, but the basic problem remains the same. For a DNS of solar convection over one dynamical time scale, it would be necessary to perform very roughly 10 42 FLOP, which would take far longer than the current age of the Universe on the fastest existing computer.
In practice, however, it is neither feasible nor useful to account for all degrees of freedom in a simulation of high-Re turbulence. To reproduce statistical properties, a much coarser sampling of the degrees of freedom can be quite sufficient. This is why LES encompass only the energy-containing scales and structures dominated by non-linear interactions, which are part of the turbulent cascade down to a cutoff scale much greater than the microscopic dissipation scale. The cutoff scale is given by grid scale ∆. The defining criterion for LES is thus L ∆ K or, equivalently, Re Re(∆) 1 .
Here, Re(∆) ∼ v (∆)∆/ν is the Reynolds number of subgrid-scale turbulence. The product v (∆)∆ can be interpreted as turbulent viscosity of the numerically unresolved eddies of size ∆. The effective Reynolds number of the numerically computed flow is therefore given by This means that LES reduces the number of degrees of freedom by replacing the microscopic viscosity ν by a turbulent viscosity of the order v (∆)∆ ν. As a result, the purely non-linear turbulent dynamics of the "large eddies" is separated from microscopic dissipation. 1 The biggest challenge when implementing this concept is to find an appropriate model for the coupling between the small-and large-scale dynamics.
A mathematical framework for LES is based on the notion of a filter, which separates large-scale ( ∆) from small-scale ( ∆) fluctuations. Filters can be used to decompose the equations of fluid dynamics into equations for smoothed variables, which have a very similar mathematical structure as the unfiltered equations, and equations for second-order moments of the fluctuations. The latter are interpreted as subgrid-scale variables. In Section 2, we will carry out the decomposition of the compressible Navier-Stokes equation by applying the filter formalism of Germano [40]. This formalism comprises the so-called Reynolds-averaged Navier-Stokes (RANS) equations as limiting case if the filter length is comparable to the integral length scale of the flow. Simulations based on the RANS equations work with low Re eff , while LES have high Re eff . In principle, second-order moments can be expressed in terms of higher-order moments. Since this would entail an infinite hierarchy of moments, the set of variables is limited by introducing closures. Usually, one attempts to find closures for the second-order moments by expressing them in terms of the filtered variables. This is what is called a subgrid-scale (SGS) model. 2 For example, a complete second-order closure model for turbulent convection is formulated in [19]. Much simpler, yet often employed is the one-equation model for the SGS turbulence energy K, i. e., the local kinetic energy of numerically unresolved turbulent eddies. For this reason, it is sometimes called the K-equation model. Closures for the transport and source terms in the SGS turbulence energy equation are presented in some detail in Section 3, followed by a discussion of how the closure coefficients can be determined (Section 4). Of particular importance is the prediction of the local turbulent viscosity, which is is given by ∆ √ K times a dimension-less coefficient. The turbulent viscosity is required to calculate the turbulent stresses, which enter the equations for the filtered variables analogous to the viscous stresses in the unfiltered Navier-Stokes equations (see Section 3.1).
Filtering the dynamical equations is usually considered to be equivalent to numerical discretization. The filter length can then be identified with the grid scale ∆. Since the numerical truncation errors of stable finite difference or finite volume schemes are more or less diffusion-like terms, they produce a numerical viscosity that effectively reduces the Reynolds number to a value comparable to equation (3). It is actually a common assumption that numerical viscosity approximates the turbulent viscosity on the grid scale. This leads to the notion of an implicit large eddy simulation (ILES) [39], which is widely used for simulating turbulent flows in astrophysics. Numerous numerical studies demonstrated that ILES is a very robust method, which reliably predicts scaling laws of compressible turbulence at sufficiently high resolution [124,60,4,113,33,61]. This is a consequence of the independence of inertial-range scaling from the dissipation mechanism, be it microscopic, turbulent or numerical viscosity, provided that the dynamical range of the simulation is large enough. In simulations of statistically stationary isotropic turbulence, however, the iner-tial subrange is very narrow for computationally feasible resolutions because the bottleneck effect distorts the spectrum over a large range of high wave numbers below the Nyquist wavenumber [30,27,114]. It appears that LES with an explicit SGS model, such as the K-equation model, can reduce the bottleneck effect to some degree and reproduce scalings from ILES or DNS at lower resolution [45,128,109]. However, more systematic studies covering the parameters space of forced compressible turbulence are necessary to confirm this effect.
There are, of course, alternative methods of scale separation and a large variety of SGS models (for a comprehensive overview, see the monographs [103,39]). An example are the Camassa-Holm equations, which follow from the incompressible Navier-Stokes equations by decomposing the trajectories of fluid elements into mean and fluctuating parts in the Lagrangian framework [21]. Since the filtered component of the velocity is defined by an inverse Helmholtz operator of the form (1 − α 2 ∇ 2 ) −1 , which is explicitly applied to determine the turbulent stresses in the filtered velocity equation, the resulting model is called Lagrangian-averaged Navier-Stokes α-model (LANS-α). Depending on the choice of α, the variables computed in LES based on LANS-α are typically smoothed over length scales somewhat larger than the grid resolution. In other words, this type of simulation partially resolves the sub-filter scales, which improves the controllability of the model. While there is no handle on the competition between the SGS model and numerical truncation errors on the grid scale in convectional LES, LANS-α can, in principle, alleviate this problem by adjusting the balance between truncation and model errors [89]. Although the idea is very elegant, the numerical studies discussed in [89,90] show that the applicability of LANSα and similar models is limited, particularly for very high Re. Moreover, the generalization to compressible turbulence is not straightforward. Models such as LANS-α are not further covered by this review, but they might be an option for magnetohydrodynamical LES [91].
Currently, LES are mainly applied to complex astrophysical systems. In simulations of cosmological structure formation, which are discussed in Section 6.3, the length scales on which turbulence is driven by gravity are varying. Although adaptive mesh refinement is applied to track down collapsing structures, it is difficult to to resolve a wide range of length scales between the smallest driving scale and the the grid scale at the highest refinement level. In this situation, SGS effects can become fairly large. However, the variable grid scale complicates the scale separation in AMR simulations because energy has to be transferred between the resolved and SGS energy variables if a region is refined or de-refined. Section 5.1 describes how to combine LES and AMR. This method, for which the acronym FEARLESS (Fluid mEchanics for Adaptively Refined Large Eddy SimulationS) was coined in [74], has been applied to galaxy clusters, the intergalactic medium, and primordial atomic cooling halos. The results from these simulations indicate that the contribution of the numerically unresolved turbulent pressure to the support against gravity is non-negligible and the turbulent viscosity tends to stabilize disk-like structures around collapsed gas clouds. Moreover, the SGS model provides indicators of turbulence production and dissipation and allows for the computation of the turbulent velocity dispersion. A difficulty is that turbulence production by cosmological structure formation is highly inhomogeneous. This entails the problem that the SGS model should dynamically adapt to conditions ranging from laminar flow to developed turbulence. Inhomogeneous and non-stationary turbulence can be treated by dynamical procedures for the calculation of closure coefficients or shear-improved SGS models, which decompose the numerically resolved flow into mean and fluctuating components. These techniques are outlined in Sections 4.2 and 5.2.
Furthermore, SGS models offer unique possibilities for modeling physical processes that are influenced by turbulence. An example is turbulent deflagration, where the turbulent diffusivity predicted by the SGS model dominates the effective flame propagation speed in underresolved numerical simulations. Turbulent deflagration plays a role at least in the initial phase of thermonuclear explosions of white dwarfs (see Section 6.1), which is one of the scenarios that are thought to produce type Ia supernovae. A recent application along similar lines are LES of iso-lated disk galaxies, where the SGS turbulence energy is a crucial parameter for calculating the star formation rate and the feedback due to supernova blast wave (see Section 6.2). Since the impact of feedback processes on the formation of galaxies and their evolution leaves many questions unanswered, galaxies are a particularly promising field of application.
While great progress has been made for compressible hydrodynamics, magnetohydrodynamical LES are still in their infancy. Several SGS models have been proposed in the context of terrestrial plasma physics [77,78,45,22,91,122], but their applicability to astrophysical plasmas is unclear. Astrophysical MHD turbulence, particularly in the interstellar medium, extends to the supersonic and super-Alfvénic regimes. Moreover, plasmas become collisionless for high temperatures and low densities. A typical example is the solar corona. It is also likely to be the case in the intracluster medium. Since the fluid-dynamical description is not applicable in this case, kinetic methods have to be employed. Nevertheless, MHD-LES could provide a reasonable approximation on length scales that are sufficiently large compared to the characteristic scales of kinetic processes. In any case, SGS models for MHD turbulence will be a very challenging problem because of the local anisotropy of turbulent fluctuations, the potentially strong back-reaction from smaller to larger scales, and complicated dissipative processes such as turbulent reconnection [10,16,134]. In this area, extensive fundamental studies will be necessary.

Scale Separation
Large eddy simulations are based on the notion of scale separation. Although turbulence is a multiscale phenomenon, with interactions among different length scales, a separation into smoothed and fluctuating components can be rigorously defined by means of filter operators. Of course, the filtering of non-linear terms gives rise to interactions between these components. Filter operators were originally applied in the context of mean-field theories, but can be generalized to LES. For incompressible hydrodynamical turbulence, Germano [40] introduced a general framework that encompasses mean field theories as limiting case.
The smoothed component of a generic field variable q(x, t) is defined by means of a spatial low-pass filter, which is a convolution of q with an appropriate filter kernel G (see Chapter 2 in [103]): A homogeneous isotropic low-pass filter has the following properties: • The filter kernel is independent of direction: • Filtering smoothes out fluctuations on length scales smaller than the filter length ∆ G . Length scales that are large in comparison to ∆ G are not affected. This implies • The filter operator is linear, conserves constants, and commutes with spatial derivatives: The simplest low-pass filter is the box or top-hat filter. For Cartesian coordinates x i , the kernel of the box filter is defined by Usually, ∆ i is assumed to be equal for all spatial dimensions. The mean value of q in a rectangular domain with periodic boundary conditions follows in the limit that ∆ i is the linear size of the domain in each dimension. The construction of a filter is particularly simple in Fourier space. For a low-pass filter, the Fourier transform of the filter kernel, the so-called transfer function G(k), drops rapidly to zero for wavenumbers k k c = π/∆ G . Since only the Fourier modesq(k, t) with k k c contribute significantly to the corresponding filtered field q G (x) in physical space. The simplest case is the sharp cutoff filter, for which The sharp cutoff filter, however, is not equivalent to the box filter, which has the Fourier representation A filter that is intermediate between these two cases is the Gaussian filter.

Decomposition of the compressible Navier-Stokes equations
The compressible Navier-Stokes equations for the mass density ρ, the momentum density ρu, and the energy density ρE of a neutral fluid subject to gravitational and mechanical accelerations g and f , respectively, are Thermal conduction is neglected here. The energy per unit mass can be expressed as where e is the internal or thermal gas energy. For a perfect gas, e is related to the gas pressure P and the temperature T via the ideal gas law: where γ is the adiabatic exponent, k B the Boltzmann constant, µ the mean molecular weight, and m H the mass of the hydrogen atom. The viscous stress tensor σ is defined by where the two coefficients η and ζ are the dynamic and bulk viscosities of the fluid, is the rate-of-strain tensor, and the trace S ii is equal to the divergence d = ∇·u. The gravitational acceleration is given by g = −∇φ, where the gravitational potential φ is determined by the Poisson equation for a constant background density ρ 0 (G is Newton's constant). Mean-field equations for compressible turbulence are derived in [20]. Much in the same way, a general low-pass filter G can be applied to the system of PDEs (9)- (11). Alternative formulations can be found in [39], Section 2.4. For brevity, we omit the subscript G in the following. Since commutes with differential operators, the smoothed mass density ρ obeys an equation of exactly the same form as the continuity equation, if we set ρu = ρ ũ. This identiy implies the definition of the Favre-filtered velocitỹ 7 Filtering the momentum equation results in Owing to the non-linearities, however, we are facing some difficulties here. To obtain a PDE with the same basic structure as the unfiltered momentum equation, the advection term on the left-hand side should read ∇ · [ ρ ũ ⊗ũ]. The solution is to split the filtered non-linear terms: Since the Poisson equation (16) is linear, the smoothed potential φ is solely determined by ρ . The self-gravity term ρ∇φ , however, has to be split by defining The specific force f , on the other hand, usually varies only over the largest scales of the system. If the filter length is small compared to these scales, ρf ρ f is a good approximation. Thus, the filtered momentum equation can be casted into the following form [39,76,131,40,20,116]: Now, what is the physical interpretation of the terms τ and γ? Let us first consider the weakly compressible limit. By assuming that ρ varies only little over the filter length, density factors can be pulled out of brackets. In this case,ũ u . By defining the fluctuation of the velocity as u = u −ũ, it follows that If we futher assume that is a Reynolds operator (see Section 3.3 in [103]), which is not generally true for filters but applies, for example, to global averages, filtered quantities can be pulled out of brackets and the above expression simplifies to Although this simple relation holds only for a Reynolds operator in the weakly compressible limit, τ is generally interpreted as the stress tensor associated with the turbulent velocity fluctuations below the filter length. For this reason, τ is called the subgrid-scale turbulence stress tensor in the context of LES. The non-linear interactions of the filtered flow (the "large eddies") with smallscale fluctuations below the grid scale ∆ are given by ∇ · τ in equation (21). Likewise, the term γ defined by equation (20) accounts for the momentum transfer due self-gravitating fluctuations in the density. The trace of τ defines the fraction of kinetic energy on length scales smaller than the filter length: If the filter length is the grid scale, ρK is called the subgrid-scale turbulence energy. The first term on the right-hand side of equation (22) is the total kinetic energy, the second term the kinetic energy on length scales greater than the filter length (i. e., the numerically resolved kinetic energy in LES).
In the limit of high Reynolds numbers, the viscous dissipation scale (also known as Kolmogorov scale) is typically much smaller than the filter length. In this case, scaling arguments for incompressible turbulence imply σ τ , i. e., the filtered viscous stresses are negligible compared to the stresses associated with the turbulent velocity fluctuations [102]. Since the scaling of compressible turbulence tends to be stiffer than for incompressible turbulence [60,113,112], one can reasonably assume that this conclusion is generally applicable. The filtered momentum equation (21) thus can be written as where τ * is the trace-free part of τ : As one can see from equation (23), the trace of τ is associated with the turbulent pressure 2 3 ρK at the filter length scale.
In contrast to the filtered momentum density, which can be expressed as ρu = ρ ũ, the energy density on length scales greater than the filter length is given by where the second equality follows from equations (12) and (22). Consequently, ρ Ẽ = ρE . 3 A PDE for ρ Ẽ follows form the contraction of equation (23) withũ plus the filtered internal energy equation. The subtraction of this PDE from the filtered equation for the total energy yields the PDE for ρK (see Section 3.3 in [103] and [76,131,40,20,116]). In the limit of high Reynolds numbers, the resulting equations are: The additional source and transport terms resulting form the scale separation of the energy are defined as follows.
• Gravitational energy injection on subgrid scales: • Rate of subgrid-scale turbulence energy production: 4 where τ ij is defined by equations (19) andS ij is the rate-of-strain tensor associated with the Favre-filtered velocity: 5S 3 In [39], ρ Ẽ is identified with ρE and a different symbol is used forẽ + 1 2ũ 2 . However, we do not follow this nomenclature here. 4 Also called turbulence energy flux, although this is not a transport term. In the incompressible limit, Σ corresponds to the energy transfer in spectral space. 5 The definition ofS ij is a consequence of integration by parts ofũ i ∂ j τ ij . The symbolS ij is used for convenience.
• Rate of viscous energy dissipation in the limit of high Reynolds numbers: 6 where S ij is defined by equation (15), |S * | 2 = 2S * ij S * ij is the squared norm of the trace-free rate-of-strain tensor S * ij = S ij − 1 3 dδ ij and d = S ii . Although the viscous stresses can be neglected in the filtered momentum equation, viscous dissipation is crucial for the energy balance of turbulent flows.
• The flux of turbulent energy diffusion on sub-grid scales: For Reynolds operators in the weakly compressible limit, F (kin) can be expressed as a thirdorder moment of the velocity fluctuation: 2F (kin) j −ρ u i u i u j [40].
• There is also a viscous flux, which can be neglected relative to other flux terms if the Reynolds number is sufficiently high.
By adding the equations (26) and (27), we obtain an equation for the filtered total energy Except for the gravitational source terms, production and dissipation rates cancel out. The fluxes on the very right are related to turbulent transport processes below filter length. In particular, the sum of F (conv) and F (press) can be expressed as convective enthalpy flux, where ρh = ρe + P , which corresponds to −ρ u h in the weakly compressible limit. For a closed system of PDEs, it is necessary to compute all terms defined above in terms of known quantities. A rigorous calculation requires further PDEs, which involve higher-order moments and so on ad infinitum. This is known as the closure problem. A subgrid-scale model truncates the closure problem by approximating moments above a given order by lower-order moments. 6 Since S ij is a velocity derivative, it is of the order of the velocity fluctuation at the smallest length scales.

Cosmological fluid dynamics
In cosmological simulations, the equations of fluid dynamics are solved in a comoving coordinate system. Coordinates of observers that are stationary relative to the Hubble expansion of the Universe are constant in this system. The expansion is characterized by the scale factor a(t), which is determined by the Friedmann equations for a homogeneous and isotropic cosmology [87]. If the proper coordinates, which include changes of position due to the expansion of the Universe, are denoted by x proper and t proper , the corresponding comoving coordinates are x = x proper /a and t = t proper . Derivative operators transform as Furthermore, the invariance of mass implies that the comoving baryonic density ρ is related to the proper density by ρ = a 3 ρ proper . It can then be shown that continuity equation for ρ in comoving coordinates assumes exactly the same form as equation (9): Here, u is the so-called peculiar velocity, which is defined as where u proper =ẋ proper is the proper velocity and H =ȧ/a the Hubble constant. This means that, in the comoving coordinate system, matter moves with velocity u relative to the Hubble flow Hx.
With some algebra, also the momentum and energy equations can be transformed to comoving coordinates. The resulting equations do not have the same form as equations (10) and (11), but include additional terms with prefactors H. However, a particularly simple representation of the momentum and energy equations is obtained if the proper peculiar velocity is used in place of u. Filtered dynamical equations for cosmological fluids were first derived in [73,74] and presented in an alternative formulation in [118]. The applied filter kernel is static in comoving coordinates, i. e., the filter length increases proportional to the cosmological scale factor a. Consequently, commutation of the filter with time derivatives is unaffected by the cosmological expansion and equations for filtered dynamical variables follow completely analogous to Section 2.1. By neglecting gravitational terms associated with fluctuations below the filter length, the following equations for the filtered mass density ρ , the filtered momentum density ρU = ρ Ũ , and energy density ρ Ẽ , whereẼ =ẽ + 1 2Ũ 2 , are obtained: Here, the filtered internal energy density is ρ ẽ = ρe = P /(γ − 1), where P = a 3 P proper , and the gravitational potential φ of baryonic and dark matter density fluctuations is given by the cosmological Poisson equation where H 0 =ȧ(t 0 ) is the Hubble constant and Ω m (t 0 ) the density parameter of matter at redshift zero. Since the mean matter density ρ m,0 is constant in comoving coordinates, the source term of the Poisson equation can expressed in terms of the density fluctuation δ m = ( ρ dm +ρ −ρ m,0 )/ρ m,0 for the local dark matter density ρ dm and baryonic mass density ρ. The density parameter is defined by Ω m (t 0 ) = ρ m,0 /ρ crit,0 , where ρ crit,0 = 3H 2 0 /(8πG) is the critical density at time t = t 0 . The kinetic energy associated with peculiar velocity fluctuations below the filter length, is given by the dynamical equation Since the prefactors of all terms except for the time derivatives in the momentum and energy equations are unity and a, respectively, the definitions of all source and transport terms in equations (41), (42) and (44) are analogous to the definitions given in Section 2.1, with u i being replaced by U i . This suggests that closures for turbulence in a static space are applicable to cosmological fluids as well. Although cosmological expansion, in principle, causes a dampening of the kinetic energy [118], this effect is subdominant for turbulent eddies even on the largest scales in galaxy clusters because turbulence is driven in gravitationally bound gas on time scales shorter than the current Hubble time 1/H 0 .

Subgrid-Scale Models
There is a beautiful correspondence between finite-volume discretization and filtering. Finitevolume methods solve an equation for the cell averages of some dynamical variable q(x): Here, (x i , y j , z k ) are the cell-centered coordinates and ∆ is the linear size of a grid cell. It is not difficult to see that Q ijk equals the box-filtered variable q G for a box filter G with filter length ∆ at the discrete points (x i , y j , z k ). 8 Also finite differences correspond to low-pass filters. Thus, numerical discretization can be interpreted as implicit filtering. The numerical errors in the approximations to Q ijk can be characterized by the truncation error of the finite-volume method. For stable schemes with flux limiters, these errors are usually associated with terms of the diffusion type, i. e., proportional to ∇ 2 q. 9 In ILES, this is what causes the dissipation of kinetic energy into heat (a detailed account of ILES is given in [39]). A similar expression follows if the Boussinesq expression for the turbulent stresses is used as explicit SGS model for the interaction between numerically resolved and unresolved turbulent eddies. An important difference, however, is that the turbulent viscosity in the resulting diffusion terms in the momentum and energy equations is controlled by the dynamical variable K, which is the kinetic energy associated with turbulent velocity fluctuations below the grid scale. In this section, we mainly discuss the computation of K in LES. The approach we follow here is known as functional modeling. The aim is to model only statistical effects of SGS turbulence on the dynamics of the filtered fields, which are identified with the numerical solution. An alternative strategy is structural modeling (see Chapter 5 in [39]), which is not covered here.

Closures for the turbulence stress tensor
The SGS turbulence stress tensor, which is associated with the non-linear energy transfer between large and small scales, is the central quantity that has to be modeled in LES. The most commonly used closure is the eddy-viscosity closure. The underlying assumption is that the form of the tracefree part τ * is analogous to the anisotropic viscous stress tensor σ * , with the correspondence 10 where S ij andS ij are defined by equations (15) and (30), respectively. The turbulent viscosity ν sgs is assumed to depend on the grid scale and the unresolved turbulent velocity fluctuation (see, for example, Section 4.3 in [103]): Hence, For brevity, we drop brackets and tildes indicating filtered and Favre-filtered quantities from now onwards, so that the turbulent stresses can be written as In the following, it is understood that all quantities are either numerically resolved variables or modeled in terms of these variables. The production rate (energy flux) corresponding to the eddyviscosity closure is where d = S ii is the divergence. The eddy-viscosity coefficient C ν is typically in the range from 0.05 and 0.1 [103,116,111].
In the incompressible case (d = 0), the eddy-viscosity closure admits only positive energy flux. However, direct numerical simulation data show that there is a certain amount of backscattering from smaller to larger scales, corresponding to a negative energy flux [116,111]. This motivated a closure for the turbulent viscosity that is constructed from the determinant of the trace-free rate-of-strain tensor [129]: 11 By substituting the above expression into equation (49) for the turbulent stresses, it follows that the production rate is given by Since the determinant can be positive under certain flow conditions, in principle, this closure accounts for backscattering (also known as inverse cascade). This phenomenon can be explained by the so-called the "tornado" topology, i. e., the alignment of vortices along a single stretching direction [129]. Then the flow is contracting in one dimension and expanding in the other two, which results in a positive determinant. A negative determinant, on the other hand, corresponds to the standard situation of a forward cascade transporting energy from larger to smaller eddies. While the determinant closure modifies only the turbulent viscosity, a different expression for the SGS turbulence stress tensor is proposed for compressible turbulence in [128]. Based on Taylor series expansions of the velocity around grid cell centers, an appropriate normalization leads to the non-linear closure τ where |∇ ⊗ u| = (2u i,k u i,k ) 1/2 is the norm of the velocity derivative. The above expression satisfies the identity τ ii = −2ρK. However, it is generally not adequate as a model for the turbulence stress tensor in LES [111]. In contrast to the eddy-viscosity closure, rotation invariance is violated because of the antisymmetric part of ∇ ⊗ u. This would cause spurious production of K in a uniformly rotating fluid. A further problem is that K = 0 would be a fixed point of equation (27) if all other sources of turbulence energy are zero. This results in unphysical behavior. With the eddy-viscosity closure, on the other hand, K can grow sufficiently fast from arbitrarily small initial values because ν sgs is proportional to √ K rather than K. For this reason, a linear combination of τ (nonlin) ij and 2ν sgs S * ij is used in [129], where ν sgs is given by equation (51). The additional determinant term, with a small coefficient C 1 , has the function of a seed term that triggers the production of turbulence energy, while the production rate vanishes for a uniformly rotating fluid.
With the standard turbulent viscosity defined by equation (47), the same idea leads to the following generalized two-coefficient closure [111]: The coefficient C 2 determines the relative contributions from the non-linear and divergence terms to the trace τ ii . The purely non-linear closure corresponds to C 1 = 0 and C 2 = 1. Equation (49), on the other hand, is obtained if C 1 = C ν / √ 2 and C 2 = 0. For the application in LES, it is necessary to calibrate the closure coefficients C 1 and C 2 . For supersonic turbulence, C 1 = 0.02 and C 2 = 0.7 appear to be robust values (see Section 4). The rate of production following from the generalized closure is The first term dominates if K 1/2 is small compared to ∆|S * |. For strong turbulence intensity, i. e., K 1/2 ∆|∇ ⊗ u|, the second term contributes significantly. The transition is further influenced by the ratio C 2 /C 1 .

The Sarkar-Smagorinsky model for weakly compressible turbulence
In the case of isotropic incompressible turbulence, the mean SGS turbulence energy K for a sharp cutoff at the length scale ∆ is obtained by integrating the Kolmogorov spectrum E(k) over wavenumbers k ≥ π/∆: The mean dissipation rate is therefore given bȳ for the Kolmogorov constant C ≈ 1.65 [94]. It is commonly assumed that an expression of this form also holds for the local dissipation rate in LES (see, for example, [103]): with C ∼ 1. The above dimensional closure for the dissipation rate basically means that the time scale of energy dissipation is given by τ ∼ ∆/ √ K. For subsonic compressible turbulence, closures for the dissipation rate and pressure dilatation are obtained by separating the pressure fluctuations into a rapid osciallatory and a slow component [104]. The resulting combined expression for + λ reads where is the turbulent Mach number associated with the SGS velocity fluctuation √ 2K. A particularly simple SGS model can be formulated by neglecting all gravitational and transport terms associated with subgrid-scale effects in equation (27). If furthermore a balance between production and dissipation is assumed, then Here, the eddy-viscosity closure (49) is substituted for τ * ij . If the above algebraic equation is solved for K, the PDEs (17), (23), and (26) form a closed system. The effect of the compressibility corrections is a reduction of the production due to anisotropic shear, ν sgs |S| * , by the factor (1 − α 2 M sgs ) and an enhancement of the solenoidal dissipation rate C K 3/2 /∆ by a factor that increaes with the square of M sgs (α 1 tends to be greater than α 3 ). An extension to a non-equilibrium model based on the dynamical equation (27) for K was exploited in [73,74] for cosmological LES of the gas in galaxy clusters. However, this model applies only if M sgs is small compared to unity, which is the case for turbulence in the intracluster medium. On the other hand, the correction factors are close to unity for small M sgs and, given the many approximations involved, it is not clear whether they have any significant effect. Apart from that, the model definitely breaks down in the vicinity of accretion shocks and in the cooler regions of the intergalactic medium, where M sgs can become large compared to unity.
In the limit M sgs → 0 and d → 0, the classical Smagorinsky model for incompressible turbulence [121] follows from equation (61). In this case, where This expression has an important implication. One could calculate from ILES data analogous to the viscous dissipation rate following from the Navier-Stokes equations, i. e., Here, ν eff = V L/Re eff is assumed to be the constant numerical viscosity, which is given by the effective Reynolds number Re eff of the simulation [86]. However, the above estimation of the dissipation rate is clearly at odds with equation (63). This can be understood as follows. If Re eff ∼ Re, where Re is the physical Reynolds number defined by equation (1), then ν|S| 2 is the physical dissipation rate in a direct numerical simulation of a fluid with microscopic viscosity ν. In an LES with Re eff Re, on the other hand, the Smagorinsky model implies a turbulent viscosity of the order ∆ 2 |S| for steady-state turbulence. The turbulent viscosity is not a constant. In this case, the dissipation rate is approximately given by equation (63). This is a consequence of equation (31), which implies that ρ cannot be expressed as the contraction of the filtered viscous stress tensor with the filtered rate-of-strain tensor. The dissipation rate is instead given by the filtered contraction of the two tensors. As shown in [111], the argument remains valid even if the dissipation rate is calculated for LES of supersonic turbulence with the advanced SGS model presented in the following section.

The compressible subgrid-scale turbulence energy model
To determine K, one can either invoke the equilibrium condition, such as in the Smagorinsky model, or numerically solve the PDE (27). The latter is called the SGS turbulence energy model [103,119,76,131,40,116,111]. If gravitational terms are negligible, the turbulence energy equation can be explicitly written as where the closure (55) for the production rate Σ, the dissipation rate defined by equation (58), and the gradient-diffusion closure for F (kin) + F (press) were substituted into equation (27). The gradient-diffusion hypothesis, which is also known as Kolmogorov-Prandtl relation, is based on the assumption that the turbulent transport of K is a diffusion process satisfying Fick's law (see, for example, [94,116]): with a turbulent diffusivity The Prandtl number of turbulent transport, C κ /C ν , is often assumed to be of the order unity. For a calibration of C κ , see Section 4. The pressure-dilatation λ is assumed to be negligible in equation (64) because no satisfactory closure is known for the highly compressible regime [128,111]. For weakly compressible turbulence, equation (59) could be used, but the applicability of this closure and an appropriate calibration of the coefficients α 1 , . . . , α 4 requires further investigation. As pointed out in Section 3.2, the contribution from λ is too small to significantly influence K in the weakly compressible regime. In this case, the equation with the eddy-viscosity closure, i. e., C 1 = C ν / √ 2 and C 2 = 0, can be regarded as a sufficient model for most applications. In particular, this variant of the SGS turbulence energy model was used for simulations of thermonuclear combustion in white dwarfs (Section 6.1) and cosmological simulations (Section 6.3).
Negligible pressure-dilatation is also a reasonable assumption at high Mach numbers because the kinetic energy is large compared to the internal energy and non-linear interactions between turbulent velocity fluctuations should be the dominant mode of energy transfer. This is manifest in the closure (50) for Σ, which is solely constructed from the gradient of the resolved velocity field, but does not depend on density or pressure gradients. There are both theoretical and numerical studies in support of this conjecture. In [2,3], it is argued that a range of length scales exists, in which the kinetic and internal energies decouple and the flux through the kinetic energy cascade becomes asymptotically constant, while ρλ is subdominant. The computation of the different contributions to the total energy flux at varying length scales from supersonic turbulence data in [61] confirms this conclusion. Based on an analytical theory for the two-point correlations of compressible turbulence [37], it is shown that the main contribution to the energy flux is where δu is the velocity difference between two points separated by a distance r, δu is the longitudinal component of the velocity difference in the direction of r, and the brackets denote the ensemble average. The closure (55) for the Σ has a similar structure, with factors of ρK 1/2 and derivatives of u corresponding to fluctuations on the grid scale. However, as pointed out in Section 2.5.2 of [39], a subtlety arises in the presence of shocks because the Rankine-Hugeniot conditions for jumps across shock fronts should be filtered in place of the PDEs. This entails SGS terms that are different from the terms in the filtered PDEs. However, it is questionable whether any attempt to model these terms would be useful. The assumption that the unmodified jump conditions apply to the numerical solution amounts to a fallback from LES with an explicit SGS model to ILES. Since shock-capturing schemes, such as PPM, fall back to stronger diffusion in the vicinity of shocks, this is probably the most reasonable thing one can do. Nevertheless, the closure (55) captures the non-linear interscale transfer of energy due to supersonic turbulent velocity fluctuations. The SGS model outlined above accounts for the statistical effect of shocks as well as vortices interacting with each other across the grid scale [113,112], while any SGS terms in the jump conditions would mainly correct geometric differences between smoothed shock fronts and the corresponding unfiltered fronts with substructure on smaller length scales (just like the turbulent flame fronts discussed in Section 6.1). Figure 1: Visualization of the SGS turbulence energy density ρK in a 512 3 LES with solenoidal forcing [111].
Indeed, equation (64) for ρK was demonstrated to work very well in the highly compressible regime [111]. As an example, Figure 1 shows a visualization of ρK from an LES of isotropic supersonic turbulence, where solenoidal stochastic forcing maintains a root mean square Mach number of about 5 in the statistically stationary regime. 12 The numerical resolution is 512 3 . In the reddish regions of the plot, K sgs is higher than the spatial average, while it is lower in the bluish regions. The structure of the numerically resolved turbulent flow is illustrated by the so-called denstrophy, in Figure 2.
Since Ω 1/2 combines density fluctuations and the rotation of the velocity, ∇ × u, it indicates both small-scale compression and eddy-like motion [60]. There is clearly a correlation between Ω 1/2 and ρK, which reflects the local interaction between resolved small-scale modes and subgrid-scale turbulence, as expressed by the production terms in equation (64). This correlation is akin to the equilibrium condition (62) following from the Smagorinsky model for incompressible turbulence. Owing to the non-local effects in the PDE (64), however, the SGS turbulence energy cannot be reliably estimated from local quantities such as Ω 1/2 [111]. In particular, turbulent diffusion smears out ρK in comparison to Ω 1/2 . Figure 2: Visualization of the the denstrophy Ω 1/2 for the same LES as in Figure 1.
A critical property is the scaling behavior of the SGS turbulence energy. For statistically stationary homogeneous turbulence, the mean value of ρK should scale as a power of the grid resolution because the fraction of unresolved kinetic energy changes as the the cutoff of the energy spectrum is shifted (see equation 56). This was verified in [111] by running LES with different grid scales ∆ and fixed forcing length L. The global spatial averages ρK in these simulations are plotted in Figure 3 (left panel) for ∆ ranging from L/256 to L/32, where the case ∆ = L/256 corresponds to the 512 3 simulation depicted in Figures 1 and 2. Although there are substantial fluctuations, one can qualitatively see that ρK decreases with ∆. Time-averaging over the statistically stationary regime yields mean values that are close to the power law with α ≈ 0.799 ± 0.009 (see Figure 4). The scaling exponent is in between the Kolmogorov and Burgers exponents and roughly comparable to the slope of the second-order structure functions with fractional mass-weighing reported in [113].
In the filtered momentum and energy equations (21) and (26), respectively, the trace of the SGS turbulence stress tensor acts as an additional turbulent pressure. The sum of the thermal and turbulent pressures is sometimes called the effective pressure: It is important to keep in mind that P eff depends on the numerical resolution and P eff → P in the limit ∆ → 0 (DNS). Figure 5 shows a phase plot of the effective pressure vs. the mass density for the highest-resolution case. One can see that the average of the effective pressure for a given   (69) vs. the mass density [111]. Both quantities are normalized by their mean values, P 0 and ρ 0 . The thin dashed line corresponds to the isothermal relation P ∝ ρ, and the dotted line indicates the mean effective pressure for given mass density.
mass density closely follows the isothermal relation P ∝ ρ. Although the mean turbulent pressure 2 3 ρK is small compared to the thermal pressure for the resolution ∆ = L/256, the intermittency of turbulence can locally produce an effective pressure that exceeds the thermal pressure by one order of magnitude. Consequently, the turbulent pressure can become important for compressible turbulence, particularly if there are other sources than the turbulent cascade. As an example, turbulent feedback in galaxy simulations is discussed in Section 6.2. In addition to the turbulent pressure, the non-diagonal turbulent stresses τ * ij act on the resolved flow. In the case of the eddy-viscosity closure, τ * ij causes a diffusion effect on top of the numerical diffusion, which occurs regardless of the compressibility of the flow. 13 For strongly diffusive numerical schemes, this effect is marginal. For high-resolution schemes, on the other hand, the explicit turbulent stresses in LES can become significant. Moreover, the non-linear term in equation (54) modifies the diffusion-like tensor in the case of supersonic turbulence. Of course, adding the turbulent stresses in LES does not merely degrade a high-resolution scheme to a more diffusive scheme because the diffusion is linked to the non-linear turbulent interactions across the grid scale. For non-turbulent flow, the turbulent stresses should vanish if the SGS model is consistent.

Two equation models and gravity
In the framework of the Renolds-averaged Navier-Stokes equations (RANS), which are equivalent to the filtered Navier-Stokes equation in the limit of a filter scale comparable the integral scale of the flow, the K-turbulence model can be used to calculate both the turbulence energy K and the dissipation rate . Two inhomogeneous PDEs of the advection-diffusion type determine K and [94]. In contrast to the simple dimensional closure (58) with a single coefficient of order unity, the diffusion and source terms in the equation for come with several additional closure coefficients. This type of model is commonly used for industrial and environmental flows. A two equation model of similar structure is proposed for buoyancy-driven flows in [26]. The model predicts the energy K and characteristic size L of the dominant eddies produced by Rayleigh-Taylor and Richtmeyer-Meshkov instabilities. The evolution of these variables is given by the following two equations (in notation adapted to this review): The production rate due to the Rayleigh-Taylor instability is basically given by Γ ∝ ρ √ 2K g eff . For each grid cell, the buoyant acceleration g eff is assumed to depend on the density contrast across the cell faces, the length scale L, and the components of gravity along the three coordinate axes. Apart from differences in the determination of g eff , the above expression for Γ is the same as in equation (128) for the enhanced production of turbulence at the interface between lowdensity and high-density material, which is applied to the propagation of turbulent flame fronts in thermonuclear supernovae (see Section 6.1).
The K-L model outlined above was adopted as an SGS model for simulations of turbulence driven by active galactic nuclei in galaxy clusters [105]. In this case, turbulence is thought to be stirred by hot bubbles rising due to their buoyancy in the ICM. These bubbles originate from the AGNs in the cluster. For this reason, production through the turbulent cascade is set to zero, i. e., τ * ij = 0. 14 In contrast to LES based on the consistent decomposition derived in Section 2.1, it follows from the very concept of the K-L model that K cannot be interpreted as the kinetic energy associated with velocity fluctuations below the grid scale. Since K is the kinetic energy associated with the dominant eddies driven by the RT instability on a length scale L, where L is a dynamical variable that can become larger than the grid scale, K generally encompasses some fraction of the numerically resolved turbulence energy on top of the SGS turbulence energy (if L falls below the grid resolution, on the other hand, K will represent only a fraction of the SGS turbulence energy). Consequently, the K-L model works as a hybrid model in these simulations, which stand somewhere between RANS and LES. As a result, the resolution-dependent small-scale structure of the RT-unstable bubbles is smeared out and only the coherent structure on large scales is captured in [105]. LES, on the other hand, resolve small-scale structure always down to the grid scale.
For the general case of self-gravitating turbulent gas, no satisfactory closure for Γ has been found yet. A conceptual difficulty is that the acceleration caused by gravity is genuinely anisotropic, while SGS models such as the turbulence energy model are based on local isotropy. The usual solution to this problem is to resolve the flow down to length scales that are not strongly affected by gravity. In AMR simulations, this is achieved by imposing a Truelove-like resolution criterion such that a sufficiently large ratio between the local Jeans length and the grid scale is maintained [125,34]. Since the density of gravitationally unstable gas would increase indefinitely, excess mass is usually dumped into sink particles at the highest refinement level [63,31,127]. Thereby, collapsing gas is decoupled form the numerically computed gas dynamics. In a certain sense, a sink particle is nothing but an SGS model for a self-gravitating overdense cloud that collapses down to scales below the minimal grid scale. Despite being a crude model, sinks particles are a reasonable approximation to collapsed clouds because they mainly interact through accretion (i. e., mass accumulation) with the numerically resolved gas dynamics. A more complex situation is encountered if the objects represented by the sink particles produce feedback onto the gas. An example is stellar feedback in galaxy simulations, which can be treated with the approach discussed in Section 6.2.

Determination of Closure Coefficients
One of the basic assumptions of the Kolmogorov theory is that turbulence is statistically selfsimilar in the inertial subrange (see, for example, [36]). With regard to subgrid-scale closures, the self-similarity of turbulence implies that dimensionless coefficients such as C ν in equation (47) should be independent of the chosen filter scale. This is not only a necessary condition for the feasability of LES, but it also allows for the calibration of closure coefficients by explicitly filtering turbulence data. Since closures do not exactly match SGS terms, an improved approximation can be achieved by so-called dynamical procedures, which estimate coefficients from properties of the numerically resolved flow under the assumption of local self-similarity.

Hierarchical filtering
As a formal framework, let us consider an infinite series of isotropic and time-independent filter operators n . Each filter is defined by a kernel G n (r) with filter length ∆ n (see Section 2). We shall assume that ∆ 0 ∼ L, where L is the integral length scale of the flow, and i. e., n for n = 0, 1, 2, . . . is a self-similar hierarchy of filters. Typical examples are the box filter defined by equation (5) or the Gaussian filter, which has the kernel [103,94] G n (r) = 6 π∆ 2 for isotropic filter lengths ∆ n = ∆ 0 /λ n . Since G n (r) → δ(r) in the limit n → ∞, ∞ is the idenity operator. Because filtering in physical space corresponds to a multiplication with the transfer function of the filter in Fourier space, it follows that For Gaussian filters, the validity of this approximation becomes immediately clear by calculating the product of the transfer functions: We can now apply the scale separation of the Navier-Stokes equations introduced in Section 2 at different levels of the filter hiearchy. In particular, the filtered density field at the n-th level is ρ n , and the Favre-filtered velocity is given bỹ By filtering twice at levels m and n, we obtaiñ If the n-th level is much coarser than the m-th level, the asymptotic relation (74) The turbulence stress tensor on the length scale ∆ n of the n-th filter is defined by The stress tensors for two filter levels m and n, where ∆ m < ∆ n , are related by the Germano identity (see Section 3.3.3 in [103] and [40,107]): The stress tensor associated with the double-filtered variables is defined by and is the Leonard stress tensor, which is associated with velocity fluctuations in the intermediate range of length scales ∆ m ≤ ≤ ∆ n . The Germano identity also holds for two arbitrary filters in the hierarchy. In the limit ∆ n ∆ m , the contribution from τ [m] n becomes negligible and where the second relation follows from equation (74). As a consequence, the turbulent stresses associated with the scale ∆ n are not sensitive to the flow structure on much smaller scales. In particular, it follows that In [107,116,111], Gaussian filters (see equation 73) are applied to data form ILES of forced compressible turbulence for the verification of closures. The following line of reasoning is of central importance for estimating closure coefficients from finite-resolution data. To begin with, let us assume that ρ and u are the physical density and velocity fields. Let us further assume that the implicit filter of the ILES correspond to the filter level m = I, i. e., ρ I andũ [I] represent the numerically computed density and velocity fields. Now, if the numerical data are coarse-grained by an explicit filter n in the inertial subrange, the turbulence stress tensor τ [I,n] defined by equation (82) can be calculated. But a closure for the turbulent stresses on the the length scale ∆ n applies to τ [n] , which is defined by equation (79). If ∆ n is reasonably large compared to ∆ I , however, one can make use of the approximation τ [n] τ [I,n] (see relation 83 for m = I). Owing to equation (78), the distinction between, on the one hand, the physical densities and velocities (or DNS data) and, on the other hand, the ILES data becomes immaterial. It is thus possible to calculate coarse-grained eddy-viscosity coefficients as: Of course, since the eddy-viscosity closure is not exact, the value of C ν varies. But the mean value turns out to be roughly 0.05 for different filter lengths and simulation parameters [107], in good agreement with other estimates in the literature. The same method was used to determine the coefficient C κ for the gradient-diffusion closure (65). The result C κ ≈ 0.4 implies a Prandtl number C κ /C ν around 10, contrary to the common assumption that the kinetic Prandtl number is of the order unity [103].

Dynamical procedures
Subgrid-scale models in their standard form apply to statistically stationary and isotropic turbulence. But turbulent flows in nature often deviate from this idealization: In terrestrial applications, flow inhomogeneities are inevitably caused by boundary conditions ("walls"). In astrophysics, one of the major energy sources is gravity. It causes matter to clump (galaxies and clusters) or to move under the action of central gravitational fields (stars), which produces inherently inhomogeneous and anisotropic flows. For example, turbulent convection in stars introduces a vertical anisotropy of the flow. Turbulence driven by violent energy release (supernovae) can also be highly inhomogeneous.
One of the solutions to this problem is to localize closures, i. e., to calculate local closure coefficients. This requires local estimators that take properties of the flow in some small region as input. Obviously, this works only if the size of this region is not significantly affected by the flow inhomogeneity on larger scales. In other words, the flow must be asymptotically homogeneous and isotropic at least on length scales of the order of the grid scale. In this case, a so-called test filter T can be applied in LES, with a filter length ∆ T that is a small multiple of the grid scale ∆. Test filters are usually implemented as discrete filters over several grid cells (see Section 2.3.2 in [39]). A multi-dimensional test filter can be composed as a succession of one-dimensional filters. 15 The test filter length ∆ T can be adjusted by varying the weights of the cells. An optimal ratio γ T = ∆ T /∆ is given by the closest match between the filter transfer functions of the discrete and analytical box filters with filter length ∆ T [126]. For instance, a test filter with γ T = 2.771 is optimal if a five-point stencil is used in each spatial dimension [107].
By identifying ∆ m with ∆ and ∆ n with ∆ T , the Germano identity (80) allows us to express the turbulence stress tensor on the length scale of the test filter as the sum of the test-filtered SGS turbulence stress tensor and the Leonard tensor for the intermediate velocity fluctuations (see also Section 4.3 in [103]): Here, the Leonard tensor L associated with the test filter is defined by where we use the simplified notation ρ and u for the density and velocity on the grid scale, as in Section 3.1. Because of the scale-invariance of turbulence, [41] proposed that the eddy-viscosity closure (49) holds for both T and τ . In the case of the Smagorinsky model (see equation 62 for ν sgs ), the corresponding tensors are: The rate-of-strain tensor (S T ) ij at the test filter level is given by the symmetrized derivative of the test-filtered numerically resolved velocity field, ∂ i u j T , analogous to equation (30). The variable C S (x, t) needs to be determined. This can be achieved by substituting the above expressions for τ * ij and T * ij into the trace-free part of the Germano identity (84), which implies Under the assumption that C S varies only little over the smoothing length of the test filter, one can set C 2 S β ij T C 2 S β ij T . Since L ij can be evaluated from equation (85), minimalization of the residual error between L * ij and the expression on the right-hand side of equation (88) yields where m ij = α ij − β ij T . This is the Germano-Lilly dynamical procedure, which was applied, for example, in LES of turbulent channel flows [92]. In principle, this procedure could also be applied to the non-equilibrium model with the turbulent viscosity defined by equation (47). In this case, the turbulence energy associated with T is given by the contracted Germano identity, where ρK = −τ ii /2 and ρ T K T = −L ii /2. However, the dynamical procedure as outlined above has several caveats. In particular, the assumption of negligible variation of C S over the the test filter length is found to be violated significantly. Moreover, C S diverges if m ij vanishes. Consequently, several attempts were made to improve the dynamical procedure [71,93,42]. A particularly simple modification was found by analyzing experimental measurements of turbulent velocity fluctuations in consecutive wave number bands [k n−1 , k n ], corresponding to a hierarchy of filters. By explicitly evaluating the turbulent stresses τ [n] associated with the wave numbers k n = π/∆ n , the correlations with localized closures were verified. Although some correlation between the turbulent stresses at different filter levels was found, the correlation of τ [n][n−1] with the Leonard stresses L [n,n−1] turned out to be significantly better. This observation can be understood as a consequence of the locality of the energy transfer [58,103], i. e., the energy transfer across a certain wave number k is mainly caused by interactions in the narrow spectral band [ 1 2 k, 2k]. With regard to test filtering in LES, this implies that the eddy-viscosity closure should be applied to L in place of T. The localized coefficient C ν (x, t) of the turbulent viscosity is then given by [57,107,116,102] where K T = −L ii /(2ρ T ) is the resolved kinetic energy on length scales ∆ ≤ ≤ ∆ T . Substitution of the above expression for C ν into equation (50) for the localized rate of production yields The above formula was used for simulations of turbulent thermonuclear combustion in white dwarfs (see Section 6.1). A generalization of the dynamical procedure to localize both C 1 and C 2 in the closure (54) would be straight-forward, but has not been applied so far. In this case, a linear system in the coefficients C 1 and C 2 has to be solved to minimalize the residual. For ILES of subsonic and transonic turbulence produced by stochastic forcing in periodic boxes [114], the enhanced fidelity of the localized closure (91) can be verified by coarse-graining the data with hierarchical Gaussian filters n as explained in Section (4.1) [116]. Let us first consider the turbulence energy flux produced by anisotropic shear on the length scale ∆ n , assuming a turbulent viscosity with a constant coefficient C ν , which is obtained by averaging equation (90) over the domain of the flow. If the closure with this coefficient were exact, we would have Here, the divergence term is added to Σ [n](eddy) to express the flux associated with the tracefree rate-of-strain tensor S [n] * . However, coarse-grained data show that equation (92) is not very well satisfied. The probability density functions of the expression on the right-hand side and the explicitly calculated energy flux on the left-hand side are compared in the left plot in Figure 6.
One can see that the latter is negative in about 20 % of the domain (purple line), corresponding to backscattering from smaller to larger scales. This is excluded by the eddy-viscosity closure with a fixed coefficient (light blue line). In this case, negative values of the total energy flux Σ [n](eddy) are solely due to the divergence term. The bias toward positive fluxes can be avoided by localizing the eddy-viscosity closure [116]. For a test filter n−1 with filter length γ T = ∆ n−1 /∆ n > 1, the energy flux is given by the following analogue of equation (91): The probability density functions that are plotted for different test filtering ratios γ T in Figure 6 (left plot) indicate a substantially improved match between the localized closure and the explicitly calculated energy flux. Indeed, distributions of the localized closure coefficients show that C ν has a negative branch (see right plot in Figure 6). This result suggests an improvement due to the dynamical procedure even in the case of homogeneous turbulence. On the other hand, the mean values of C ν appear to be fairly robust for different forcing parameters.

Global least squares method
Closures can also be tested by analyzing correlations. This allows for the calibration of the closure coefficients by least squares minimalization of the integrated residual [111]. For example, let us consider a generic closure with a single coefficient C 1 at the n-th filter level: The global residual of the explicitly computed turbulence energy flux Σ [n] = τ ij is defined by equation (79), can be quantified by the squared error integrated over the whole domain V of the turbulent flow: The minimum of err 2 (C 1 ) is obtained by setting the derivative with respect to C 1 equal to zero: In contrast to the dynamical procedure, the resulting closure coefficients are constants. The method of least squares described above is applied in [111] to various ILES of supersonic isothermal turbulence produced by stochastic forcing [112,33]. To coarse-grain the data, a Gaussian filter with a smoothing length ∆ 4 = L/16 = 32∆ I is used, where ∆ I = ∆ 9 ≡ ∆ is the grid resolution of the ILEs. 16 Since the filter length is large compared to the grid resolution in this case, it is advantageous to apply the filter operation in Fourier space. For the eddy-viscosity closure, we have The closure coefficients following from equation (96) are, for instance, C 1 ≈ 0.102 for the 1024 3 ILES with purely solenoidal (divergence-free) forcing, and C 1 ≈ 0.092 in the case of compressive (rotation-free) forcing. The corresponding value of C ν = √ 2C 1 is about 0.14. When comparing this value to the results in [107,116] (see Section 4.1), one has to bear in mind that not only a different method is applied to determine the coefficients, but also that the turbulence properties differ substantially.
The correlation diagram for Σ [4](cls) , with the least-squares coefficient C 1 , versus Σ [4] in the case of solenoidal forcing is shown in Figure 7 where std denotes the standard deviation and the angle brackets indicate an average over the whole domain. The correlation coefficients of the eddy-viscosity closure are found to be 0.95 and 0.93 for solenoidal and compressive forcing, respectively [111]. However, it becomes apparent that the closure breaks down for negative fluxes. This corresponds to the bias of the probability density function for the static closure with an averaged coefficient in Figure 6 (left plot). Rather than applying the dynamical procedure, it is shown in [111] that the determinant closure (52) results in a largely improved approximation of negative turbulence energy flux. This is a consequence of the varying sign of the determinant, det S * , while ∆K 1/2 |S * | 2 is positive. However, the scatter of the determinant closure is high, particularly for large positive flux. This is remedied by the non-linear closure (53) for the turbulence stress tensor, which produces an excellent correlation between Σ [n](cls) and Σ [n] , as demonstrated by the right plot in Figure 7. Since the formulae are analogous to the eddy-viscosity closure, we refer to [111] for details. The correlation coefficients are 0.991 for both solenoidal and compressive forcing.
However, as explained in Section 3.1, the purely non-linear closure is not suitable for an SGS model. This is why the least squares method was applied to the generalized closure with two coefficient, C 1 and C 2 . For  and non-linear (right) closures for supersonic isothermal turbulence produced by solenoidal forcing [111]. The applied filter length is 32∆, where ∆ is the grid resolution. The average prediction of the closure for each bin is indicated by the blue dots.
the closure coefficients are given by the linear system of equations where The solution is C 1 ≈ 0.02 and C 2 ≈ 0.7, with a correlation coefficient 0.990 [111]. These coefficients also yield good approximations to the turbulence energy flux for isothermal and adiabatic turbulence simulations at lower Mach numbers [112,111]. Moreover, the coefficients appear to vary only little with the filter length, at least in the range from ∆ 3 = 64∆ to ∆ 5 = 16∆. Choosing higher or lower filter lengths is not sensible because of the influence of the forcing (∆ n must be small compared to ∆ 0 = L) and numerical dissipation (∆ n ∆).

Adaptive Methods
The most powerful technique for finite-volume codes to resolve localized and anisotropic structures in a flow is adaptive mesh refinement (AMR) [6,5]. Even with AMR, however, it is generally not possible to fully resolve turbulence. This entails the problem that the numerically resolved and unresolved turbulence energy fractions vary as regions are refined or de-refined. In the following, it shown how to address this problem in adaptively refined LES. In principle, global energy and momentum conservation can be achieved, while reducing the need for artificial changes in the internal energy, which is the standard method to restore energy consistency between different refinement levels in AMR simulations. Apart front that, localized and anisotropic flow structures pose the problem that SGS models with constant coefficients introduce systematic errors because they are usually calibrated for statistically stationary and isotropic turbulence. Shear-improved SGS models alleviate this problem by adjusting the non-linear energy transfer across the grid scale to local flow conditions. This is possible by applying an adaptive temporal filter, the so-called Kalman filter.

Energy-and momentum conservation in AMR simulations
In AMR simulations, data have to be transferred between different refinement levels by conservative interpolation or averaging. For example, if a region is refined, data from coarser grids are interpolated to finer grids. The same operation is used for filling ghost cells at the boundaries between a finer and a coarser level, which is required to compute fluxes through the faces of adjacent finer and coarser cells. Moreover, block-structured AMR codes usually average down the data from the highest-level grid to the all coarser levels. The mass density, momentum, and energy variables at two levels, say, l and l + 1, are in the simplest case related by where N grid cells at level l+1 are summed up to a single value in a coarse cell at level l. Obviously, these relations guarantee mass, momentum, and energy conservation. For more sophisticated interpolation schemes, the fine-grid values have different weights w n in the above sums. Without loss of generality, we assume w n = 1 in the following. Now, if we work out the internal energy at the coarser level, we obtain is generally non-zero because of equation (103). This implies that To put it another way, the kinetic energy differences between refinement levels have to be compensated by numerical cooling or heating in order to maintain energy conservation. This can be alleviated by incrementing the SGS turbulence energy by ∆(ρK) when the cutoff scale is shifted from ∆ l+1 to ∆ l = r∆ l+1 , where r > 1 is the refinement ratio: We then have conservation of the total (resolved plus unresolved) kinetic energies, and e crs = ρe. As demonstrated in [118], equation (107) works very well for fully developed homogeneous turbulence, but leads to erroneous projections of the SGS turbulence energy to coarser grids if the flow structure is strongly inhomogeneous. A particular problem is posed by non-turbulent bulk flows such as gas accretion into gravitational wells. In this case, applying the increment ∆(ρK) defined by equation (105) can greatly overestimate the energy difference associated with turbulent velocity fluctuations. A tentative solution to this problem follows along similar lines as the rescaling procedure used in [74]. By extrapolating the turbulent energy on the length scale ∆ l+1 via a power law to ∆ l , we have Substitution into equation (106) As a result, we have and the total energy is conserved because (ρE) crs + (ρK) crs = ρE + ρK .  (105), For grid refinement, the data from the parent grid are interpolated to fine-grid values ρ * n , (ρU n ) * , etc. such that One can set ρ n = ρ * n , (ρU ) n = (ρU ) * n , and (ρE) n = (ρE) * n , but then the kinetic energy difference resulting from the conservative interpolation of the momenta has to be compensated. By defining the energy corrections the interpolated values of the SGS turbulence and internal energies can be adjusted as follows: The conservation of total energy follows from where the last equality follows by substituting equations (114) and (115) for ∆(ρe) and ∆(ρK), respectively. Equation (116) is formally the same as the correction rule for the SGS turbulence energy in [74]. In contrast to the method outlined above, however, momentum conservation is not fulfilled by [74] because the full kinetic energy difference between the finer and coarser levels is compensated by internal energy and then the velocities are re-scaled to compensate the power-law correction of the SGS turbulence energy (see also Section 6.3).

Shear-improved model
For inhomogeneous and non-stationary flows, dynamical procedures can be applied to adjust SGS model coefficients to local flow conditions (see Section 4.2). A completely different idea was put forward by Lévêque to calculate the turbulent stresses in LES of wall-bounded turbulence [69]. Rather than adjusting the eddy-viscosity coefficient C ν , the numerically resolved velocity field is decomposed into a mean flow u and turbulent fluctuations u . It is then argued that the turbulence energy flux for the Smagorinsky model should linearly depend on the shear associated with the fluctuating component, i. e., where | S | is the rate of strain of the mean flow. If the flow is laminar, u u and turbulence production vanishes because |S|−| S | 0. For developed isotropic turbulence, on the other hand, u u. In this case, the standard Smagorinsky model with Σ ∝ |S| 3 applies. For intermediate cases, the model corrects the energy flux Σ by taking into account interactions of the grid-scale fluctuations with the mean shear. This is why it is called the shear-improved model.
As proposed in [118], this idea can be carried over to the SGS turbulence energy model (see Section 3.3) by defining the shear-improved eddy-viscosity closure as where S * ij is the trace-free part of the rate-of-strain tensor associated with the fluctuating component of the flow: The turbulence energy flux is then given by Apart from the divergence term, the expression on the right-hand side is motivated by the generalized Kármán-Howarth equation in [69]. A difficulty of implementing the shear-improved model is the computation of the mean flow u , which is an ensemble average over the flow velocity u. A practical solution is to find an approximation to the mean flow by smoothing u with a temporal low-pass filter. For statistically stationary turbulence, it is possible to use an exponentially weighted recursive average. In component notation, an estimate for the mean velocity at time t n+1 is calculated as weighted sum of the estimate at the previous time step and the local velocity u (n+1) i for each grid cell (cell indices are omitted): The weighing coefficients are defined by where T c is the constant integral time scale of the flow [17]. Changes occurring on time scales smaller than the smoothing scale T c are suppressed in [u i ]. Compared to dynamical procedures, this algorithm is very easy to implement and computationally much cheaper. In [69], it is demonstrated that the shear-improved Smagorinsky model with exponential smoothing performs well in LES of plane-channel flow and reproduces data from direct numerical simulations. However, the simple exponential smoothing algorithm produces an estimate [u i ] that lags behind the ensemble average u i if the mean flow evolves. To address this problem, the so-called Kalman filtering technique is introduced in [17]. The Kalman filter adapts itself to an unsteady mean flow by dynamically adjusting the weights of the recursive relation (121), depending on the variances of the mean flow evolution, [u i ] (n) −[u i ] (n−1) , and the detected deviation from the mean, u , which is defined by the ratio between the error variance of the smoothed component and the total error variance, including the fluctuating component. Since the error variances have to be evaluated at time t (n+1) , a predictor-corrector scheme is used: 1. Given the error variance P (n) i at time t (n) , the prediction for t (n+1) is where σ (n) Here it is assumed that the typical correction of the mean flow,

The Kalman gain is then given by
where σ 2 (n) is the contribution of the fluctuating component δu 3. The corrected error variance for the next step is given by In a statistically stationary state, the velocity fluctuations should be of the order σ (n) δui u c . In this case the Kalman filter corresponds to simple exponential smoothing: The two filter parameters, T c and u c , have to be chosen such that u c is roughly the integral velocity of turbulence if the flow enters a steady state and T c is the characteristic time scale over which the flow evolves. In [18], LES of turbulence produced by the flow past a cylinder were shown to agree well with experimental data if Kalman filtering is applied with T c and u c being set to the inverse of the expected vortex-shedding frequency and upstream velocity, respectively. A similar computational problem in astrophysics is an isothermal spherical cloud, which is bound by a static gravitational potential, in a homogeneous wind. The initial density profile of the cloud is determined by hydrostatic equilibrium. As mass is stripped from the cloud by the wind, a turbulent wake forms in the downstream direction. This problem was originally investigated as a simple model for the infall of small subcluster into the ICM of a much more massive cluster, computed in the frame of reference attached to the center of mass of the subcluster [51]. LES with the shear-improved SGS model are presented in [118], where the Kalman filter parameters are given by the velocity of the wind and the turn-over time of the largest eddies. Figure 8 compares the gas density and flow structure for two runs, one with the standard SGS model and the other one with Figure 8: Slices of the mass density (left), squared vorticity (middle), and specific SGS turbulence energy (right) for AMR simulations of a gravitationally bound gas cloud in a wind after 3 Gyr of evolution [118]. The box size is 4 Mpc.
the shear-improved model. In both cases, a 128 3 root grid and three levels of refinement are used. The AMR control variables are the squared vorticity and the compression rate [51,112]. The latter is defined by the substantial time derivative of the velocity divergence d and tracks down the bow shock in front of the cloud. The top panels in Figure 8 show that the strain associated with the bow shock produces substantial SGS turbulence energy if the standard SGS model is applied. One can also discern the down-scaling of the SGS turbulence energy that is transported with the wind, as it enters the refined regions around the cloud. This is based on the algorithm explained in the previous section. Kelvin-Helmholtz instabilities between the low-density wind and the high-density gas in the cloud cause vortex shedding, which in turn produces turbulence. This becomes manifest in large values of ω and K in the turbulent wake that extends from the cloud towards the right. In contrast to the turbulent wake, however, the shear experienced by the gas when it is passing the bow shock is not associated with turbulence. As a result, the steep increase of K is largely a spurious effect induced by the eddy-viscosity closure (49). Indeed, the SGS turbulence energy in the shocked wind is substantially reduced if the shear-improved closure (118) is applied (bottom right panel in Figure 8). In comparison to the standard SGS model, production is also suppressed at the front side of the cloud and, hence, vortex shedding is the dominant source of turbulence production. This has a clearly visible impact on the structure of the turbulent wake.
In general, an appropriate prior choice of the filter parameters might not be obvious. Nevertheless, they can be calibrated a posteriori by performing low-resolution test runs. This is shown for cosmological simulations in [118]. In this case, the Kalman filter also has the merit of separating the turbulent flow in clusters from the gravity-driven bulk flows (see Section 6.3.3).

Thermonuclear combustion in white dwarfs
Among the various possible scenarios for supernovae of type Ia are thermonuclear explosions of gas accreting white dwarfs in close binary systems [50,98,49]. If the mass of a white dwarf approaches the Chandrasekhar limit, explosive carbon and oxygen burning is ignited [81]. Owing to the degeneracy of white dwarf matter, the thermonuclear reaction zones propagate as thin flame fronts, whose thickness δ f and propagation speed s f are determined by the very high thermal conductivity of the fuel. This mode of burning is called deflagration. Since the burned material has lower density than the fuel, it rises because of its buoyancy. Consequently, the energy released by thermonuclear deflagration drives convection. Since eddies exert strong shear on rising bubbles of burning material, they are deformed into mushroom-like shapes and Kelvin-Helmholtz instabilities at the surface are rapidly producing turbulence [75]. Eventually, this results in a very complex flame front with a fractal structure that cannot be resolved in numerical simulations over the full range of length scales.
This problem was addressed by performing LES with an effective propagation speed of the turbulent flame front ( [95,115,101], to mention just a few examples). In these simulations, the flame front is tracked by means of the level set method [82,96], which is able to follow complex topological changes by determining the interface between fuel and burned material as the spatial surface for which a signed distance function is zero. The evolution of the distance function from given initial conditions is given by the advection velocity relative to the fluid plus the flame propagation speed. In a fully resolved simulation, the flame propagation speed would be given by the microscopic flame speed s f (also called laminar flame speed). If the flame front is underresolved, however, the wrinkling and folding of the front by turbulent eddies below the grid scale leads to an enhanced rate of energy release. In this case, the relevant time scale is not the microscopic diffusion time scale but the turn-over time of numerically unresolved eddies. Simple dimensional reasoning implies that s f has to be replaced by the turbulent flame speed s t = √ 2K [79,80,88], where K is the SGS turbulence energy given by equation (64). In [117], the formula is proposed for a smooth transition between laminar and turbulent flame propagation. A further complication comes from the pronounced anisotropy of turbulence at the flame front [23,75]. While the burned material inside the flame is highly turbulent, there is little or no turbulence in the fuel just outside the flame. In a way, this is similar to walls in terrestrial flows. For this reason, it is important to apply the dynamic produced explained in Section 4.2 to locally calculate the eddy-viscosity coefficient, which determines the rate of production of K. The minimal length scale for which the flame front is affected by turbulence is given by the Gibson scale G . If v ( ) is the mean turbulent velocity fluctuation on the length scale , then G is implicitly given by the condition i. e., the turbulent velocity fluctuation on the Gibson scale equals the microscopic flame speed. For G , the turn-over time associated with v ( G ) is much longer than the crossing time of the flame front through an eddy of size . On these scales, the flame front is virtually unaffected by turbulence. If G , on the other hand, eddies will significantly distort the flame front on length scales between G and . In LES of a thermonuclear supernovae, the Gibson scale G is much smaller than computationally feasible grid resolutions ∆ during most of the deflagration phase and, consequently, a turbulent flame speed model has to be applied. 17 As an illustration of the level set method with the turbulent flame speed (125), Figure 9 shows snapshots of the flame front for a thermonuclear supernova simulation from [115]. The turbulent velocity fluctuations on subgrid scales, √ 2K, are shown as color shades at the flame surface. The asymptotic value of turbulent flame speed is s t C t √ 2K if s t s f . In this simulation, a Poisson process is used to randomly place small ignition spots in the high-density core of the white dwarf. The statistics of this process is based on a simple model for temperature fluctuations produced by convection in pre-ignition phase. Although the number of ignition resulting from this model appears to be by far too large in the light of recent numerical studies of the ignition process [81], the simulation nevertheless demonstrates in an exemplary manner how the turbulent deflagration progresses. At early time (a), one can see a large number of small bubbles generated by the stochastic ignition process at distances of the order 100 km from the center of the white dwarf. As the burning bubbles are rising from the centre, they begin to form the typical Rayleigh-Taylor mushroom shapes (b). At this point, the effective flame propagation speed is already dominated by turbulence. After about half a second (c), the turbulent flames mostly have merged into a single structure of about 1000 km diameter. Then the burning front rapidly expands to much larger radii and causes the white dwarf to explode after roughly one second (d).
Since the the asymptotic rising velocity of a perturbation of size l due to its buoyancy is given by the Sharp-Wheeler relation [120] v RT ( ) = 0.5 g eff , where g eff is the effective gravity associated with the density contrast at the interface between burned and unburned material, s t v RT (∆) was used as a simple turbulent flame speed model (e. g., [38]). In [117], buoyancy is included as an ad-hoc source term in the SGS turbulence energy equation (27): Here, Γ is defined such that it vanishes everywhere except for the close vicinity of the flame front. Moreover, Γ is assumed to be non-zero only if the so-called fire polishing length λ fp = 4πs 2 f /g eff is smaller than the grid resolution ∆. The fire polishing length is the smallest length scale on which perturbations in the flame front are Rayleigh-Taylor unstable. If the buoyancy term dominates over the turbulent cascade, √ 2K ∼ v RT (∆) is obtained as asymptotic solution. However, different numerical studies indicate Kolmogorov scaling on small scales [133,101,23], corresponding to Γ Σ in equation (27). Although the impact of the turbulent flame speed model on the energy release was demonstrated in various simulations [79,117], it tends to become less significant if the burning is resolved down to very small scales. This has become possible with the application of of AMR and the power of contemporary computing facilities [72].
However, since the predictions from pure deflagration models are not consistent with observational properties of the majority of type Ia supernovae, a transition from the deflagration phase to a supersonic detonation is now considered as the most likely explosion mechanism of Chandrasekharmass white dwarfs [98,49]. A theoretical explanation of such a transition is a long standing problem. One possibility is that very strong local velocity fluctuations at the onset of distributed burning might trigger the transition to a detonation [80,56,70,130]. The statistical distribution of turbulent velocity fluctuations following from the SGS turbulence energy K in a high-resolution simulation (see Figure 10) shows that such events might indeed occur [97]. The very wide tail indicates the strong intermittency of turbulence at the flame front. On the basis of this result, a criterion was recently proposed to set off detonations in LES of thermonuclear supernovae [24]. . A similar simulation is presented in [115]. The spatial scale is adjusted to the bulk expansion of the white dwarf by means of a co-moving grid technique [99].  Figure 10: Flame front in a high-resolution simulation of a thermonuclear supernova [97]. High turbulent velocity fluctuations occur in the reddish regions. Extremely strong fluctuations that could trigger a transition to a detonation are indicated by the green arrows. By courtesy of Friedrich Röpke.

Galaxy simulations
Numerical simulations of galaxies, particularly from cosmological initial conditions, cannot fully resolve processes such as star formation and feedback from supernova (SN) explosions. As suggested in [55,106,12], a production term Σ due to feedback can be incorporated as an additional source term in the SGS turbulence energy equation (27). The simplest model for this term is where e SN is the typical energy released by a core-collapse supernova and τ ff ∝ ρ −1/2 is the free-fall time scale. The parameter C controls the effective time scale of the feedback. The above expression for Σ follows from the simple Kennicutt-Schmidt relationρ ∝ ρ −3/2 for the star formation rate, where the factor C ρ/τ ff is some fraction ofρ . Apart from stellar feedback, turbulence is produced by instabilities in the ISM, such as gravitational, thermal, and hydrodynamical instabilities. Typically, energy is injected on numerically resolved length scales by these instabilities and transferred to smaller scales by the turbulent cascade. If we simply assume that a turbulent velocity of magnitude V is produced on the length scale L, the energy flux through the turbulent cascade is of the order Σ ∼ ρV 3 /L [11]. In LES, Σ is defined in terms of the grid scale and the shear of the numerically resolved flow via the closure (55). As a very crude model, let us consider the local equilibrium between production and dissipation, i. e., Σ + Σ ∼ ρ . By substituting equations (55), (129), and ∼ ρK 3/2 /∆, the equilibrium condition reads This relation implies the turbulent pressure where is the dynamical time scale associated with energy transfer through the cascade. Depending on the energy ratio K/e SN and the ratio of the dynamical and feedback time scales, τ /τ ff , the production of SGS turbulence energy is dominated by the turbulent cascade or by supernovae. The star formation rate, which in turn determines the feedback rate, can be parameterized by the gas density, temperature, and turbulence intensity [62,84,47,32,83]. The two main parameters appearing in these parameterizations are the turbulent Mach number M and the virial parameter α [7], which correspond to the ratios of the turbulence energy to the internal and gravitational energies: Here, σ( ) is the turbulent velocity dispersion on the length scale , c s the speed of sound, and G Newton's constant. A major problem in galaxy simulations is the choice of a characteristic scale of star formation. In [11] it is argued that is given by the Jeans length for gravitational instability, i. e., = λ J := c s π γGρ where γ is the adiabatic exponent of the gas. Since the SGS model predicts the specific turbulence energy K = 3σ 2 (∆)/2 on the grid scale ∆, it is possible to estimate the turbulent velocity dispersion in gravitationally unstable, star-forming clouds as where the scaling exponent η is in the range between 1/3 and 1/2, depending on the compressibility of the gas and the intermittency of turbulence [60,113,112,48]. The star formation rate is then given byρ The mass fraction SFR(α , M ) that is converted into stellar mass per free-fall time is modeldependent. Basically, this factor accounts for the gravo-turbulent fragmentation of star-forming clouds. 18 For example, a power-law function based on data from numerical simulations is proposed in [62]. In [84] it is assumed that SFR(α , M ) is determined by the log-normal probability density function, whose width depends on the turbulent Mach number M [85,33], and a critical density that is controlled by the virial parameter α . The different parameterizations are compared and calibrated by numerical data in [32]. A simple star formation law, for which SFR(α ) is solely determined by the virial parameter, is proposed in [83].
A star formation rate proportional to the density of molecular hydrogen, f H2 ρ instead of the total gas density ρ, was suggested on grounds of the observed tight correlation between the star formation rate and the column density of molecular hydrogen [64,43]. Since molecular hydrogen forms only in the cold phase of the interstellar medium, the star formation law (135) is further modified by replacing ρ and c s with the mean mass density and speed of sound, respectively, in the cold phase [11]. Since the separation into cold and warm phases cannot be fully resolved in galaxy simulations, a model such as [123] has to be employed to estimate the fractional densities of the two phases. This entails additional turbulence energy production by cooling instabilities on length scales below the grid resolution [11,54]. If the resulting internal driving due to feedback and cooling instabilities, Σ int = Σ + Σ TI , along with the closure (55) for the compressible turbulence energy flux Σ are incorporated into equation (64), a full model for the numerically unresolved turbulence energy budget in galaxy simulations is obtained [12].
In combination with the two-phase model for the gas and a particular flavor of the star formation and feedback models outlined above, this SGS turbulence energy model has recently been applied to adaptively refined LES of isolated disk galaxies. As schematically shown in Figure 11, the computation of the SGS turbulence energy ρK plays a central role. In Figure 12 (left), an equatorial slice of the gas density illustrates the fragmentation of the inner gaseous disk into dense starforming clumps (bluish regions) in one of these simulations. The dilute gas with densities below 1 cm −3 (reddish regions) is partially heated by supernovae. The fragmentation of the disk produces turbulence, as indicated by the slice of ρK in Figure 12 (right). Particularly large turbulent energies are found in regions with strong feedback, which suggest that the non-thermal feedback from supernovae plays a significant role in the production of SGS turbulence energy. This is indeed confirmed by the statistics of the production terms plotted in Figure 13 (left). This plot shows profiles of Σ int Σ tot = Σ + Σ TI Σ + Σ + Σ TI and Σ Σ tot = Σ Σ + Σ + Σ TI , 18 The assumption = λ J does not imply that a cloud of size collapses and is turned into stars as a whole. The internal structure of star-forming clouds with strong density fluctuations due to supersonic turbulence is usually not resolved in galaxy simulations. For this reason, the Jeans length for the gas density in a grid cell serves only as a reference length scale for the mean properties of the numerically unresolved star-forming clouds. If the clouds are partially resolved, however, this assumption needs revision. Figure 11: The SGS turbulence energy ρK has several functions in the multiphase model for the turbulent ISM in disk galaxies [12]. It determines the shape and width of the mass density PDF of cold clumps and the effective pressure equilibrium between the cold and warm phases. Via the density PDF, it influences the star formation rate, which acts back on ρK through feedback. Other sources of ρK are the turbulent cascade and cooling instabilities.
which are calculated by averaging Σ, Σ int , and Σ tot over bins of SGS turbulence energy per unit mass. The maximum of Σ/Σ tot around K ∼ 100 (km/s) 2 implies that the turbulent cascade maintains a ground level corresponding to a turbulent velocity dispersion around 10 km/s, while the internal driving caused by supernovae excites much stronger turbulent velocity fluctuations. Remarkably, the average turbulence energy is quite close to the equilibrium value implied by the balance between production and dissipation. This can be seen in Figure 13 (right), where the ratio of the equilibrium pressure to the dynamical turbulent pressure, P eq,tot P sgs is plotted against P sgs = 2 3 ρK. Also shown is the asymptotic equilibrium pressure P eq,int P sgs if Σ tot Σ int Σ, i. e., internal driving dominates. This is the case for high turbulence intensity. Since P eq,int /P sgs is only slightly above unity, turbulence is close to equilibrium in this regime. For lower SGS turbulence energies, equation (130) tends to overestimate the turbulent pressure, which indicates that turbulence production exceeds dissipation. In this case, the turbulent cascade contributes significantly to the production. For comparison, also profiles of the thermal pressures of the cold and warm phases are plotted. On the one hand, P eq,tot is small compared to the warm-phase pressure, except for the highly turbulent regions with strong feedback, where P w ∼ P eq,tot ∼ P sgs . On the other hand, P eq,tot P c . The pressure equilibrium between the cold and warm phases, which is one of the basic assumptions of the model, is therefore dominated by the turbulent pressure in the cold clumps and the thermal pressure in the warm medium. Figure 12: Projections of the number density (left) and SGS turbulence energy (right) on the equatorial plane in an isolated disk galaxy simulation [12]. Figure 13: Left: profiles of the rate of production through the turbulent cascade, Σ, and internal driving due to numerically unresolved feedback from supernovae and thermal instabilities, Σ int , relative to Σ tot = Σ + Σ int for the disk galaxy simulation shown in Figure 12 [12]. Right: profiles of the pressure ratios P X /P sgs , where P sgs = 2 3 ρK is the numerically unresolved turbulent pressure, for the thermal pressures P c and P w of the cold and warm phases and the equilibrium pressures P eq,tot and P eq,int defined by equations (136) and (137), respectively. By courtesy of Harald Braun.

Cosmological simulations
Since cosmological scale structure formation produces a strongly clumped medium through gravitational contraction or collapse, the numerical simulation of turbulence on cosmological scales is particularly challenging. Potential production mechanisms of turbulent flows in the baryonic gas are mergers between dark-matter halos and the accretion of gas into halos, but also feedback form active galactic nuclei and winds produced by strongly star-forming galaxies [28,59]. A largely open question concerns whether the gaseous component of halos is in a state of developed turbulence. While there is no doubt about turbulence in the interior of galaxies, there is no direct observational evidence yet for turbulence on larger scales, particularly in the intracluster medium (ICM). Theoretical and numerical studies suggest that magnetic fields play a key role in the physical dissipation mechanism and, possibly, the onset of instabilities, but the associated length scales are highly uncertain [35,13]. Notwithstanding these uncertainties, turbulent flows resulting from cosmological structure formation are numerically investigated by computing the evolution of an ideal fluid subject to the gravitational potential of dark matter, which is modeled as a collisionless N -body system [8]. To follow gravitational collapse, AMR is an essential if Eulerian grid codes are used. As pointed out in Chapter 5, turbulence is generally underresolved in a clumpy medium, even at very high refinement levels.
This problem was addressed for the first time by combining LES and AMR in [74]. 19 The main idea is to solve the filtered Euler equations (40)(41)(42) and the SGS turbulence energy equation (44) in co-moving coordinates. For grid refinement and de-refinement, [74] apply the power-law relation (108) between the SGS turbulence energies at different refinement levels. In contrast to the algorithms outlined in Section 5.1, however, their implementation is not conservative by construction. The SGS turbulence energy decrement −∆(ρK) in the case of refinement from a coarser level, for example, is simply compensated by adjusting the interpolated momenta (ρU ) * n such that which follows by setting (ρK) n = r −2η (ρK) * n and (ρU ) n = (ρU ) * n 1 + While the first relation is identical to equation (116), the second relation, when substituted into the total energy balance, implies This is just the standard method of compensating the difference of the kinetic energies between refinement levels entirely by internal energy, which is applied in addition to the kinetic energy transfer (138). Analogous relations are applied to average the data from a finer to a coarser level. The rationale of equation (138) is that the resolved turbulent velocity fluctuations should increase from a coarser to a finer level. However, as argued in Section 5.1 and in [118], conservative interpolation of the momenta to a finer level entails already an increase of the resolved kinetic energy that encompasses and, in some cases, even overestimates the scale-dependence of the numerically resolved turbulence energy. Both energy and momentum conservation is guaranteed by compensating this increase with ∆(ρK) and ∆(ρe) defined by equations (114) and (115), respectively. Despite its shortcomings, the method based on [74] has its merits as a first approximation. For example, Figure 14 shows the density of the baryonic gas and the small-scale turbulence in an adiabatic cosmological simulation of a cluster from [74]. As indicator of turbulence, the magnitude of the numerically unresolved turbulent velocity fluctuation √ K is scaled down from refinement level l to the minimal cell size at the maximal level l max via the power law factor r (l−lmax)η . It is particularly interesting that the infall of a subhalo, which is marked by the small box in slice (a), produces a pronounced turbulent wake, quite similar to the idealized scenario discussed in Section 5.2. However, a shear-improved model was not used in this simulation. Even at z = 0, the turbulent velocity slice (d) in Figure 14 shows a clear trace of the minor merger, which is not discernible in the density slice (c).

Turbulence production and support against gravity
Adiabatic simulations of cosmological structure formation with Enzo [15] show that the SGS turbulence energy traces different production mechanisms of turbulence in the intracluster medium (ICM) and the warm-hot intergalactic medium (WHIM) [52]. While the former is found at high densities and temperatures and is dominated by weakly compressible, subsonic turbulence, the WHIM is constituted by gas of lower density, which is undergoing compression by shocks. The plot in Figure 15 shows different histories of the mean thermal and SGS turbulence energies in the ICM and WHIM. While the SGS turbulence energy in the ICM peaks around redshift z = 1, which can be interpreted as a consequence of mergers between galaxy clusters, the energy in the WHIM gradually rises as a result of the continuous turbulence production by accretion shocks. These are caused by the infall of low-density gas into the potential wells of clusters and filaments, which accelerates the gas to supersonic speed.  Figure 15: Evolution of the mean internal (solid) and SGS turbulence (dashed) energies in the ICM and the warm-hot intergalactic medium (WHIM) of a galaxy cluster [52]. By courtesy of Luigi Iapichino.
In addition to other indicators, [52] analyses the relative importance of the turbulent and thermal pressures for the support of the gas against gravity. As shown in [110], the local support of self-gravitating gas due to the thermal pressure P is given by The above expression is derived by taking the divergence of the equation for the gas velocity, which results in an equation for the rate of compression of a fluid parcel, −Dd/Dt, where D/Dt signifies the substantial time derivative and d the divergence of the flow. In the weakly compressible limit, the thermal support reduces to Λ therm −(∇ 2 P )/ρ. By considering equation (23), it immediately follows that the influence of turbulence below the grid scale can be expressed as where P sys = 2 3 ρK is the turbulent pressure on the grid scale. However, Λ sgs does not account for the effect of the non-diagonal stresses τ * ij . The contribution of numerically resolved turbulence, on the other hand, is given by where ω = ∇ × u is the vorticity and |S| the rate of strain (see equation 15). The first term on the right-hand side is associated with the pressure-like support caused by turbulent eddies, the second term with the compression of the gas by shocks. Figure 15 shows a mass-weighted histogram of the ratio for different overdensities of the baryonic gas. The ratio r tp is an approximation to the ratio of the turbulent and thermal support functions, (Λ turb + Λ sgs )/Λ therm , which was used for the sake of comparability with [132]. The distribution of r tp plotted in Figure 16 shows two peaks, one at intermediate densities and the other at higher densities, which can be associated with the WHIM and the ICM. There is a trend of decreasing r tp toward high densities. This reflects the smaller Mach numbers of turbulence in the ICM, but a caveat is that only positive contributions to the support are taken into account. Since it is demonstrated in [110,67] that Λ turb can be predominantly negative in the presence of shocks, the question of the turbulent support of the gas in clusters needs to be revisited. Figure 16: Mass-weighted histogram of the turbulent-to-thermal support ratio vs. baryon overdensity (the solid line indicates the mean value for given overdensity) [52]. By courtesy of Luigi Iapichino.

Gravitational collapse of gas in primordial halos
The Enzo implementation of the SGS model for cosmological fluid dynamics was also applied to the gravitational collapse of gas clouds in primordial atomic cooling halos. This scenario is potentially relevant for the formation of intermediate-mass black holes through direct collapse and subsequent accretion [66,65,68]. The resulting black holes are possibly the seeds for supermassive black holes at the centers of galaxies. Direct collapse into a supermassive star, which subsequently collapses into an intermediate-mass black hole, can occur if the molecular hydrogen formation is suppressed through photodissociation by a Lyman-Werner radiation background. To follow the collapse of the primordial gas down to AU scales, deep-zoom in simulations with many level of refinements have to be performed. In [65], both LES and ILES of several different halos of masses of the order 10 7 M are presented. One of the most remarkable results of this study is that the additional turbulent viscosity produced by the SGS model favors disk-like structures around collapsed objects in LES. This can be seen from the comparison of density slices in the innermost regions around the density peaks. Figure 17 shows that more or less compact collapsed structures are produced in ILES. In contrast, more extended and disk-like structures are found in the corresponding LES, as shown in Figure 18. Apart form morphological differences, the SGS model also has a significant influence on the accretion of mass by the protostars that are formed through the collapse of gas clouds in the halos. As mentioned in Section 3.4, it is possible to follow the evolution of gravitationally bound dense objects, such as protostars, by inserting sink particles [31,127]. This method is applied in [68] to follow the accretion history of protostars in atomic cooling halos. The resulting time evolution of the accretion rate is plotted in Figure 19 for simulations of three different halos, using both ILES and LES. As one can see, the accretion rates reach peak values around 10 M yr −1 roughly within 10 4 yr. This value agrees with the theoretical expectation for Bondi-Hoyle accretion. A comparison of ILES and LES suggests a systematically higher accretion rate for LES. This trend is confirmed by calculating the cumulative masses of the sink particles (see right plot in Figure 19), which reach masses above 10 5 M . As a result, the SGS model favors the formation of higher black hole masses.

Turbulent velocity dispersion
A key question that arises in connection with analyzing turbulence production by cosmological structure formation concerns the turbulent velocity dispersion σ turb . In [53], σ turb is defined as σ turb = √ 2K, which implies that σ turb is identified with the magnitude of the turbulent velocity fluctuations on the grid scale of the simulation. Although this variable obviously depends on numerical resolution, it is nevertheless useful to infer the dependence on various factors that influence the production of turbulence. For example, Figure 20 (left plot) shows the ratio of the turbulent and thermal pressures, where P t = P sgs = 2 3 ρK, and the Doppler broadening parameter b t = σ turb / √ 3 = 2K/3 [29] for a cosmological simulation with radiative background and cooling in a box of 10 Mpc h −1 comoving size. The analysis is carried out for data cubes at redshift z = 2.0. For the intergalactic medium (IGM), which is usually defined by moderate baryonic overdensities ρ/ρ 0 in the range from 1 to about 50, b t is found to increase with density from roughly 1 to 10 km/s. The WHIM, on the other hand, has a flat turbulent velocity dispersion, corresponding to P t /P ∼ 0.1. The phase plot in Figure 20 shows that the ratio b t /b varies over several orders of magnitude for different densities and temperatures and reaches peak values ∼ 1. The WHIM is associated with gas that is heated by accretion shocks to temperatures between 10 5 and 10 7 K. For this reason, the WHIM tends to be more turbulent then the diffuse gas in the IGM. The phase diagram also shows that the distinction between gas in the IGM (1 ≤ ρ/ρ 0 ≤ 50) and the WHIM (10 5 K ≤ T ≤ 10 7 K) is not mutually exclusive and somewhat arbitrary.
A resolution-dependent turbulent velocity dispersion is avoided by the K-L hybrid model for ] Figure 18: The same halos as in Figure 17 computed with an explicit SGS model [65]. By courtesy of Muhammad Latif.  Rayleigh-Taylor-driven turbulence (see Section 3.4), however, at the cost of smearing out turbulent structures over most of the numerically resolved dynamical range. This model was used to simulate the production of turbulence by AGN feedback in galaxy clusters [105,14]. Another caveat of both the K-L model and the SGS model based on [74] are the constant closure coefficients. Strictly speaking, constant-coefficient models are applicable only to statistically homogeneous and stationary turbulence. Turbulence produced by cosmological structure formation and AGNs, however, is highly inhomogeneous and non-stationary.
To ameliorate this problem as well as the dependence of σ turb on the grid scale, the shearimproved SGS model outlined in Section 5.2 has recently been applied to cosmological structure formation [118]. By running simulations of the Santa Barbara cluster [46] with the Nyx code [1], it is demonstrated that the application of the standard SGS model for homogeneous turbulence can produce biases, particularly in regions of active turbulence production, such as the WHIM. In Figure 21, the growth of the cluster is illustrated by slices at different redshifts from a simulation performed with the shear-improved SGS model. There clearly is a correlation between the high vorticity and SGS turbulence energy in the cluster and a sharp drop at the outer accretions shocks. Although a meaningful comparison to the results in [53] cannot be made at this point because the Santa Barbara cluster is based on a matter-dominated universe without cosmological constant and the gas dynamics is adiabatic, some general conclusions regarding the nature of turbulence in clusters can be drawn. In [118], the turbulent velocity dispersion is defined by where U is the fluctuating component of the velocity computed with the Kalman filter (see Section 5.2). The average values of σ turb for logarithmic bins of the baryonic overdensity as well as the corresponding radial profiles are plotted for different numerical resolutions in Figure 22. Except for the lowest-resolution case, the radial profiles of σ turb show little sensitivity to numerical resolution. This is an important property of the turbulent velocity dispersion defined by equation (144). While the profiles are nearly flat for the ICM, there is a sharp drop around 10 Mpc, which is about the radius of the outer accretion shocks in Figure 21. Since the gravitational pull of the cluster causes a large non-turbulent bulk flow toward the center, the total flow velocity beyond the accretion shocks is much higher than σ turb . As a function of the overdensity ρ/ρ 0 , σ turb gradually increases towards the center of the cluster. Although the simulations in [53] differ in important aspects, a roughly similar trend can be seen for the Doppler broadening parameter b t in Figure 20. The turbulent kinetic energy 1 2 ρσ 2 turb and the energy flux Σ through the turbulent cascade (see Section 2.2) define the dynamical time scale The radial profiles of 1/τ plotted in Figure 23 show that the dynamical time scale in the ICM is several Gyr. The pronounced peak at 10 Mpc radius is a further indication of turbulence production by the accretions shocks. Analogous to equation (145), the dissipation time scale can be defined as For statistically stationary turbulence, the balance between turbulence production and dissipation, Σ ∼ ρ , implies τ ∼ τ . This can indeed be seen in Figure 23 for the central region of the cluster.
In the vicinity of the accretion shocks and outside the cluster, however, τ τ . In this case, the flow is far from equilibrium. These results demonstrate how (nearly) scale-invariant quantities computed with the shear-improved SGS model can be utilized to investigate statistical properties of turbulence. Figure 21: Slices of the baryonic mass density (left), vorticity modulus (middle), and specific SGS turbulence energy (right) at different redshifts for a simulation of the Santa Barbara cluster [118]. The box size is 64 Mpc h −1 .  (144) vs. baryonic overdensity (left) and radius from the center (right) for simulations of the Santa Barbara cluster with different numerical resolutions [118]. Starting from a uniform root-grid with N 0 grid cells, refinement by overdensity and vorticity modulus up to l max levels is applied. N 0 =256 3 , l max =2 N 0 =256 3 , l max =1 N 0 =128 3 , l max =2 N 0 =128 3 , l max =1 Figure 23: Radial profiles of the inverse dynamical (left) and dissipation (right) time scales for the same simulations as in Figure 22 [118].