- Email: [email protected]

A ﬁnite volume dynamic large-eddy simulation method for buoyancy driven turbulent geophysical ﬂows 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 ﬁnite volume approach is developed and evaluated for buoyancy driven turbulence in the absence of rotational eﬀects. The eﬀect of the sub-grid scale motions is modelled by using the dynamic scaling formulation, without introducing any ad hoc assumption for the vertical stratiﬁcation. This way, the eddy coeﬃcients for momentum and temperature equations are independently computed, while avoiding to introduce unjustiﬁed buoyancy production terms in the constitutive equations. As a result, the computational cost is considerably reduced with respect to the classical stratiﬁcation formulation. The method is presented in detail, by stressing the particular features of the ﬁnite 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 ﬁnite-depth stably stratiﬁed 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 stratiﬁed ﬂuid. Ó 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 ﬂuxes of heat and momentum as well as by surface waves. The shear due to wind-driven currents can eﬃciently 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 eﬀects). 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 ﬂows in wide horizontal layers like, for instance, Rayleigh–Be´nard ﬂow. The radically diﬀerent 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 ﬂow 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 stratiﬁcation (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 diﬀusion as a function of the vertical stratiﬁcation 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 ﬂow 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 ﬁltering operation is applied on the ﬂow 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 speciﬁc 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 stratiﬁcation 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 suﬃciently small to be considered almost isotropic. Actually, in the case of strong stratiﬁcation it appears to be practically impossible to use a ﬁlter 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 stratiﬁcation. The so-called resolved LES, that is a solution in which the anisotropic scales close to a boundary and in the stratiﬁcation region are all resolved, is unfeasible for real ocean parameters since this would require practically the DNS grid resolution. This diﬃculty is often faced by introducing ad hoc assumptions into the SGS models for taking into account the eﬀects of strong stratiﬁcation, 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 speciﬁc dynamic

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

201

SGS procedure, based on a second test ﬁltering, has been proposed for resolving neutral atmosphere with the aim of taking into account the variation of the coeﬃcients 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 stratiﬁcation. 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 diﬃcult 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 stratiﬁed ﬂows. 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 ﬁnite 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 deﬁnition 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 stratiﬁed 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 ﬂow 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 ﬂow 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 coeﬃcient assumed constant. The set of governing equations is completed by the internal energy balance, written in terms of the temperature variable, where a constant speciﬁc 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 ﬂuid, respectively. 2.2. LES formulation The above mathematical model is reformulated in the framework of LES approach based on the FV numerical methodology. The diﬀerential 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 ﬁltering, the corresponding ﬁltered variable being deﬁned 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 ﬁlter 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 ﬁlter width. When applying the ﬁltering 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 ﬁeld 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 ﬁltered velocity is imposed to be numerically divergence-free, that is the second order approximation q ﬃ 0 holds (Sagaut, 2002). As to the ﬁltered 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 ﬁltered 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 ﬂux SGS term is deﬁned according to k r# r# þ v# v#: ð10Þ t¼ q0 c p The SGS terms diﬀer from the corresponding terms in the ﬁnite diﬀerence (FD) approach due to the appearance of the viscous and thermal ﬂuxes 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 ﬃ 2mt r v; and the new pressure variable is deﬁned as p p Tr(T)/3. Similarly, the SGS heat ﬂux vector-ﬁeld t is approximated as a Fourier-like ﬂux kt tﬃ 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 deﬁned the total diﬀusion coeﬃcients 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 coeﬃcients, mt and kt, and the balanced ﬁltered ﬁelds, v and #, must be speciﬁed. Typically, the smallest ﬂow 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 ﬂow conditions, one usually considers the standard Smagorinsky eddy viscosity model, mt ¼ CD2 2S ij S ij . For buoyancy driven ﬂows, the temperature ﬁeld being not simply a passive scalar but actively entering the ﬂow dynamics, it is more complex to deﬁne a suitable eddy viscosity. In fact, the turbulence theory for thermally stratiﬁed ﬂows is still less developed than for shear ﬂows and there are no deﬁnitive research results. In particular, according to Lilly (1962), it appears doubtful that one can express the energy transfer only in terms of the mean ﬂow, regardless of the buoyancy scales. Actually, the turbulent kinetic energy budget is altered by the buoyancy eﬀects and, even considering still valid the equilibrium hypothesis, the eddy viscosity deﬁnition must be modiﬁed according to 1=2 o# mt ¼ CD2 2S ij S ij a : ð16Þ oz The coeﬃcient a in the above eddy-viscosity model measures the importance of the coupling with the temperature stratiﬁcation as a function of the Richardson number (Lilly, 1962). Some authors adopted the so-called stratiﬁcation formulation, e.g., (Sullivan and Moeng, 1992; Cabot, 1992), by directly expressing the eddy viscosity mt, along with the eddy diﬀusivity 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 coeﬃcients. In order to overcome the problem, a diﬀerent 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-ﬁlter width D and the dissipation rate e as mt ¼ C 2=3 D4=3 e1=3 ;

ð17Þ

where C stands for the Smagorinsky model coeﬃcient. This way, the dissipation rate e is no more assumed to balance the SGS energy production due to the buoyancy eﬀect 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 eﬀects. 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 ﬃ 2C e D4=3 rs v;

ð18Þ

where Ce C2/3e1/3 stands for the model coeﬃcient 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 coeﬃcients 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) ﬁlter 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 ﬁltering is deﬁned 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-ﬁlter 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 ﬁltered diﬀerential 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 ﬁlter 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 deﬁnition (8), the residual stress tensor with respect to the test ﬁlter is R ¼ 2m rs v rs ev þ vv ev ev:

ð21Þ

ð22Þ

The deﬁnitions (8) and (22) allow us to generalize the Germano identity as R T ¼ L; with the following deﬁnition 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 ﬃ 2m0t rs~v; m0t being the eddy viscosity related to the test ﬁltering procedure. This way, the identity (23) becomes 2mt rs v 2m0t rs ev ﬃ 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 ﬁltering and, congruently, the deﬁnition of the known Leonard tensor (24), apart from the presence of the molecular diﬀusion terms, is rather diﬀerent. This speciﬁc consequence of the FV integral formulation appears to be very appealing as it allows us to avoid the usual approximation with which the model coeﬃcient is extracted from the ﬁltering operation. Starting from Eq. (25), the modelling coeﬃcient 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 ﬁlters, that is m ¼ C D4=3 and m0 ¼ C D t

e

t

e

e

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

ð26Þ

Eq. (26) stands for a system of ﬁve 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 diﬀusivity kt. In fact, the integral temperature Eq. (9) is rewritten in the following grid ﬁltered diﬀerential form: o# k þ r v# ¼ r r# r t; ð28Þ ot q0 c p and, when applying the test ﬁlter, 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 ﬂux corresponding to test ﬁltering 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 ﬁeld e k r# r # e l ¼ v# ev # q0 c p is deﬁned. k0 ~ with Similarly to (12), let us consider the approximation r ﬃ 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 ﬃl C d r# /4=3 r #

ð34Þ

C e 4=3 where C d Pr D stands for a new modelling coeﬃcient. 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 coeﬃcient 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 ﬁeld by q0U2. The non-dimensional temperature variable is deﬁned as ð# #0 Þ=D#; D# being the temperature diﬀerence between bottom plane and interface corresponding to the initial stratiﬁcation. 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 coeﬃcients, (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 ﬁnite-depth stably stratiﬁed horizontal layer is considered. This ﬂow can be seen as a simpliﬁed 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 ﬁxed in order for a DNS solution to be feasible. A real-world set of non-dimensional numbers will be used in the ﬂow 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 diﬀerential 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 ﬁeld and a weak initial vertical stratiﬁcation. A negative heat ﬂux is prescribed at the air–sea interface, providing a loss of sea internal energy and the onset of motion. Therefore, the ﬂow motion is generated solely by inhomogeneities in the temperature ﬁeld. The actual ﬂow dynamics can be considered in the framework of the Rayleigh–Be´nard convective regime. Namely, after a ﬁrst stage dominated by molecular temperature diﬀusion, that causes the development of a cool thermal boundary layer, the denser ﬂuid layer becomes unstable and convective plumes are generated (Zikanov et al., 2002; Leighton et al., 2003). For instance, this ﬂow 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 ﬂux. 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 ﬁxed as U (B0H)1/3, B0 being the buoyancy ﬂux 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 ﬁxed heat ﬂux at the interface is such that Nu = 600, while the ﬂuid Prandtl number is assumed to be Pr = 1. Such ﬂow 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 veriﬁcation of the proposed numerical method, rather than in the prediction of realistic ocean ﬂows. 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 ﬁeld 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 ﬂow, 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 ﬁrst stage dominated by the thermal diﬀusion, that is followed by the onset of temperature plumes. That corresponds to the potential energy falling down, that is not reported among the ﬁgures 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 ﬁeld, are slightly diﬀerent and the numerical errors act diﬀerently. Thus, it is not surprising that the initial ﬂow instability starts at diﬀerent 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 eﬀectively 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 coeﬃcients in the other one. However, owing to the diﬀerent initial transition, in absence of a statistical steady state, there is some diﬃculty in individuating a speciﬁc 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 ﬁeld, 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 ﬁeld 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 ﬁgure 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 diﬀerent initialization. Anyway, the comparison can be considered fully satisfactory and the FV-based numerical procedure successfully validated. In order to visualize the ﬂow evolution, in Fig. 5 the temperature iso-surfaces at two diﬀerent values h = 3 (cold ﬂuid zones in blue) and h = 0.5 (hot ﬂuid 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 ﬁx a suitable LES ﬁlter 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-ﬁlter ratio /, is ﬁxed to 2, the discrete test-ﬁlter being constructed by averaging in the horizontal plane over a 3 3 FV-stencil. Note that the SGS viscosity and diﬀusivity coeﬃcients are clipped so that they stay non-negative. In order to highlight the eﬀectiveness 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 stratiﬁcation 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 ﬁltering is applied, periodically in time, in order to obtain a numerically stable solution. For a discussion of why and how this explicit ﬁltering 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 ﬁltering eﬀect of the numerical FV discretization. The temporal evolution of the volume-averaged kinetic energy is shown in Fig. 6 for the diﬀerent LES solutions. Apart a time delay for the onset of ﬂow 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 eﬃcient 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 ﬁgures, as a natural reference. Note that, for a better comparison, the DNS spectra should have been a-posteriori ﬁltered 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-ﬁlter corresponding to the LES solution. The eﬀect of the implicit smooth FV-ﬁlter is evident when considering the diminished energy content with respect to the DNS solution. The diﬀerence in the maximum resolved wavenumbers for the diﬀerent LES solutions trivially reﬂects the diﬀerent 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 ﬁeld. Also, as to comparison between FV-based dynamic and static Smagorinsky models, it appears evident that the prescribed model coeﬃcients for the static version are not well suited to the actual transitional ﬂow 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 eﬀect of the additional explicit low-pass ﬁltering on the spectral solution appears evident at highest resolved wavenumbers. As pointed out by the authors, this explicit ﬁlter must be considered as an eﬀective 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 ﬁlter 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 ﬂow 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 ﬁeld is approached. However, such issues will be further explored in the future. 5. An example of application to geophysical ﬂows: 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 ﬂows the characteristic numbers are several magnitude orders greater, a DNS calculation is unfeasible because of the unaﬀordable computational cost. Therefore, after having used the reduced scale model for validation purposes, a ﬂow problem more representative of the oceanic turbulence is now addressed. The present FV-based LES model is applied to the simulation of a stratiﬁed 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 stratiﬁed quiescent ﬂuid 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 stratiﬁcation N 20 ¼ gb o#ozh and the superﬁcial 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 ﬂow 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 ﬂow conﬁguration 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 ﬂows, i.e. for time-scales shorter than the Coriolis parameter f1 while, for larger times, the rotational eﬀects become signiﬁcant and the deepening process tends to arrest when a stratiﬁed 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 ﬁrst few hours. The prediction of the OML deepening is an important requirement for an ocean turbulence model, the OML depth deﬁnition in (41) being however not obvious. Although diﬀerent deﬁnitions exist in the literature the OML is generally deﬁned as the region of uniform ﬂuid properties just below the ocean surface. This definition implies the examination of tracer values, e.g., the ﬂuid density. For instance, one can deﬁne 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 deﬁnition consists in identifying the lower boundary of the region of high-turbulence intensity (turbocline). However, for shear turbulence in stratiﬁed ﬂuids 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 ﬂow instability and the formation of turbulence and the OML depth can be identiﬁed 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 deﬁne Ri = 1 as the critical Richardson number. Diﬀerent FV-based dynamic LES solutions of the Kato–Phillips test-case have been obtained for diﬀerent numerical parameters. Namely, the numerical solution has been analysed for diﬀerent grid-to-ﬁlter 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 reﬁnement is commonly used for ﬂows of engineering interest. A further comparison is made with a static version of the Smagorinsky model, that is by a-priori prescribing the model coeﬃcients as C = 0.1 and Prt = 0.4. All the solutions have been obtained for the ﬁrst 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 diﬀerent criteria (case FGR = 4, uniform grid). Upper left panel: vertical proﬁle of Ri. Upper right panel: vertical proﬁle 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 proﬁle of Ri, computed upon the horizontally averaged ﬁelds, 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 proﬁle 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 proﬁles after 2 hours are shown in Fig. 12. The LES results are reported for diﬀerent solutions, on uniform and non-uniform grids with ﬁlter-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 diﬀerent LES solutions, on uniform and non-uniform grids with ﬁlter-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 proﬁle for diﬀerent LES solutions, on uniform and non-uniform grids with ﬁlter-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 ﬂux and shear-stress. However, the ﬁlter-to-grid ratio plays a very important role. For a given numerical resolution, that is for a ﬁxed 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 aﬀects 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 stratiﬁed, 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 diﬃcult, if not impossible, to a-priori prescribe the model coeﬃcient in such a way that the SGS model well behaves for the diﬀerent stages of the ﬂow evolution. In fact, the dynamic model allows us to follow the onset of the physical ﬂow 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 aﬀected by non-physical energy scaling that shows an increasing of the energy content close to the cut-oﬀ. The method is further applied to the simulation of a well-known sheared and stably-stratiﬁed 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 ﬂows.

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 ﬂows, 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 proﬁle data and a proﬁle-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 ﬁltering for control volume ﬂow 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 stratiﬁed ﬂuid. J. Fluid Mech. 37, 643–655. Kravchenko, A.G., Moin, P., 1997. On the eﬀect of numerical errors in large eddy simulations of turbulent ﬂows. 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 ﬁnite diﬀerence 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. Inﬂuence 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 eﬀects 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 ﬂows. 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 ﬂows. In Proc. 10th Symp. Turbulence and Diﬀusion. 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.

Copyright © 2020 KUNDOC.COM. All rights reserved.