A finite volume dynamic large-eddy simulation method for buoyancy driven turbulent geophysical flows

Ocean Modelling 17 (2007) 199–218 www.elsevier.com/locate/ocemod

A finite volume dynamic large-eddy simulation method for buoyancy driven turbulent geophysical flows Filippo M. Denaro a,*, Giuliano De Stefano a, Daniele Iudicone b, Vincenzo Botte b a

Dipartimento di Ingegneria Aerospaziale e Meccanica, Seconda Universita` di Napoli, 81031 Aversa, Italy b Stazione Zoologica ‘‘Anton Dohrn”, 80125 Napoli, Italy Received 15 March 2006; received in revised form 9 February 2007; accepted 9 February 2007 Available online 1 March 2007

Abstract A large-eddy simulation method based upon the finite volume approach is developed and evaluated for buoyancy driven turbulence in the absence of rotational effects. The effect of the sub-grid scale motions is modelled by using the dynamic scaling formulation, without introducing any ad hoc assumption for the vertical stratification. This way, the eddy coefficients for momentum and temperature equations are independently computed, while avoiding to introduce unjustified buoyancy production terms in the constitutive equations. As a result, the computational cost is considerably reduced with respect to the classical stratification formulation. The method is presented in detail, by stressing the particular features of the finite volume large-eddy simulation approach. The resulting numerical code is validated against a direct numerical simulation. Numerical experiments are conducted by simulating the thermal buoyancy driven turbulence near the water surface generated in a finite-depth stably stratified horizontal layer with an isothermal bottom surface. Diagnostics include time evolution of kinetic and thermal energy as well as energy spectral distribution. Interesting results are obtained by making a comparison with a reference spectral solution. The method is further validated for a turbulent ocean test case, that is the deepening of the mixed layer in a stable stratified fluid. Ó 2007 Elsevier Ltd. All rights reserved. Keywords: Buoyancy driven turbulence; Large-eddy simulation; Dynamic subgrid-scale model

1. Introduction Buoyancy driven turbulence is of fundamental importance in a large-variety of physical phenomena. The strong mixing associated to turbulence can have a profound impact on the surrounding environment. This is to a degree the case for the oceanic mixed layer (OML), a region of the upper open-ocean that is almost vertically uniform due to vigorous turbulent mixing processes (e.g., de Boyer Monte´gut et al., 2004). This layer

*

Corresponding author. Tel.: +39 0815010296. E-mail addresses: [email protected] (F.M. Denaro), [email protected] (G. De Stefano), [email protected] (D. Iudicone), [email protected] (V. Botte). 1463-5003/$ - see front matter Ó 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.ocemod.2007.02.002

200

F.M. Denaro et al. / Ocean Modelling 17 (2007) 199–218

is responsible of the transfer of mass, momentum, and energy promoting almost all oceanic motions and it is the region that directly interacts with the atmosphere. Turbulence in the OML is produced by surface fluxes of heat and momentum as well as by surface waves. The shear due to wind-driven currents can efficiently mix the upper ocean via, e.g., Kelvin–Helmholtz instabilities. The interaction between the wind-driven shear currents and the mean residual Stokes drift of surface waves produces the Langmuir circulation which in most cases is highly turbulent (e.g., McWilliams et al., 1997). Thermal convection is also common in the open ocean OML where it is forced by the ocean losses of heat via long-wave radiation and evaporative cooling. More in general, buoyancy-driven mixing is important in the ocean for it is involved in externally forced convection (e.g., Marshall and Schott, 1999) as well as in deep internal convection (e.g., due to thermo-baric effects). Buoyancy-driven turbulence is characterized by the occurrence of localized convective plumes (e.g., Leighton et al., 2003) that present a large-spatial and temporal variability. The basic mechanism is analogous to the physics of other thermo-convective turbulent flows in wide horizontal layers like, for instance, Rayleigh–Be´nard flow. The radically different redistribution of energy into turbulent motions for the thermal convection, wind-driven shear and wind/wave-driven Langmuir turbulence cases has been recently discussed in Li and Garrett (1997), Li et al. (2004). In particular, there is shown that in a fully-developed sea state the OML is in a Langmuir turbulent regime and that even in typical convective cases (heat losses) Langmuir turbulence is rarely negligible. This very recent and important result point out on the need of further studies on the interplay among forcing mechanisms in mixing the OML. In a numerical simulation of the OML, one would ideally resolve all the flow scales, from the integral one all the way down to the smallest dissipative one, without introducing any turbulence model. Such method, called direct numerical simulation (DNS), requires using a grid size comparable to the smallest turbulence scales. It is well known that the number of grid points depends on the Reynolds and stratification (e.g., Rayleigh or Richardson) numbers. Unfortunately, the physical parameters that characterize the OML are typically such that these numbers are too high for allowing such an approach with the present or foreseeable computers. Thus, simulating the OML requires some modelling hypotheses on the resolved scales of turbulence. When modelling the OML, simple parameterizations of the vertical diffusion as a function of the vertical stratification have been found to provide good qualitative results on relatively large-scales (e.g., Umlauf and Burchard, 2005). However, such an approach does not allow for an in-depth description of the turbulent phenomena involved in the mixing, and it is unable to provide information of flow components at small scales. Therefore, several researchers have started to investigate the potentiality of large-eddy simulation (LES) (Sagaut, 2002) in simulating this type of phenomena (McWilliams et al., 1997; Wang et al., 1998; Skyllingstad et al., 1999; Skyllingstad et al., 2000; Wang, 2001; Zikanov et al., 2003; Noh et al., 2004; Min and Noh, 2004). When using an LES approach, a filtering operation is applied on the flow variables in such a way that only the large-scales, that are energetic and anisotropic, are resolved whereas the subgrid-scale (SGS) motions are modelled under general hypothesis of energy equilibrium. A substantial improvement of LES results is generally obtained by exploiting the high-accuracy of specific numerical methods (e.g., Ghosal, 1996; Kravchenko and Moin, 1997; De Stefano et al., 2001; Brown et al., 2001), along with the adoption of the dynamic SGS modelling procedure (Germano et al., 1991). However, in the case of oceanic turbulence, where strong vertical stratification can occur, the development of a suitable SGS model is still an open task. The idea behind LES is that the residual SGS motions to be modelled are sufficiently small to be considered almost isotropic. Actually, in the case of strong stratification it appears to be practically impossible to use a filter width in the vertical direction able to separate all the anisotropic large-scales of motion. As a consequence, the SGS model must account also, in some amount, of the anisotropic part of the motion that is not resolvable with a reasonable computational grid. At present, LES seems capable of reproducing results that are in good agreement with DNS data only for weak stratification. The so-called resolved LES, that is a solution in which the anisotropic scales close to a boundary and in the stratification region are all resolved, is unfeasible for real ocean parameters since this would require practically the DNS grid resolution. This difficulty is often faced by introducing ad hoc assumptions into the SGS models for taking into account the effects of strong stratification, that is using parameterizations based on the vertical stability, so that the resulting governing equations are more similar to those ones of the Reynolds Average Navier–Stokes (RANS) formulation, e.g., see (Wang et al., 1998; Skyllingstad et al., 1999; Skyllingstad et al., 2000; Wang, 2001). Conversely, a specific dynamic

F.M. Denaro et al. / Ocean Modelling 17 (2007) 199–218

201

SGS procedure, based on a second test filtering, has been proposed for resolving neutral atmosphere with the aim of taking into account the variation of the coefficients with the scale without such arbitrary assumptions (Porte´-Agel et al., 2000). This latter approach seems much more close to the spirit of the LES idea. The goal of this study is to develop, implement and analyze the performances of a LES formulation based on the integral form of the governing equations along with the dynamic SGS modelling, without using any ad hoc assumptions concerning the stratification. The general aim of the present study is to develop a numerical methodology for the study of the response of the marine upper boundary layer to time dependent internal and external forcing (e.g., the diurnal cycle). Thus, we focus on the reproduction of the onset of buoyancy-driven turbulent motions that is, indeed, a difficult exercise in LES models. In fact, the dynamic procedure is generally applied with eddy-viscosity SGS models that are based on the assumption that turbulent kinetic energy dissipation and production are in equilibrium. As already observed in Lilly (1962), it remains an open question to assess if such hypothesis is still valid in stratified flows. Being the aim here to validate a new SGS approach to turbulent buoyancy-driven motions and thus to focus on the small scale turbulent features, rotation is neglected in the present study. According to the scaling formulation for the SGS modelling of buoyancy-driven turbulence (Wong and Lilly, 1994), our study can be seen as an extension to the finite volume (FV) approach of the LES spectralbased method presented in Zikanov et al. (2002). Therein, the authors presented an extensive analysis of the spatial structure and statistical properties of turbulent thermal convection as it results from dynamic LES calculations at moderate Reynolds numbers. To validate the SGS closure, they performed a comparison with the corresponding DNS solution. Since the integral formulation leads us to reformulate the classical dynamic SGS model the same validation exercise is repeated in this work. In particular, the detailed comparison between the present FV method and that spectral one is obtained thanks to the availability of the original code kindly provided by the authors. The paper is organized as follows. In Section 2 the adopted LES methodology based on the FV numerical approach is presented. The particular definition of the SGS residual terms in both momentum and temperature equations is discussed, together with the eddy-viscosity modelling procedure. The original FV-based dynamic SGS model is presented in Section 3. The results of the application of the method to the OML onset simulation at moderate Reynolds and Rayleigh numbers values are shown in Section 4, along with a detailed spectral analysis. Section 5 is devoted to illustrate the performances of the present method when applied to the simulation of a stratified shear turbulence case, that is related to the classical laboratory experiment of Kato and Phillips (1969). Finally, in Section 6, some concluding remarks are provided.

2. FV-based LES method In this section, details about the mathematical model as well as the LES approach and turbulence modelling are given. 2.1. The mathematical model The physical model consists of a three-dimensional unsteady incompressible flow in a Cartesian domain with periodic boundary conditions along the horizontal directions (x, y). In the vertical direction (z), the upper boundary (i.e., the zero-depth level corresponding to z = H) is the interface between air and water along which heat transfer onsets the flow motion. Such interface is seen as a plane at constant atmospheric pressure, supposed also non-deformable, that corresponds to the assumption of a vanishing Froude number. The lower boundary (z = 0) stands for the seawater bottom where solid wall conditions apply. The mathematical model is based on the Navier–Stokes equations along with the Boussinesq approximation. Therefore, the density is expressed as only a function of temperature as q = q0[1 + b(#0  #)], b being the thermal expansion coefficient assumed constant. The set of governing equations is completed by the internal energy balance, written in terms of the temperature variable, where a constant specific heat (cp) is assumed and the viscous dissipation term is neglected. Thus, the governing equations are

202

F.M. Denaro et al. / Ocean Modelling 17 (2007) 199–218

rv¼0 ov 1 þ r  ðvvÞ þ rp ¼ r  ð2mSÞ  kgb½#h ðzÞ  # ot q0   o# k þ r  ðv#Þ ¼ r  r# ; ot q0 c p

ð1Þ ð2Þ ð3Þ

where S stands for the symmetric part of the velocity gradient tensor $v, that is the rate-of-strain tensor, and m = l/q0 and k are the kinematic viscosity and the thermal conductivity of the fluid, respectively. 2.2. LES formulation The above mathematical model is reformulated in the framework of LES approach based on the FV numerical methodology. The differential governing Eq. (1)–(3) are rewritten in integral form, that corresponds to perform a local volume averaging. This can be interpreted as a sort of implicit grid filtering, the corresponding filtered variable being defined as Z 1 f ðx; tÞ  f ðx0 ; tÞdx0 ; ð4Þ jXj XðxÞ with X(x) the local FV built around the point x and jXj its measure. The characteristic filter width can be as1=3 sumed to be, for instance, DðxÞ ¼ jXj . When adopting a non-uniform FV numerical grid, that corresponds to perform a LES with a variable filter width. When applying the filtering operation (4) to the continuity Eq. (1), by virtue of the divergence theorem, i.e., Z Z 1 1 rv r  v dx ¼ n  v dS ¼ 0; ð5Þ jXj X jXj oX it holds Z

n  v dS ¼ q:

ð6Þ

oX

The term q at the RHS of Eq. (6) results from the fact that averaging and divergence operators do not commute and can be interpreted as a sort of SGS term. The particular form of the LES continuity equation directly comes from the adoption of the FV approach. In fact, the volume-averaged velocity stands for the variable actually solved and each equation must be consistently expressed in terms of this one. According to Eq. (6), the resolved FV velocity field is not divergence-free, unless one uses a uniform FV-grid. In this work, a second order accurate FV discretization is exploited so that the commutation error terms result of same magnitude of the local truncation error. This way no SGS model is introduced in the continuity equation while the filtered velocity is imposed to be numerically divergence-free, that is the second order approximation q ffi 0 holds (Sagaut, 2002). As to the filtered counterpart of momentum Eq. (2), it holds Z Z Z Z     ov 1 1 1 1 þ n  v v dS þ np dS ¼ n  2mS dS  kgb #h  #  n  T dS; ð7Þ ot jXj oX q0 jXj oX jXj oX jXj oX where #h ðzÞ is the vertical temperature corresponding to the hydrostatic equilibrium, S stands for the resolved rate-of-strain tensor and the residual SGS stress term is given by T ¼ 2mðrs v  rs vÞ þ vv  v v: Finally, the filtered temperature equation is   Z Z Z o# 1 1 k 1 þ n  v# dS ¼ n r# dS  n  t dS; ot jXj oX jXj oX q0 c p jXj oX

ð8Þ

ð9Þ

F.M. Denaro et al. / Ocean Modelling 17 (2007) 199–218

203

where the heat flux SGS term is defined according to  k  r#  r# þ v#  v#: ð10Þ t¼ q0 c p The SGS terms differ from the corresponding terms in the finite difference (FD) approach due to the appearance of the viscous and thermal fluxes in Eqs. (8) and (10). The SGS model, therefore, must mimic also part of the unknown viscous and conduction terms. 2.3. SGS modelling According to the classical Bousinnesq eddy viscosity model, the SGS tensor T is modelled as an additional viscous stress-like term, where the molecular viscosity is replaced by a point-wise time-dependent eddy viscosity, noted mt. Namely, the deviatoric part of the SGS stress tensor, henceforth denoted with the superscript star, i.e., T  T  TrðTÞ I; I being the identity tensor, is modelled according to 3  s ð11Þ T ffi 2mt r v; and the new pressure variable is defined as p  p  Tr(T)/3. Similarly, the SGS heat flux vector-field t is approximated as a Fourier-like flux kt tffi r#; q0 cp kt being the SGS thermal conductivity. This way, the FV-based LES equations to be solved become Z n  v dS ¼ 0; oX

ov 1 þ ot jXj and o# 1 þ ot jXj

Z

n  ðv vÞ dS þ

oX

Z

n  v # dS ¼ oX

1 q0 jXj

1 jXj

Z

Z oX

np dS ¼

oX

 n

1 jXj

ktot r# q0 cp

Z

  n  ð2mtot rs vÞ dS  kgb #h  # ;

ð12Þ

ð13Þ ð14Þ

oX

 dS;

ð15Þ

having defined the total diffusion coefficients mtot = m + mt and ktot = k + kt. In order to close the LES problem, governed by Eqs. (13)–(15), a functional relation between the unknown model coefficients, mt and kt, and the balanced filtered fields, v and #, must be specified. Typically, the smallest flow scales are supposed to have shorter characteristic time-scales with respect to the large-eddies, so that it is possible to assume the equilibrium between production and dissipation of turbulent kinetic energy. For instance, under homothermal  1=2 flow conditions, one usually considers the standard Smagorinsky eddy viscosity model, mt ¼ CD2 2S ij S ij . For buoyancy driven flows, the temperature field being not simply a passive scalar but actively entering the flow dynamics, it is more complex to define a suitable eddy viscosity. In fact, the turbulence theory for thermally stratified flows is still less developed than for shear flows and there are no definitive research results. In particular, according to Lilly (1962), it appears doubtful that one can express the energy transfer only in terms of the mean flow, regardless of the buoyancy scales. Actually, the turbulent kinetic energy budget is altered by the buoyancy effects and, even considering still valid the equilibrium hypothesis, the eddy viscosity definition must be modified according to  1=2 o# mt ¼ CD2 2S ij S ij  a : ð16Þ oz The coefficient a in the above eddy-viscosity model measures the importance of the coupling with the temperature stratification as a function of the Richardson number (Lilly, 1962). Some authors adopted the so-called stratification formulation, e.g., (Sullivan and Moeng, 1992; Cabot, 1992), by directly expressing the eddy viscosity mt, along with the eddy diffusivity kt, in a coupled way with the buoyancy term.

204

F.M. Denaro et al. / Ocean Modelling 17 (2007) 199–218

Actually, the presence of the buoyancy term under square-root can cause non-real solutions to the iterative procedure necessary for computing the model coefficients. In order to overcome the problem, a different procedure, in the framework of the so-called scaling formulation, was proposed in Wong and Lilly (1994). In fact, according to the Kolmogorov scaling, the eddy viscosity can be expressed in terms of the grid-filter width D and the dissipation rate e as mt ¼ C 2=3 D4=3 e1=3 ;

ð17Þ

where C stands for the Smagorinsky model coefficient. This way, the dissipation rate e is no more assumed to balance the SGS energy production due to the buoyancy effect as well as to other possible forcing terms. The dissipation rate is dynamically evaluated according to the procedure discussed in the following. The main goal of the present work is to analyze the scaling formulation in the framework of a FV-based LES approach, by exploring its capability in providing accurate results for ocean turbulence, in case of buoyancy effects. It is hoped that this model formulation can automatically take into account for the presence of other forcing terms in the momentum equation, like the salinity buoyancy as well as the Langmuir circulation (Noh et al., 2004). Thus, based upon (17), the SGS model (11) rewrites as T ffi 2C e D4=3 rs v;

ð18Þ

where Ce  C2/3e1/3 stands for the model coefficient to be actually determined. 3. FV-based dynamic procedure The method is here completed by developing an original dynamic procedure for determining the SGS model coefficients suited to our FV-based LES approach. In fact, the integral form of the governing equations requires a suitable treatment of the SGS terms, as it is illustrated in the following. The dynamic procedure is based on the application of a second (test) filter with a greater width, say e ¼ /Dð/ > 1Þ. In the framework of the FV methodology that corresponds to apply a local average over D an extended spatial domain. Namely, the test filtering is defined according to Z 1 ~  ð19Þ f ðx; tÞ ¼ f ðx0 ; tÞ dx0 ; e e j Xj X ðxÞ S e XðxÞ ¼ k Xk ðxÞ being a suitable extended local volume surrounding x. When using a non uniform vertical grid, a two-dimensional uniform test-filter is applied in the horizontal homogeneity plane. Let us start illustrating the dynamic procedure when applied to the momentum equation. For the sake of brevity, let us consider the LES Eq. (7) in the following filtered differential form   ov 1 þ r  v v þ rp ¼ r  ð2mrs vÞ  r  T  kgb #h  # : ð20Þ ot q0 On the other hand, when applying the test filter one gets:   1     g oev e ; f ¼ r g þ r  ev ev þ r 2mrs ev  rg  R  kgb f #h  # p ot q0 where, in analogy to definition (8), the residual stress tensor with respect to the test filter is   R ¼ 2m rs v  rs ev þ vv  ev ev:

ð21Þ

ð22Þ

The definitions (8) and (22) allow us to generalize the Germano identity as R  T ¼ L; with the following definition for the known Leonard stress term   L  2m rs v  rs ev þ v v  ev ev:

ð23Þ

ð24Þ

F.M. Denaro et al. / Ocean Modelling 17 (2007) 199–218

205

In analogy to model (11), let us use a similar approximation for R, that is, R ffi 2m0t rs~v; m0t being the eddy viscosity related to the test filtering procedure. This way, the identity (23) becomes 2mt rs v  2m0t rs ev ffi L :

ð25Þ

It is worthwhile to stress the particular form assumed by the Germano identity, when compared to the standard one, corresponding to the adoption of the FD numerical approach. In fact, the SGS stress tensor T does not appear under test filtering and, congruently, the definition of the known Leonard tensor (24), apart from the presence of the molecular diffusion terms, is rather different. This specific consequence of the FV integral formulation appears to be very appealing as it allows us to avoid the usual approximation with which the model coefficient is extracted from the filtering operation. Starting from Eq. (25), the modelling coefficient Ce is dynamically determined as it follows. By assuming the e 4=3 ¼ C /4=3 D4=3 , it holds same rate of dissipation for both filters, that is m ¼ C D4=3 and m0 ¼ C D t

e

t

e

e

h i 2mt rs v  /4=3 rs ev ffi L :

ð26Þ

Eq. (26) stands for a system of five independent equations for the unique scalar unknown Ce. By projecting Eq. (26) along some particular direction, one can obtain one equation for Ce. For instance, similarly to Zikanov et al. (2002), one can consider 2mt ¼ h

L : rs ev i : rs v  /4=3 rs ev : rs ev

ð27Þ

As far as the temperature equation is concerned, again one can adopt a dynamic approach for determining the eddy diffusivity kt. In fact, the integral temperature Eq. (9) is rewritten in the following grid filtered differential form:     o# k þ r  v# ¼ r  r#  r  t; ð28Þ ot q0 c p and, when applying the test filter, g  e g o# k e e g e þ r  ð v #Þ ¼ r  r#  r  r: ot q0 c p

ð29Þ

In Eq. (29), the SGS heat flux corresponding to test filtering is  k  e þ v#  ev #: e r¼ r#  r # q0 c p

ð30Þ

Also, based upon (10) and (30), one obtains the new Germano identity r  t ¼ l;

ð31Þ

where the following known vector field     e  k r#  r # e l ¼ v#  ev # q0 c p is defined. k0 ~ with Similarly to (12), let us consider the approximation r ffi  q0ct p r#, scaling formulation, e 4=3 k0t Ce D ¼ ; p0 cp Prt where Prt = q0mtcp/kt is the turbulent Prandtl number.

ð32Þ k0t q0 cp

m0

¼ Prtt . That is, according to the

ð33Þ

206

F.M. Denaro et al. / Ocean Modelling 17 (2007) 199–218

This way, owing to Eq. (31), it holds h i e ffil C d r#  /4=3 r #

ð34Þ

C e 4=3 where C d  Pr D stands for a new modelling coefficient. t Actually, Eq. (34) represents a system of three independent equations for the unknown Cd. Again, the solution can be determined by projecting Eq. (34) along some particular direction. For instance, the model coefficient can be determined according to

Cd ¼ h

e l  r# i : e  r# e r#  /4=3 r #

ð35Þ

The LES governing equations are made dimensionless by choosing proper length (L) and velocity (U) reference quantities. As a consequence, time is non-dimensionalized by L/U and pressure field by q0U2. The non-dimensional temperature variable is defined as ð#  #0 Þ=D#; D# being the temperature difference between bottom plane and interface corresponding to the initial stratification. This way, a set of characteristic non-dimensional numbers is introduced. Namely, the Reynolds number, Re = UL/m, the molecular Prandtl number, Pr = q0mcp/k, the Rayleigh number, Ra = gbLD#/U2 and the Nusselt number, Nu = (q*L)/(kD#). Finally, the non-dimensional FV-based LES equations – for the sake of simplicity written without introducing new variables – are recast as Z n  v dS ¼ 0; ð36Þ oX

ov 1 þ ot jXj

Z

1 n  ðvvÞ dS þ jXj oX

Z

1 np dS ¼ jXj oX

Z



 2 2 þ n rs v dS  kRað#h  #Þ; Re Ret oXðxÞ

ð37Þ

and o# 1 þ ot jXj

Z

n  ðv#Þ dS ¼ oX

1 jXj



Z

n oX

 1 1 þ r# dS: RePr Ret Prt

ð38Þ

In the above equation, Ret = UL/mt is the turbulent Reynolds number. This way, the non-dimensional counterparts of the SGS model coefficients, (27) and (35), become 2 L : rs~v i ¼h Ret rs v  /4=3 rs~v : rs~v

ð39Þ

e 1 1  r# i ¼h : e  r# ~ Ret Prt r#  /4=3 r #

ð40Þ

and

As usual in dynamic LES calculations, both the numerator and denominator of Eqs. (39) and (40) are spatially averaged in the homogeneous directions. 4. Numerical validation In order for the proposed FV-based LES methodology to be validated, the simulation of turbulent thermal convection generated by surface cooling in a finite-depth stably stratified horizontal layer is considered. This flow can be seen as a simplified model of the OML formation in shallow water, due to a steady surface cooling at the air–sea interface, in absence of wind stress and in a non-rotating framework. The same phenomenology has been examined in Zikanov et al. (2002), where an hybrid numerical method has been exploited, by applying a pseudo-spectral method in the horizontal homogeneity plane and a FD discretization along the vertical

F.M. Denaro et al. / Ocean Modelling 17 (2007) 199–218

207

direction. According to this work, no-slip velocity condition and constant temperature are assigned at the bottom boundary (that is, z = 0). In this section, only moderate values of the non-dimensional characteristic numbers are fixed in order for a DNS solution to be feasible. A real-world set of non-dimensional numbers will be used in the flow problem illustrated in the next section. As to the numerical implementation, the time integration of the LES Eqs. (36)–(38), supplied with the dynamic model previously described, is based on a second order time-accurate projection-method, (e.g., see Denaro, 2005). The non-dimensional computational domain is partitioned by means of a Cartesian structured non-uniform grid. Namely, the grid is uniform in the horizontal plane and stretched along the vertical direction, according to a cosine-hill distribution. Furthermore, a non-staggered cell-centered collocation of the variables is adopted. Thus, the differential operators are approximated by means of second-order discrete operators, while the surface integrals are discretized by exploiting the mean value formula based on linear interpolation. 4.1. Case settings The simulations start with a zero-velocity field and a weak initial vertical stratification. A negative heat flux is prescribed at the air–sea interface, providing a loss of sea internal energy and the onset of motion. Therefore, the flow motion is generated solely by inhomogeneities in the temperature field. The actual flow dynamics can be considered in the framework of the Rayleigh–Be´nard convective regime. Namely, after a first stage dominated by molecular temperature diffusion, that causes the development of a cool thermal boundary layer, the denser fluid layer becomes unstable and convective plumes are generated (Zikanov et al., 2002; Leighton et al., 2003). For instance, this flow can be considered a simple prototypal of convection occurring during nighttime under calm clear-sky open-ocean conditions, for which the long-wave radiation dominates the surface heat flux. As reference to the non-dimensional parameters, the depth H is chosen as the reference length so that the non-dimensional depth has unitary measure and the prescribed horizontal extension is p/2  p/2. In the present case of purely buoyancy driven convection, the reference velocity is fixed as U  (B0H)1/3, B0 being the buoyancy flux at the interface, agjq*j/(q0cp) (e.g., see Maxworthy, 1997). The simulations have been carried out for moderate Reynolds and Rayleigh numbers, namely Re = 1200 and Ra = 2. The fixed heat flux at the interface is such that Nu = 600, while the fluid Prandtl number is assumed to be Pr = 1. Such flow parameters, quite far from realistic marine values, have been chosen in order to be able to realize DNS solutions on a reasonable numerical grid. Actually, at this stage, our main interest is in the verification of the proposed numerical method, rather than in the prediction of realistic ocean flows. However, one can think of possible corresponding physical parameters as the following: a sea depth H = 10 m and a reference velocity of about 104 m/s, that leads to a characteristic time-scale of about 27 h. The initial temperature field is characterized by D# 106 K. Finally, the performances of the FV-based LES method will be explored in the following for a more realistic geophysical flow, that is the OML deepening. The FV-based DNS and LES solutions are compared to those obtained with the pseudo-spectral code of Zikanov et al. (2002), kindly provided by the authors, that have been used as reference. 4.2. DNS solution A computational grid of 1003 cells is adopted for obtaining a FV-based DNS solution. In order for the thermal boundary layer to be fully resolved, the computational cells are clustered along the vertical direction close to the air–sea interface. The reference spectral solution is obtained with a 1922 uniform grid in the horizontal plane and 150 grid-points non-uniformly distributed in the vertical direction, that are clustered on both sides. The temporal evolution of the volume-averaged kinetic energy for both simulations is reported in Fig. 1 up to non-dimensional time t = 9, the reference time-scale being H/U = 105 s. The curves show the first stage dominated by the thermal diffusion, that is followed by the onset of temperature plumes. That corresponds to the potential energy falling down, that is not reported among the figures for sake of brevity.

208

F.M. Denaro et al. / Ocean Modelling 17 (2007) 199–218 0.7 DNS FV DNS SM 0.6

0.5

0.4

0.3

0.2

0.1

0 0

1

2

3

4

5

6

7

8

9

10

Time Fig. 1. DNS: temporal evolution of kinetic energy.

Clearly, the two DNS solutions cannot be compared in strictly a deterministic meaning. Indeed, the initial conditions, that are obtained by randomly perturbing the field, are slightly different and the numerical errors act differently. Thus, it is not surprising that the initial flow instability starts at different time instants. Anyway, the plots seems substantially similar, from both a qualitative and a quantitative point of view. In particular, the energy magnitudes appear absolutely comparable during the evolution. Of course, for this type of boundary conditions, there is no statistical equilibrium and the internal energy is continuously converted into kinetic one, as it clearly appears from the picture.

DNS SM DNS FV Inertial slope

0.1

0.01

0.001

0.0001

1e-005

1e-006

1e-007

1e-008 10

100

Fig. 2. DNS: spectral distribution of kinetic energy for a horizontal velocity component (v) versus the one-dimensional wavenumber kx.

F.M. Denaro et al. / Ocean Modelling 17 (2007) 199–218

209

In fact, a fundamental tool for effectively comparing the results requires performing statistical analysis of the solution. Here, a low order statistics is considered by simply analyzing the spectral distribution in the horizontal plane of the instantaneous kinetic and thermal energy. Namely, the 1-D spectra are obtained for a prescribed depth, i.e., z = 0.75, by transforming in one horizontal direction and averaging the Fourier coefficients in the other one. However, owing to the different initial transition, in absence of a statistical steady state, there is some difficulty in individuating a specific time for the analysis. Here, the time instant t = 6 has been chosen because, as it appears from the inspection of Fig. 1, both the energy levels and the energy rates are comparable.

DNS SM DNS FV Inertial slope

0.1

0.01

0.001

0.0001

1e-005

1e-006

1e-007

1e-008 10

100

Fig. 3. DNS: spectral distribution of kinetic energy for the vertical velocity (w) versus the one-dimensional wavenumber kx.

DNS SM DNS FV Inertial slope

0.1

0.01

0.001

0.0001

1e-005

1e-006

1e-007

1e-008 10

100

Fig. 4. DNS: spectral distribution of thermal energy versus the one-dimensional wavenumber kx.

210

F.M. Denaro et al. / Ocean Modelling 17 (2007) 199–218

The instantaneous 1-D kinetic energy spectra for one horizontal velocity component and the vertical one, along one horizontal direction, are reported in Figs. 2 and 3, respectively, for both FV and spectral solutions. As to the temperature field, the thermal energy spectra at the same time-instant and depth are shown in Fig. 4. The theoretical slope (5/3) of the inertial range is also reported. Both kinetic and thermal energy spectra show a quite short inertial range owing to the moderate Reynolds and Rayleigh number values. However,

Fig. 5. DNS: visualization of temperature field in terms of iso-surfaces at t = 2 (top), t = 4 (middle) and t = 6 (bottom). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

F.M. Denaro et al. / Ocean Modelling 17 (2007) 199–218

211

in spite of the lower resolution, the FV-based DNS results appear similar to spectral ones, especially in the large-scale range. Some oscillations are caused by the absence of a time-averaging procedure while smaller discrepancies also depend upon the different initialization. Anyway, the comparison can be considered fully satisfactory and the FV-based numerical procedure successfully validated. In order to visualize the flow evolution, in Fig. 5 the temperature iso-surfaces at two different values h = 3 (cold fluid zones in blue) and h = 0.5 (hot fluid zones in red), are reported at the three time instants t = 2, 4 and 6, showing the typical development of thermal plumes. 4.3. LES solution The above DNS solution provides an indication of the inertial range extension so that one can correspondingly fix a suitable LES filter width. In fact, the FV-based LES solution is obtained on a 322  100 numerical grid, by supplying the simulation with the FV-based dynamic model introduced in Section 3. The grid-to-filter ratio /, is fixed to 2, the discrete test-filter being constructed by averaging in the horizontal plane over a 3  3 FV-stencil. Note that the SGS viscosity and diffusivity coefficients are clipped so that they stay non-negative. In order to highlight the effectiveness of the dynamic SGS model, a FV simulation without any SGS model is also performed at the same numerical resolution. For sake of completeness, a static Smagorinsky LES solution is obtained by prescribing C = 0.044 in the stratification formulation (16) as well as Prt = 0.4 for the LES temperature equation, that are the same values used in Wong and Lilly (1994). The reference spectral LES solution is obtained on a 482  150 numerical grid. This way, considered the enlarged horizontal extension of the computational domain, the spatial resolution is practically the same as in Zikanov et al. (2002), where an additional low-pass filtering is applied, periodically in time, in order to obtain a numerically stable solution. For a discussion of why and how this explicit filtering is performed, one can see the original work (Zikanov et al., 2002). On the other hand, as demonstrated by using no SGS modelling procedure, the FV-based solution is stable by itself, thanks to the built-in low-pass filtering effect of the numerical FV discretization. The temporal evolution of the volume-averaged kinetic energy is shown in Fig. 6 for the different LES solutions. Apart a time delay for the onset of flow instability, the FV-based and spectral solutions substantially 0.7 LES FV dynamic LES FV static LES FV no-model LES SM

0.6

0.5

0.4

0.3

0.2

0.1

0 0

1

2

3

4

5

Time Fig. 6. LES: temporal evolution of kinetic energy.

6

212

F.M. Denaro et al. / Ocean Modelling 17 (2007) 199–218

agree. The no-model solution shows a greater initial energy peak, corresponding to a non-physical potential energy dropping, followed by a less energized evolution. This behavior demonstrates that, despite of the moderate values of the Reynolds and Rayleigh numbers, this case is meaningful for LES testing. In particular, the use of an efficient SGS model is unavoidable. The 1-D kinetic energy spectra along one horizontal direction are reported in Figs. 7 and 8 for both a horizontal velocity component and the vertical one, respectively. The thermal energy spectra at the same timeinstant and depth are shown in Fig. 9. The FV DNS spectra are also reported in the figures, as a natural reference. Note that, for a better comparison, the DNS spectra should have been a-posteriori filtered with DNS FV LES FV dynamic LES FV static LES FV no-model LES SM Inertial slope

0.1

0.01

Log(Ek)

0.001

0.0001

1e-005

1e-006

1e-007

1e-008 10

100

Log(k) Fig. 7. LES: spectral distribution of kinetic energy for a horizontal velocity component (v) versus the one-dimensional wavenumber kx.

DNS FV LES FV dynamic LES FV static LES FV no-model LES SM Inertial slope

0.1

0.01

Log(Ek)

0.001

0.0001

1e-005

1e-006

1e-007

1e-008 10

100

Log(k) Fig. 8. LES: spectral distribution of kinetic energy for the vertical velocity (w) versus the one-dimensional wavenumber kx.

F.M. Denaro et al. / Ocean Modelling 17 (2007) 199–218

213

DNS FV LES FV dynamic LES FV static LES FV no-model LES SM Inertial slope

0.1

0.01

Log(Ek)

0.001

0.0001

1e-005

1e-006

1e-007

1e-008 10

100

Log(k) Fig. 9. LES: spectral distribution of thermal energy versus the one-dimensional wavenumber kx.

the same FV-filter corresponding to the LES solution. The effect of the implicit smooth FV-filter is evident when considering the diminished energy content with respect to the DNS solution. The difference in the maximum resolved wavenumbers for the different LES solutions trivially reflects the different adopted numerical grid. When no SGS model is used the energy spectrum shows the typical pile-up corresponding to the smallest resolved scales. That causes a clear deviation from the inertial scaling, again highlighting the necessity of an SGS model able to drain energy from the resolved field. Also, as to comparison between FV-based dynamic and static Smagorinsky models, it appears evident that the prescribed model coefficients for the static version are not well suited to the actual transitional flow conditions. In fact, the energy spectra demonstrate that the static Smagorinsky model is over-dissipative whereas the dynamic procedure provides a good estimation for the required SGS dissipation. Furthermore, the FV-based dynamic LES is similar to the spectral one, even if the former is computed with a lower grid-resolution. The effect of the additional explicit low-pass filtering on the spectral solution appears evident at highest resolved wavenumbers. As pointed out by the authors, this explicit filter must be considered as an effective tool of de-aliasing (Lele, 1992) since no explicit de-aliasing procedure has been performed in order to save computational resources. Actually, the spectral LES revealed numerically unstable when this additional filter is suppressed, whereas the FV-based LES remains however stable. In conclusion, the FVbased LES method is able to better reproduce the kinetic energy content of the flow at highest resolved wavenumbers. By inspection of the thermal energy content illustrated in Fig. 9, a slight excess of dissipation produced by the SGS model is evident. That could be due to the adoption of a pure eddy-viscosity model instead of a mixed one (that is with an additional scale-similarity term). Moreover, the adopted clipping procedure could lead to over-estimate the energy dissipation as the energy backscatter is not allowed. The simulation of this reverse energy cascade would be important mainly close to the air–sea interface, where a two-dimensional vorticity field is approached. However, such issues will be further explored in the future. 5. An example of application to geophysical flows: the Kato–Phillips test case In the previous section the proposed LES methodology has been validated via a comparison with a DNS solution for a convective turbulence case at low-Reynolds and Rayleigh numbers. Of course, the illustrated

214

F.M. Denaro et al. / Ocean Modelling 17 (2007) 199–218

problem is not exhaustive of typical phenomena occurring in the ocean turbulence, for both the reduced scale model and the absence of other driving force (e.g., wind stress, rotation, etc.). However, owing to the fact that in real geophysical flows the characteristic numbers are several magnitude orders greater, a DNS calculation is unfeasible because of the unaffordable computational cost. Therefore, after having used the reduced scale model for validation purposes, a flow problem more representative of the oceanic turbulence is now addressed. The present FV-based LES model is applied to the simulation of a stratified shear turbulence case, that is related to the classical laboratory experiment of Kato and Phillips (1969), now almost a standard reference for the ocean modeling community. In that experiment, the surface shear-stress caused by a driven belt mixes an initially linearly stratified quiescent fluid in a laboratory tank thus developing a turbulent boundary layer that mimics the initial deepening of the OML, e.g., (Burchard, 2001; Kundu, 1980; Denman, 1973). As to numerical experiment setting, the computational domain has an horizontal extension of 30 m  30 m and a depth of 50 m, for a set of typical ocean parameters, that are the initial stratification N 20 ¼ gb o#ozh and the superficial shear-stress s ¼ q0 u2 ; u ¼ 0:0127 m=s being the friction velocity. Furthermore, adiabatic boundary condition is assumed at the air–sea interface. In this case, as Ra is inappropriate for shear-driven turbulence, the initial Richardson number based on the friction velocity is used, that is Ri0 ¼ N 20 L2 =u2 . The following flow parameters have been prescribed for the numerical experiments: Re = 6.35  105, Ri0 = 1.55  103, Pr = 6.7 and Nu = 0. For such a flow configuration a well-known analytical correlation exists, i.e., according to Price (Price, 1979) (see also in Umlauf and Burchard (2005)), the OML depth h varies in time according to 1

h ¼ 1:047u ðt=N 0 Þ2 :

ð41Þ

In fact, Eq. (41) is valid only for non-rotating flows, i.e. for time-scales shorter than the Coriolis parameter f1 while, for larger times, the rotational effects become significant and the deepening process tends to arrest when a stratified Ekman turbulent layer has developed to some extent (e.g., Zikanov et al., 2003). However, since at mid-latitudes it holds f = 104 s1, the non-rotating solution can be however applied for the first few hours. The prediction of the OML deepening is an important requirement for an ocean turbulence model, the OML depth definition in (41) being however not obvious. Although different definitions exist in the literature the OML is generally defined as the region of uniform fluid properties just below the ocean surface. This definition implies the examination of tracer values, e.g., the fluid density. For instance, one can define the OML depth as the location corresponding to the maximum vertical density or as the depth where the density reaches a given threshold value, e.g., (de Boyer Monte´gut et al., 2004). Another possible definition consists in identifying the lower boundary of the region of high-turbulence intensity (turbocline). However, for shear turbulence in stratified fluids the mostly used parameter for describing the turbulence state is the gradient Richardson number, that is the dimensionless ratio of the buoyant production/consumption of turbulence and the turbulence shear production Ri ¼ N 2

 2 du : dz

ð42Þ

This ratio can be therefore used to characterize the flow instability and the formation of turbulence and the OML depth can be identified as the depth where Ri takes a prescribed critical value, even if the choice of the latter is still an open issue. Here, following a common criterion, we define Ri = 1 as the critical Richardson number. Different FV-based dynamic LES solutions of the Kato–Phillips test-case have been obtained for different numerical parameters. Namely, the numerical solution has been analysed for different grid-to-filter ratios, /, and for both uniform and non-uniform vertical grids made of 100 grid-points, while maintaining a 32  32 horizontal grid resolution. The choice of the vertical grid is crucial for such an application to shear turbulence simulation. For realistic oceanic conditions a uniform grid is usually preferred, in order to save computational resources, even though local grid refinement is commonly used for flows of engineering interest. A further comparison is made with a static version of the Smagorinsky model, that is by a-priori prescribing the model coefficients as C = 0.1 and Prt = 0.4. All the solutions have been obtained for the first 2 h of the OML deepening.

215

Z

Z

F.M. Denaro et al. / Ocean Modelling 17 (2007) 199–218

Ri

Z

T

dT/dz Fig. 10. OML depth estimated with different criteria (case FGR = 4, uniform grid). Upper left panel: vertical profile of Ri. Upper right panel: vertical profile of T. Lower left panel: vertical gradient of T. Superimposed are the theoretical OML depth (solid line) and the estimated OML depth (dashed line) calculated for Ri = 1 (upper left panel) or the position of maximum vertical gradient of T (other panels).

First, some possible criteria for determining the OML depth are illustrated in Fig. 10. The top left picture shows the vertical profile of Ri, computed upon the horizontally averaged fields, for an example solution. As expected, it stays approximately constant in the upper region and then sharply increases at a certain depth. The horizontal dashed line corresponds to the OML depth based upon the Ri = 1 criterion that well matches the analytical one (horizontal solid line) calculated from (41). In contrast, by choosing Ri = 0.25, as generally accepted, one gets into the core of the OML, as shown by inspection of the temperature profile illustrated in the top right picture. Finally, the criterion based upon the maximum temperature (or, equivalently, density) gradient clearly over-estimates the OML depth (see the bottom left plot). The theoretical OML deepening is well captured by the dynamic LES solutions, as reported in Fig. 11, while the static Smagorinsky model does not work properly. The vertical temperature profiles after 2 hours are shown in Fig. 12. The LES results are reported for different solutions, on uniform and non-uniform grids with filter-to-grid ratio FGR  / = 2 and 4. The existence of a well mixed region below the air–sea surface is evident, while it is worthwhile highlighting that the run with / = 2 produces a more rapid cooling, irrespective of the vertical grid resolution. As to spectral analysis, the 1-D energy spectra for the longitudinal velocity versus the wavenumber kx are depicted in Fig. 13.

216

F.M. Denaro et al. / Ocean Modelling 17 (2007) 199–218

FGR=4, Uniform FGR=4, Non-uniform FGR=2, Uniform Static, Uniform

Fig. 11. OML deepening: time evolution of the depth for different LES solutions, on uniform and non-uniform grids with filter-to-grid ratio FGR  / = 2 and 4, compared to the static Smagorinsky model.

1 uniform grid - FGR = 4 non uniform grid - FGR = 4 uniform grid - FGR = 2

0.8

z

0.6

0.4

0.2

0 0

0.2

0.4

0.6

0.8

1

Temperature Fig. 12. OML deepening: temperature profile for different LES solutions, on uniform and non-uniform grids with filter-to-grid ratio FGR  / = 2 and 4.

One can generally conclude that the FV LES method provides quite good results, even for the adopted coarse horizontal grid and second-order accuracy of the numerical method. In particular, the adoption of non-uniform grid-spacing at the air–sea interface makes it possible to better represent the boundary heat flux and shear-stress. However, the filter-to-grid ratio plays a very important role. For a given numerical resolution, that is for a fixed numerical grid-spacing, a better evaluation of the Germano identity can be reached by increasing /, that

F.M. Denaro et al. / Ocean Modelling 17 (2007) 199–218

217

0.1 uniform grid - 4 non uniform grid - 4 uniform grid - 2 -5/3 slope

Ek11

0.01

0.001

0.0001

1e-005 10

100

Kx Fig. 13. OML deepening: spectral distribution of kinetic energy for the longitudinal velocity component (v) versus the one-dimensional wavenumber kx.

leads to a better determination of the unknown SGS stress terms. In fact, the local truncation error of the numerical FV scheme strongly affects the highest part of the resolved spectrum and the numerical errors become unacceptable also in the range of resolved large-scales. For larger /, this unavoidable contamination of the LES solution by the numerics is drastically reduced. 6. Concluding remarks A new dynamic FV-based LES method has been introduced and validated for numerical simulation of the OML onset. A novel parameterization of turbulent small-scale motions has been tested when coupled to the FV approach. A simple physical model, useful for a deeper understanding of buoyancy-driven mixing occurring in the ocean, has been considered where the turbulent convection is generated by cooling the upper surface of a layer initially stratified, while maintaining the lower surface at constant temperature. That is a rude but insightful approach to modelling the marine night-time cooling in calm weather. Several simulations have been performed to verify the FV approach. In particular, DNS and LES solutions have been obtained at moderate Reynolds and Rayleigh numbers and compared to the spectral reference ones (Zikanov et al., 2002). Since no spectral analysis was presented in that original work, the spectral DNS and LES energy spectra presented here can be seen as a fundamental result of the present study. An interesting result is obtained when comparing the proposed dynamic SGS model to the corresponding static version. It appears difficult, if not impossible, to a-priori prescribe the model coefficient in such a way that the SGS model well behaves for the different stages of the flow evolution. In fact, the dynamic model allows us to follow the onset of the physical flow instabilities and their further development without ad hoc prescriptions. On the other hand, the LES no-model clearly highlights that without any modelling the results are affected by non-physical energy scaling that shows an increasing of the energy content close to the cut-off. The method is further applied to the simulation of a well-known sheared and stably-stratified turbulence case, that is the OML deepening. Finally, it is worth emphasizing that, although the results presented have to be considered preliminary, nevertheless one can conclude that the proposed model is very promising for the study of geophysical flows.

218

F.M. Denaro et al. / Ocean Modelling 17 (2007) 199–218

Acknowledgements The authors express their gratitude to Prof. Oleg Zikanov for having provided the spectral code with which reference solutions have been obtained. Also, the authors thank the HPC center CINECA for granted calculation time on IBM SP Power 5. In addition, we wish to thank the two reviewers for constructive recommendations and important hints. References Brown, D.L., Cortez, R., Minion, M.L., 2001. Accurate projection methods for the incompressible Navier–Stokes equations. J. Comput. Phys. 168, 464–499. Burchard, H., 2001. Note on the q2l equation by Mellor and Yamada [1982]. J. Phys. Oceanogr. 31, 1377–1387. Cabot, W., 1992. Large eddy simulations of time-dependent and buoyancy-driven channel flows, in: Proc. of CTR Annual Research Briefs, Center for Turbulence Research 1992. pp. 111–122. de Boyer Monte´gut, C., Madec, G., Fischer, A.S., Lazar, A., Iudicone, D., 2004. Mixed layer depth over the global ocean: An examination of profile data and a profile-based climatology. J. Geophys. Res. 109, C12003. Denaro, F.M., 2005. Time-accurate intermediate boundary conditions for large-eddy simulation based on projection methods. Int. J. Numer. Meth. Fluids 48, 869–908. Denman, K.L., 1973. A time-dependent model of the upper ocean. J. Phys. Oceanogr. 3, 173–184. De Stefano, G., Denaro, F.M., Riccardi, G., 2001. High order filtering for control volume flow simulations. Int. J. Numer. Methods Fluids 37, 797–835. Germano, M., Piomelli, U., Moin, P., Cabot, W.H., 1991. A dynamic subgrid-scale eddy viscosity model. Phys. Fluids A 3 (7), 1760–1765. Ghosal, S., 1996. An analysis of numerical errors in large eddy simulations of turbulence. J. Comput. Phys. 125, 187–206. Kato, H., Phillips, O.M., 1969. On the penetration of turbulent layer into stratified fluid. J. Fluid Mech. 37, 643–655. Kravchenko, A.G., Moin, P., 1997. On the effect of numerical errors in large eddy simulations of turbulent flows. J. Comput. Phys. 131, 310–322. Kundu, P.K., 1980. A numerical investigation of mixed layer dynamics. J. Phys. Oceanogr. 10, 220–236. Leighton, R.I., Smith, G.B., Handler, R.A., 2003. Direct numerical simulations of free convection beneath an air–water interface at low Rayleigh numbers. Phys. Fluids 15, 3181–3193. Lele, S.K., 1992. Compact finite difference schemes with spectral-like resolution. J. Comput. Phys. 103, 16–42. Li, M., Garrett, C., 1997. Mixed layer deepening due to Langmuir circulation. J. Phys. Oceanogr. 27, 121–132. Li, M., Garrett, C., Skyllingstad, E., 2004. A regime diagram for classifying turbulent large eddies in the upper ocean. Deep-Sea Res. 52, 259–278. Lilly, D.K., 1962. On the numerical simulation of buoyant convection. Tellus 14 (2), 148–172. Marshall, J., Schott, F., 1999. Open-ocean convection: observations, theory and models. Rev. Geophys. 37, 1–64. Maxworthy, T., 1997. Convection into domains with open boundaries. Annu. Rev. Fluid Mech. 29, 327–371. McWilliams, J.C., Sullivan, P.P., Moeng, C., 1997. Langmuir turbulence in the ocean. J. Fluid Mech. 334, 1–30. Min, H.S., Noh, Y., 2004. Influence of the surface heating on langmuir circulation. J. Phys. Oceanogr. 34, 2630–2641. Noh, Y., Min, H.S., Raasch, S., 2004. Large-eddy simulation of the ocean mixed layer: the effects of wave breaking and Langmuir circulation. J. Phys. Oceanogr. 34, 720–735. Porte´-Agel, F., Meneveau, C., Parlange, M.B., 2000. A scale-dependent dynamic model for large-eddy simulation: application to neutral atmospheric boundary layer. J. Fluid Mech. 415, 261–284. Price, J.F., 1979. On the scale of stress driven entrainment experiment. J. Fluid Mech. 90, 509–529. Sagaut, P., 2002. In: Large Eddy Simulation for incompressible flows. An introduction. Springer. Skyllingstad, E.D., Smyth, W.D., Mourn, J.N., Wijesekera, H., 1999. Upper-ocean turbulence during a westerly wind burst: a comparison of large-eddy simulation results and microstructure measurements. J. Phys. Oceanogr. 29, 5–28. Skyllingstad, E.D., Smyth, W.D., Crawford, G.B., 2000. Resonant wind-driven mixing in the ocean boundary layer. J. Phys. Oceanogr. 30, 1866–1890. Sullivan, P., Moeng, C.H., 1992. An evaluation of the dynamic subgrid scale model in buoyancy driven flows. In Proc. 10th Symp. Turbulence and Diffusion. Portland, Oregon. 1992. Umlauf, L., Burchard, H., 2005. Second-order turbulence closure models for geophysical boundary layers, a review of recent work. Cont. Shelf Res. 25, 795–827. Wang, D., 2001. Large-eddy simulation of the diurnal cycle of oceanic boundary layer: sensitivity to domain size and spatial resolution. J. Geophys. Res. 106, 13959–13974. Wang, D., McWilliams, J.C., Large, W.G., 1998. Large-eddy simulation of the diurnal cycle of deep equatorial turbulence. J. Phys. Oceanogr. 28, 129–148. Wong, V.C., Lilly, D.K., 1994. A comparison of two dynamic subgrid closure methods for turbulent thermal convection. Phys. Fluids 6, 1016–1023. Zikanov, O., Slinn, D.N., Dhanak, M.R., 2002. Turbulent convection driven by surface cooling in shallow water. J. Fluid Mech. 464, 81–111. Zikanov, O., Slinn, D.N., Dhanak, M.R., 2003. Large-eddy simulations of wind-induced turbulent ekman layer. J. Fluid Mech. 495, 343–368.