- Email: [email protected]
Energy-conserving discretisation of turbulent shear and buoyancy production
Ocean Modelling 4 (2002) 347–361 www.elsevier.com/locate/omodol
Energy-conserving discretisation of turbulent shear and buoyancy production Hans Burchard
*
Institut f€ur Meereskunde, Universit€at Hamburg, Troplowitzstraße 7, D-22529 Hamburg, Germany
Abstract An energy-conserving discretisation of the shear and buoyancy production terms for turbulent kinetic energy is derived. In contrast to ‘ad hoc’ numerical schemes, this guarantees that all mean kinetic and potential energy which is lost or gained due to vertical mixing is exactly subtracted or added to the turbulent kinetic energy budget. It is further shown in a numerical wind-entrainment experiment that this new methods results in significantly more stable numerical solutions. Ó 2002 Published by Elsevier Science Ltd.
1. Introduction The physical equations on which ocean models are based are derived from basic conservation principles, see e.g. Haidvogel and Beckmann (1999). In order to obtain physically meaningful model simulations by means of numerical models, the most important conservation properties of the physical equations need to be contained in the discrete solutions as well. This is straightforward for the conservation of mass, momentum, salt, heat and dynamically passive tracers if the so-called flux-formulations of the relevant advection–diffusion equations are discretised with standard finite-volume methods. Other conservation laws such as kinetic energy or entrophy conservation do require the construction of more sophisticated schemes such as those discussed by Arakawa and Lamb (1977). A property which should at least be approximately conserved in large scale models is potential vorticity, see M€ uller (1995). One conservation principle which has so far not been considered for the numerical discretisation is the conservative energy-exchange between mean kinetic energy, potential energy and
*
Now at: Baltic Sea Research Institute Warnem€ unde, Seestraße 15, D-18119 Rostock-Warnem€ unde, Germany. Tel.: +49-381-5197111; fax: +49-381-5197440. E-mail address: [email protected] (H. Burchard). 1463-5003/02/$ - see front matter Ó 2002 Published by Elsevier Science Ltd. PII: S 1 4 6 3 - 5 0 0 3 ( 0 2 ) 0 0 0 0 9 - 4
348
H. Burchard / Ocean Modelling 4 (2002) 347–361
turbulent kinetic energy. This has maybe not been thought of as a major problem, since the production of turbulent kinetic energy is approximately balanced by the dissipation into heat, which is thus a sink for mechanic energy. In this paper, an energy-conserving discretisation method for the shear and buoyancy production terms will be developed, which guarantees that all mean kinetic and potential energy which is lost or gained due to vertical mixing is exactly subtracted or added to the turbulent kinetic energy budget. It will be shown for a simple windentrainment experiment, how this new method considerably stabilises the numerical scheme in comparison to a ‘straight-forward’ scheme.
2. Physical equations For this study, a simplified set of one-dimensional dynamic equations for momentum u, buoyancy b and turbulent kinetic energy k is considered for a infinitely deep water column. Thus, horizontal transports, pressure gradients, and rotation are neglected. In this context, the momentum equation reads as follows: ot u oz ðmt oz uÞ ¼ 0
ð1Þ
with the surface boundary condition mt oz u ¼ ss ¼ u ju j;
z ¼ 0;
ð2Þ
where mt is the eddy viscosity and u the surface friction velocity. The latter is usually depending on the wind speed and direction and can thus be positive, zero or negative. The vertical coordinate z is pointing upwards with z ¼ 0 at the surface. Time is denoted by t. By multiplying (1) with u, a dynamic equation for the mean kinetic energy ekin ¼ u2 =2 is obtained: 2
ot ekin oz ðmt oz ekin Þ ¼ mt ðoz uÞ ¼ P ;
ð3Þ
where P is a sink term denoting the shear production of turbulent kinetic energy. The flux of mean kinetic energy through the surface is then mt oz ekin ¼ uu ju j;
z ¼ 0:
The buoyancy b, which is defined as q q0 b ¼ g q0
ð4Þ
ð5Þ
with gravitational acceleration g ¼ 9:81 m s2 , density q and reference density q0 (which is here set to q0 ¼ 1027 kg m3 ), is calculated by means of the following dynamic equation: ot b oz m0t oz b ¼ 0 ð6Þ with eddy viscosity m0t . Here, no flux of buoyancy through the surface is considered m0t oz b ¼ 0;
z ¼ 0:
ð7Þ
By multiplying (6) with the coordinate z, a dynamic equation for the potential energy epot ¼ zb is obtained:
H. Burchard / Ocean Modelling 4 (2002) 347–361
ot epot oz ð zÞm0t oz b ¼ m0t oz b ¼ B:
349
ð8Þ
In (8), the so-called buoyancy production term B is a source term for oz b > 0 and a sink term for oz b < 0. It should be noted that N 2 ¼ oz b is the Brunt–V€ais€al€a frequency squared. Thus, potential energy is increased for stable stratification and decreased for unstable stratification. The second term on the left hand side of Eq. (8) is a turbulent transport term with a no-flux condition at the surface (z ¼ 0), and thus vertically redistributes potential energy. The sinks and sources for mean kinetic energy ekin and potential energy epot , P and B, respectively, are dissipating or producing turbulent kinetic energy k: ot k oz ðmt oz k Þ ¼ P þ B e;
ð9Þ
see e.g. Burchard and Bolding (2001). Here, e > 0 denotes the dissipation rate of turbulent kinetic energy into heat. According to the law of the wall, a zero flux of turbulent kinetic energy through the surface is assumed: mt oz k ¼ 0;
z ¼ 0:
ð10Þ
The total energy budget E of the water column with Z 0 E¼ ðekin þ epot þ kÞdz
ð11Þ
1
is then ot E ¼ uð0Þu ju j
Z
0
e dz:
ð12Þ
1
Thus, total energy is injected or extracted at the surface (depending on the sign of uð0Þu ) and dissipated into heat throughout the water column.
3. Discretisation The equations for momentum and buoyancy will be discretised on a grid in time and space, see Fig. 1. In order to obtain a finite number of grid points, a lower bound z ¼ H for the vertical coordinate has to be chosen such that the solution near the surface is not affected. The onedimensional space is divided into M intervals hj ¼ zjþ1=2 zj1=2 with H ¼ z1=2 < z3=2 < < zM1=2 < zMþ1=2 ¼ 0. The discrete values for momentum, uj , and buoyancy, bj , are located in the centres of the intervals at zj ¼ 12ðzjþ1=2 þ zj1=2 Þ, the turbulent quantities for eddy viscosity and eddy diffusivity (for convenience both denoted by mjþ1=2 ), turbulent kinetic energy, kjþ1=2 , and turbulent dissipation rate, ejþ1=2 , are located on the interval interfaces at zjþ1=2 , see Fig. 1. The constant time increment is denoted by Dt. For a fixed time t, uj denotes the velocity at time t, u^j the velocity at time t þ Dt and uj the velocity at time t þ rDt with 0 6 r 6 1, see Fig. 1. After these preliminaries, the idealised momentum equation (1) can be readily discretised as follows:
350
H. Burchard / Ocean Modelling 4 (2002) 347–361
Fig. 1. Spatial (upper panel) and temporal (lower panel) organisation and indexing of the numerical grid. Here, a time stepping slightly more implicit than the Crank–Nicholson scheme with r ¼ 0:6 is shown.
uM uM1 u^M uM ss mM1=2 zM zM1 ¼ 0; Dt zMþ1=2 zM1=2 ujþ1 uj
uj uj1
u^j uj mjþ1=2 zjþ1 zj mj1=2 zj zj1 ¼ 0; Dt zjþ1=2 zj1=2
for j ¼ 2; . . . ; M 1
ð13Þ
u u
j 2 u^1 u1 m3=2 z2 z1 sb ¼0 Dt z3=2 z1=2
with the surface stress ss and the bed stress sb (which is set to zero here). For r > 0, this discrete equation for u^j is implicit, such that a linear system with a tri-diagonal matrix has to be solved. In analogy to the derivation of the mean kinetic energy equation (3), a discrete mean kinetic uj þ uj Þ. After some algebraic operenergy equation can be derived by multiplying (13) with 12ð^ ations, the following equation results for the discrete mean kinetic energy e^j ¼ 12 u^2j :
H. Burchard / Ocean Modelling 4 (2002) 347–361 ejþ1 ej
351
ej ej1
e^j ej mjþ1=2 zjþ1 zj mj1=2 zj zj1 1 ¼ mjþ1=2 2 Dt zjþ1=2 zj1=2 rð^ ujþ1 u^j Þð^ ujþ1 uj Þ þ ð1 rÞðujþ1 uj Þðujþ1 u^j Þ ðzjþ1=2 zj1=2 Þðzjþ1 zj Þ 1 mj1=2 2 rð^ uj u^j1 Þð^ uj uj1 Þ þ ð1 rÞðuj uj1 Þð^ uj uj1 Þ ðzjþ1=2 zj1=2 Þðzj zj1 Þ l u Pj1=2 ; ¼: Pjþ1=2
for j ¼ 2; . . . ; M 1:
ð14Þ
The equations for the surface and the bottom layer are omitted here. In the surface layer, the first term in the nominator on the left hand side has to be replaced by the surface energy flux, 1 ð^ u þ uM Þss and in the bottom layer, the second term in this nominator will be 12 ð^ u1 þ u1 Þsb . 2 M l u ¼ P1=2 ¼ 0. The terms on the right hand side are discrete sinks of discrete Furthermore, PMþ1=2 mean kinetic energy. They have to be source terms for the discrete turbulent kinetic energy equation: k
k
k
k
jþ3=2 jþ1=2 jþ1=2 j1=2 k^jþ1=2 kjþ1=2 mjþ1 zjþ3=2 zjþ1=2 mj zjþ1=2 zj1=2 ¼ Pjþ1=2 þ Bjþ1=2 ejþ1=2 ; Dt zjþ1 zj
for j ¼ 2; . . . ; M 1: ð15Þ
For the surface and the bottom layers, Eq. (15) would have to be modified such that fluxes of TKE at z ¼ zM and z ¼ z1 are given. For the shear dominated problem investigated here, these fluxes are set to zero in accordance with the law of the wall. Due to the staggering of the grid with the discrete values of momentum being located between the discrete values of turbulent kinetic energy, the discrete shear production terms must be of the following form in order to guarantee conservation of total energy: Pjþ1=2 ¼
l u Pjþ1=2 ðzjþ1=2 zj1=2 Þ þ Pjþ1=2 ðzjþ3=2 zjþ1=2 Þ
¼ mjþ1=2
zjþ1 zj ð ujþ1 uj Þð~ ujþ1 u~j Þ ðzjþ1 zj Þ2
;
for j ¼ 2; . . . ; M 1;
ð16Þ
where u~j ¼ 12 ð^ uj þ uj Þ. For deriving (16), the discrete Eq. (14) was first multiplied with ðzjþ1=2 zj1=2 Þ in order to obtain layer-integrated discrete energy fluxes. When inserting the discrete production terms Pju and Pjl into the discrete turbulent kinetic energy equation, they have to be divided by the thickness of the reference layers for turbulent kinetic energy. Such a weighted averaging is specifically important in three-dimensional models where the TKE equation and the momentum equations are usually located at different horizontal grid points and thus at different water depths. For terrain-following coordinates, this has the consequence that the layer thicknesses for the momentum equations and for the TKE equation can be substantially different, even when vertically equidistant r coordinates are used. However, using also the weighting in time
352
H. Burchard / Ocean Modelling 4 (2002) 347–361
included in (16) for three-dimensional models would require to additionally store the old velocity fields, which is usually not needed in two-level models. This formulation differs from the form which would be obtained by directly discretising the definition of the turbulent shear production given in (3). Such a direct ‘ad hoc’ discretisation would be for example: " #2 u^jþ1 u^j : ð17Þ Pjþ1=2 ¼ mjþ1=2 zjþ1 zj It can easily seen that (16) and (17) are identical only for the steady-state (^ uj ¼ uj ) case. Another ‘ad hoc’ discretisation could be the semi-implicit calculation for Pjþ1=2 : " #2 u~jþ1 u~j Pjþ1=2 ¼ mjþ1=2 ð18Þ zjþ1 zj which is identical to (16) for r ¼ 0:5. This discretisation is often applied when a leap-frog timestepping is used, see e.g. the Princeton Ocean Model (POM; Blumberg and Mellor, 1987). This semi-implicit choice for r is known to cause numerical instabilities, see e.g. Patankar (1980). Investigating this further here, is however beyond the scope of this paper. An equivalent formalism can be applied to the buoyancy equation (6). With the formalism given above, the discretisation is of the following form: b
b
b b
j j jþ1 j1 b^j bj mjþ1=2 zjþ1 zj mj1=2 zj zj1 ¼ 0: Dt zjþ1=2 zj1=2
ð19Þ
After multiplication with zj and some algebraic operations, the following equation for the discrete potential energy pj ¼ zj bj is obtained: b
b
b b
j j jþ1 j1 1 1 p^j pj 2 ðzjþ1 þ zj Þmjþ1=2 zjþ1 zj 2 ðzj þ zj1 Þmj1=2 zj zj1 þ Dt zjþ1=2 zj1=2 1 1 bjþ1 bj bj bj1 ¼ mjþ1=2 þ mj1=2 ¼: Bljþ1=2 Buj1=2 : 2 zjþ1=2 zj1=2 2 zjþ1=2 zj1=2
ð20Þ
The terms on the right hand side are the discrete expressions for the sinks and sources of potential energy due to vertical mixing of density. In order to achieve total energy conservation, these have to be used for the discretisation of the turbulent buoyancy production. The discrete buoyancy production for kjþ1=2 would thus be of the following form: bj bjþ1 ð21Þ Bjþ1=2 ¼ mjþ1=2 zjþ1 zj which again differs from a direct discretisation which would be formulated as: Bjþ1=2 ¼ mjþ1=2
b^jþ1 b^j : zjþ1 zj
ð22Þ
Formulations (21) and (22) are identical only for steady-state flow or fully implicit time stepping with r ¼ 1 and, consequently, bj ¼ b^j .
H. Burchard / Ocean Modelling 4 (2002) 347–361
353
Since the production terms involve the new velocity u^j and the new buoyancy b^j , the velocity and buoyancy equation are solved before the turbulent kinetic energy equation in each time step. Thus, the old eddy viscosity and diffusivity value have to be used in the velocity and buoyancy equation. Due to the completely independent derivation of the energy-conserving shear and buoyancy production terms, different values of r can of course be used for momentum and buoyancy. For temperature and salinity, which are usually calculated first as the two quantities defining buoyancy, the same r-values should be used. The physical model has to be completed by formulations for eddy viscosity mt , eddy diffusivity m0t and turbulent dissipation rate, e. Here, a model as simple as possible has been chosen in order to reduce the complexity of the calculation, but still retain basic physical properties. Eddy viscosity and diffusivity are calculated as mt ¼ cl
k2 ; e
m0t ¼
mt Pr
ð23Þ
with the empirical constant cl ¼ 0:09, see Rodi (1980) and the turbulent Prandtl number Ri Ri Pr ¼ Pr0 exp 0 1 þ 1 Pr Ri Ri
ð24Þ
2
with the gradient Richardson number Ri ¼ oz b=ðoz uÞ , the neutral turbulent Prandtl number Pr0 ¼ 0:74 and Ri1 ¼ 0:25, see Schumann and Gerz (1995). The dissipation rate is calculated as follows: e ¼ c3=4 l
k 3=2 L
ð25Þ
with the macro length scale L¼
1 1 þ k 1=2 jðz þ z0 Þ c N
!1 ;
ð26Þ
where the first term is a length limitation near the surface given by the law of the wall and the second term is a length limitation due to the effect of stable stratification. The surface roughness parameter is here set to z0 ¼ 0:1 m and the von Karman constant is given by j ¼ 0:4. The stability parameter c should depend on the stability functions cl and c0l and has been calibrated here to c ¼ 0:4 such that it reasonably reproduces the mixed layer depth evolution for the windentrainment experiment described below, see Eq. (30). This length scale parameterisation is similar to that of Axell and Liungman (2001) who used slightly different values for c. In order to guarantee positivity for the turbulent kinetic energy, the sink terms on the right hand side of the discrete k-equation have to be treated quasi-implicitly, see Patankar (1980). This leads to the following discrete form of the eddy diffusivity in the term Bjþ1=2 from (21): 8 k2 < cl jþ1=2 for bjþ1 < bj ; ejþ1=2 mjþ1=2 ¼ ð27Þ ^jþ1=2 kjþ1=2 k :c else: l
ejþ1=2
354
H. Burchard / Ocean Modelling 4 (2002) 347–361
An equivalent discretisation is applied when (22) is used for the buoyancy production term. Thus, the energy flux from the turbulent kinetic energy to the potential energy is not conservative for stable stratification. There is one other reason why the budget of the turbulent kinetic energy is not closed: In order to avoid too low values for the turbulent kinetic energy (which would lead to an underflow in numerical calculations), a lower limiting value is imposed, which is here: k > kmin ¼ 1010 J m2 :
ð28Þ
This has the consequence of a small artificial positive energy flux in areas with net turbulence dissipation. With (25), the dissipation term on the right hand side of the turbulent kinetic energy equation will be discretised according to Patankar (1980) as follows: ejþ1=2 ¼
c3=4 l
1=2 k^jþ1=2 kjþ1=2
Ljþ1=2
:
ð29Þ
This does however not violate the total energy conservation, since the turbulent dissipation rate is a sink of the total energy, see Eq. (12). The wind-mixing experiment used here is motivated by the laboratory experiment described by Kato and Phillips (1969). There, a mixed layer induced by a constant surface stress penetrates into a stably stratified fluid with density increasing linearly down from the surface. The water depth is assumed to be infinite. Price (1979) suggested a solution for the evolution of the mixed-layer depth Dm based on a constant Richardson number 1=2 1=2
Dm ðtÞ ¼ 1:05u N0
t
;
ð30Þ
where u is the surface friction velocity and N0 the constant initial Brunt–V€ais€al€a frequency. Following several authors, see e.g. Deleersnijder and Luyten (1994); Burchard et al. (1998) this laboratory experiment is transformed to ocean dimensions with u ¼ 102 m s1 and N0 ¼ 102 s1 . The empirical estimate (30) results in a mixed layer depth of Dm ¼ 34:5 m after 30 h of entrainment, which will be the run time of this numerical experiment. A total depth of H ¼ 50 m is chosen here which is sufficiently deep for simulating an infinite depth for this study.
4. Numerical performance study In this section, the numerical performance of the energy-conserving formulations for turbulent shear production and buoyancy production, (16) and (21) is investigated by means of the simple wind-mixing experiment given above and compared to the non-conservative schemes (17) and (22), respectively. The performance of the relatively simple turbulence closure model consisting of Eqs. (9), (23), (24), (25), and (26) is first compared to a more sophisticated closure model consisting of an additional transport equation for the dissipation rate e instead of using Eqs. (25) and (26). The stability functions in that model are based on an algebraic second-moment closure recently suggested by Canuto et al. (2001). The complete model has been described and tested in detail by Burchard and Bolding (2001). Fig. 2 shows profiles of eddy viscosity mt , turbulent kinetic energy k,
H. Burchard / Ocean Modelling 4 (2002) 347–361
355
Fig. 2. Profiles of eddy viscosity mt (left), turbulent kinetic energy k (middle), and macro length scale L (right) after 30 h of wind-entrainment. Calculations have been carried out with the one equation k model with algebraic length scale (––), see Eq. (26) and with the two-equation k–e model with algebraic second-moment closure by Canuto et al. (2001) (- - -).
and macro length scale L after 30 h of wind-entrainment calculated with both models. Despite the structural differences between the two turbulence closure models, the resulting turbulence profiles do not differ that much. This similarity to the well-tested, more complex model, gives some credibility to the simpler model with algebraic length scale. It is thus concluded that this model is sufficiently general in order to represent turbulence closure models for this numerical performance test. Three different resolutions in time and space have been used for this study: 1. A non-equidistant grid spacing with strong zooming towards the surface: tanh ðdl þ du Þ Mj dl þ tanhðdl Þ 1 ; j ¼ 0; . . . ; M zjþ1=2 ¼ H tanhðdl Þ þ tanhðdu Þ
ð31Þ
with the zooming coefficients du ¼ 3, dl ¼ 0 and the number of layers, M ¼ 200. This results in a near-surface layer depth of hM ¼ 0:008 m and a near bed layer depth of h1 ¼ 0:754 m. The time step chosen here is Dt ¼ 200 s. 2. An equidistant grid-spacing with M ¼ 200 (resulting in layer thicknesses of 0.25 m) and Dt ¼ 200 s. 3. An equidistant grid-spacing with M ¼ 200 (resulting in layer thicknesses of 0.25 m) and Dt ¼ 20 s. The level of implicitness has been set to r ¼ 0:6 for the momentum and the buoyancy equation, a value slightly larger than the semi-implicit value 0.5 ensuring second-order accuracy in time, but which is also known for causing numerical instabilities. For the turbulent kinetic energy equation, a fully implicit time-stepping with r ¼ 1 has been chosen in order to further stabilise the model. Fig. 3 shows profiles of eddy viscosity mt obtained with the new energy-conserving discretisation of turbulence production terms (16) and (21) and the old not energy-conserving discretisation (17) and (22), using the three spatial and temporal resolutions 1–3 given above.
356
H. Burchard / Ocean Modelling 4 (2002) 347–361
Fig. 3. Profiles of eddy viscosity mt obtained by applying the new ( ) and the old ( ) discretisations of turbulence production terms. The three spatial and temporal discretisations 1 (left), 2 (middle) and 3 (right) have been used.
It is clearly seen that the old discretisation is numerically unstable in the region of maximum eddy viscosity for non-equidistant grid-spacing with long time steps of Dt ¼ 200 s. In contrast to this, the new discretisation results in a stable and smooth solution which is not significantly deviating (by optical inspection) from the equidistant grid-spacing case with short of Dt ¼ 20 s, see Fig. 2 (left panel) and Fig. 3 (right panel). For equidistant grid-spacing and large time steps of Dt ¼ 200 s, the old method is numerically stable at the end of the run time, but significantly underestimates the mixed layer depth predicted obtained with high resolution. Only for equidistant grid spacing and small time steps of Dt ¼ 20 s, the old method results in a final viscosity profile optically not distinguishable from the profile obtained with the new method. It is also instructing to inspect some aspects of the temporal evolution of the numerical solutions. In Fig. 4, time series of the near surface velocity u at z ¼ 0:125 m are given. This coincides with the discrete velocities at M ¼ 200 for equidistant and M ¼ 187 for non-equidistant gridspacing given by (31) with du ¼ 3. In Fig. 4, the time axis is given in logarithmic scale in order to emphasise the instabilities occurring directly after the sudden onset of surface forcing over a laminar water column. These instabilities are visible for all model runs with long time steps of Dt ¼ 200 s. However, the near-surface velocity converges towards the result for Dt ¼ 20 s after 4 h for the non-equidistant case and after 1 hour for the equidistant case, when the new method is used. With the old method, some approximation of the high temporal resolution solution after about 6 h is visible, but even at the end of the run time, visible differences are present, see also Fig. 3. Finally, balances for the vertically integrated energy equations are plotted in Fig. 5 (new discretisation) and 6 (old discretisation) for non-equidistant grid-spacing and Dt ¼ 200 s. The total energy budget, Eq. (12), the mean kinetic energy budget, vertical integral of Eq. (3), the potential energy budget, vertical integral of Eq. (8), and the turbulent kinetic energy budget, vertical integral of Eq. (9) are shown. For the new discretisation, the total energy balance is basically an equilibrium between energy input due to surface stress, increase of mean kinetic energy due to mixing and turbulent dissipation, similarly to the mean kinetic energy budget where the energy input by wind is balanced by an increase of mean kinetic energy and shear production of turbulent kinetic energy. The turbulent kinetic energy budget is mainly a balance between shear production
H. Burchard / Ocean Modelling 4 (2002) 347–361
357
Fig. 4. Time series of near surface velocity u at z ¼ 0:125 m computed with the three spatial and temporal discretisations given above. Upper panel: new discretisation of turbulence production terms (16) and (21); lower panel: old discretisation of turbulence production terms (17) and (22). For time axis is given in logarithmic scale for clarity.
and dissipation, slightly modified by loss of turbulent kinetic energy to potential energy, for which the buoyancy production term fully balances the increase of potential energy. Some numerical instabilities are visible during the first three hours of the simulation for the shear production P, the dissipation rate e and the tendency of turbulent kinetic energy, ot k. The residuals for the mean kinetic and the potential energy are exactly zero (down to machine accuracy), but the turbulent kinetic energy budget and the total energy budget are not fully closed due to the two model features discussed above, the quasi-implicit treatment of sink terms due to Patankar (1980) and the lower limit set for the turbulent kinetic energy. The vertically integrated energy budgets for the old discretisation are shown in Fig. 6. It can be seen that the shear production is some orders of magnitude larger than the energy input due to surface wind stress. This is merely due to the non-conservative discretisation of the shear production term. This leads to an excessive production of turbulent kinetic energy, which is not fully compensated by dissipation, which depends on k and N, two variables which are lagging behind the shear production term. It seems that a positive feedback occurs in the mean kinetic energy budget during the first three hours of the simulation. Too large shear production causes too large turbulent kinetic energy, which in turn causes too large eddy viscosity. This again increases the
358
H. Burchard / Ocean Modelling 4 (2002) 347–361
Fig. 5. Results of the wind-entrainment experiment using the new, energy-conserving discretisations for the turbulence production terms, (16) and (21). Balance terms in the vertically integrated energy budget equations for total mechanic energy (upper left), mean kinetic energy (upper right), potential energy (lower left) and turbulent kinetic energy (lower right).
shear production. Also the potential energy budget is out of balance during the first five hours of the experiment due to the non-conservative discretisation of the buoyancy production term. It can be assumed that the unbalanced energy budgets cause the prevailing numerical instabilities seen in Figs. 3 and 4. However, a strict mathematical reasoning for this could not be found yet.
H. Burchard / Ocean Modelling 4 (2002) 347–361
359
Fig. 6. Results of the wind-entrainment experiment using the old, not energy-conserving discretisations for the turbulence production terms, (17) and (22). Balance terms in the vertically integrated energy budget equations for total mechanic energy (upper left), mean kinetic energy (upper right), potential energy (lower left) and turbulent kinetic energy (lower right).
5. Conclusions The new energy-conserving discretisations for the shear and buoyancy production terms of turbulent kinetic energy presented here differ from ‘ad hoc’ discretisations by the level of implicitness which is used for the vertical differences. While the procedure for the linear buoyancy
360
H. Burchard / Ocean Modelling 4 (2002) 347–361
production term seems to be straight-forward, the non-linear shear production term is based on a product of differences which are taken on different levels of implicitness. Although the differences between the ‘new’ and the ‘old’ formulations are small, their consequences for the overall performance of the numerical model are dramatic. After some initial instabilities due to the shock-like model initialisation, the ‘new’ discretisations compute even for strong zooming and large time steps numerically stable results for time series of near-surface velocity and also for profiles of eddy viscosity, a parameter known for its sensitivity to numerical instabilities, see Burchard and Deleersnijder (2001). In contrast to this, the initial instabilities do not vanish when the ‘old’ discretisation is used. Further experiments (not shown here) using twoequation k–e models with complex algebraic second-moment closure proved a similar stability behaviour than the simpler one-equation model discussed in this paper. It is assumed that also other turbulence models such as the k–kL two-equation model by Mellor and Yamada (1982) or the k–x model with x ¼ e=k, see Wilcox (1998) would profit from the new production term discretisations. It is satisfying that a new discretisation which ensures basic physical principles also for the numerical solution is indeed numerically more stable and thus also of practical relevance. An implementation of this new method into three-dimensional ocean models is straight-forward. It should however be noted that it is usually not the buoyancy b which is prognostically calculated in ocean models but the potential temperature and salinity from which the buoyancy is then diagnostically calculated by means of a non-linear equation of state. Thus, the energy flux between the potential energy and the turbulent kinetic energy is not completely closed, which is a consequence of the Boussinesq approximation usually made in ocean modelling.
Acknowledgements This study has been funded through a Habilitation grant by the Deutsche Forschungsgemeinschaft (German Research Foundation). Further support has been provided by the CARTUM (Comparative Analysis and Rationalisation of Second-Moment Turbulence Models) project, a concerted action of the MAST-III program of the European Commission (MAS3-CT98-0172). The author is indebted to Eric Deleersnijder (Universite Catholique de Louvain, Institut d’Astronomie et de Geophysique George Lema^ıtre, Belgium) for detailed discussions on this topic. Parts of the work presented in this paper have been carried out during a stay of Hans Burchard at this institute. The comments of two anonymous referees were greatly appreciated. For the numerical calculations in this study, of the General Ocean Turbulence Model (GOTM, http:// www.gotm.net) has been used.
References Arakawa, A., Lamb, V.R., 1977. Computational design of the basic dynamical processes of the UCLA general circulation model. Meth. Comput. Phys., 173–263. Axell, L., Liungman, O., 2001. A one-equation turbulence model for geophysical applications: comparison with date and the k–epsilon model. Environ. Fluid Mech. 1, 71–106.
H. Burchard / Ocean Modelling 4 (2002) 347–361
361
Blumberg, A.F., Mellor, G.L., 1987. A description of a coastal ocean circulation model. In: Heaps, N.S. (Ed.), Three dimensional Ocean Models. American Geophysical Union, Washington D.C, pp. 1–16. Burchard, H., Bolding, K., 2001. Comparative analysis of four second-moment turbulence closure models for the oceanic mixed layer. J. Phys. Oceanogr. 31, 1943–1968. Burchard, H., Deleersnijder, E., 2001. Stability of algebraic non-equilibrium second-order closure models. Ocean Modell. 3, 33–50. Burchard, H., Petersen, O., Rippeth, T.P., 1998. Comparing the performance of the Mellor–Yamada and the k–e twoequation turbulence models. J. Geophys. Res. 103, 10543–10554. Canuto, V.M., Howard, A., Cheng, Y., Dubovikov, M.S., 2001. Ocean turbulence I: one-point closure model. Momentum and heat vertical diffusivities. J. Phys. Oceanogr. 31, 1413–1426. Deleersnijder, E., Luyten, P., 1994. On the practical advantages of the quasi-equilibrium version of the Mellor and Yamada level 2.5 turbulence closure applied to marine modelling. Appl. Math. Modell. 18, 281–287. Haidvogel, D.B., Beckmann, A., 1999. Numerical Ocean Circulation Modelling. Series on Environmental Science and Management, vol. 2. Imperial College Press, London. Kato, H., Phillips, O.M., 1969. On the penetration of a turbulent layer into stratified fluid. J. Fluid Mech. 37, 643–655. Mellor, G.L., Yamada, T., 1982. Development of a turbulence closure model for geophysical fluid problems. Rev. Geophys. 20, 851–875. M€ uller, P., 1995. Ertel’s potential vorticity theorem in physical oceanography. Rev. Geophys. 33, 67–98. Patankar, S.V., 1980. Numerical Heat Transfer and Fluid flow. McGraw-Hill, New York. Price, J.F., 1979. On the scaling of stress-driven entrainment experiments. J. Fluid Mech. 90, 509–529. Rodi, W., 1980. Turbulence models and their application in hydraulics. Tech. Rep., Int. Assoc. for Hydraul. Res., Delft, The Netherlands. Schumann, U., Gerz, T., 1995. Turbulent mixing in stably stratified shear flows. J. Appl. Meteorol. 34, 33–48. Wilcox, D.C., 1998. Turbulence Modelling for CFD, 2nd ed. DCW Industries.
Energy-conserving discretisation of turbulent shear and buoyancy production Hans Burchard
*
Institut f€ur Meereskunde, Universit€at Hamburg, Troplowitzstraße 7, D-22529 Hamburg, Germany
Abstract An energy-conserving discretisation of the shear and buoyancy production terms for turbulent kinetic energy is derived. In contrast to ‘ad hoc’ numerical schemes, this guarantees that all mean kinetic and potential energy which is lost or gained due to vertical mixing is exactly subtracted or added to the turbulent kinetic energy budget. It is further shown in a numerical wind-entrainment experiment that this new methods results in significantly more stable numerical solutions. Ó 2002 Published by Elsevier Science Ltd.
1. Introduction The physical equations on which ocean models are based are derived from basic conservation principles, see e.g. Haidvogel and Beckmann (1999). In order to obtain physically meaningful model simulations by means of numerical models, the most important conservation properties of the physical equations need to be contained in the discrete solutions as well. This is straightforward for the conservation of mass, momentum, salt, heat and dynamically passive tracers if the so-called flux-formulations of the relevant advection–diffusion equations are discretised with standard finite-volume methods. Other conservation laws such as kinetic energy or entrophy conservation do require the construction of more sophisticated schemes such as those discussed by Arakawa and Lamb (1977). A property which should at least be approximately conserved in large scale models is potential vorticity, see M€ uller (1995). One conservation principle which has so far not been considered for the numerical discretisation is the conservative energy-exchange between mean kinetic energy, potential energy and
*
Now at: Baltic Sea Research Institute Warnem€ unde, Seestraße 15, D-18119 Rostock-Warnem€ unde, Germany. Tel.: +49-381-5197111; fax: +49-381-5197440. E-mail address: [email protected] (H. Burchard). 1463-5003/02/$ - see front matter Ó 2002 Published by Elsevier Science Ltd. PII: S 1 4 6 3 - 5 0 0 3 ( 0 2 ) 0 0 0 0 9 - 4
348
H. Burchard / Ocean Modelling 4 (2002) 347–361
turbulent kinetic energy. This has maybe not been thought of as a major problem, since the production of turbulent kinetic energy is approximately balanced by the dissipation into heat, which is thus a sink for mechanic energy. In this paper, an energy-conserving discretisation method for the shear and buoyancy production terms will be developed, which guarantees that all mean kinetic and potential energy which is lost or gained due to vertical mixing is exactly subtracted or added to the turbulent kinetic energy budget. It will be shown for a simple windentrainment experiment, how this new method considerably stabilises the numerical scheme in comparison to a ‘straight-forward’ scheme.
2. Physical equations For this study, a simplified set of one-dimensional dynamic equations for momentum u, buoyancy b and turbulent kinetic energy k is considered for a infinitely deep water column. Thus, horizontal transports, pressure gradients, and rotation are neglected. In this context, the momentum equation reads as follows: ot u oz ðmt oz uÞ ¼ 0
ð1Þ
with the surface boundary condition mt oz u ¼ ss ¼ u ju j;
z ¼ 0;
ð2Þ
where mt is the eddy viscosity and u the surface friction velocity. The latter is usually depending on the wind speed and direction and can thus be positive, zero or negative. The vertical coordinate z is pointing upwards with z ¼ 0 at the surface. Time is denoted by t. By multiplying (1) with u, a dynamic equation for the mean kinetic energy ekin ¼ u2 =2 is obtained: 2
ot ekin oz ðmt oz ekin Þ ¼ mt ðoz uÞ ¼ P ;
ð3Þ
where P is a sink term denoting the shear production of turbulent kinetic energy. The flux of mean kinetic energy through the surface is then mt oz ekin ¼ uu ju j;
z ¼ 0:
The buoyancy b, which is defined as q q0 b ¼ g q0
ð4Þ
ð5Þ
with gravitational acceleration g ¼ 9:81 m s2 , density q and reference density q0 (which is here set to q0 ¼ 1027 kg m3 ), is calculated by means of the following dynamic equation: ot b oz m0t oz b ¼ 0 ð6Þ with eddy viscosity m0t . Here, no flux of buoyancy through the surface is considered m0t oz b ¼ 0;
z ¼ 0:
ð7Þ
By multiplying (6) with the coordinate z, a dynamic equation for the potential energy epot ¼ zb is obtained:
H. Burchard / Ocean Modelling 4 (2002) 347–361
ot epot oz ð zÞm0t oz b ¼ m0t oz b ¼ B:
349
ð8Þ
In (8), the so-called buoyancy production term B is a source term for oz b > 0 and a sink term for oz b < 0. It should be noted that N 2 ¼ oz b is the Brunt–V€ais€al€a frequency squared. Thus, potential energy is increased for stable stratification and decreased for unstable stratification. The second term on the left hand side of Eq. (8) is a turbulent transport term with a no-flux condition at the surface (z ¼ 0), and thus vertically redistributes potential energy. The sinks and sources for mean kinetic energy ekin and potential energy epot , P and B, respectively, are dissipating or producing turbulent kinetic energy k: ot k oz ðmt oz k Þ ¼ P þ B e;
ð9Þ
see e.g. Burchard and Bolding (2001). Here, e > 0 denotes the dissipation rate of turbulent kinetic energy into heat. According to the law of the wall, a zero flux of turbulent kinetic energy through the surface is assumed: mt oz k ¼ 0;
z ¼ 0:
ð10Þ
The total energy budget E of the water column with Z 0 E¼ ðekin þ epot þ kÞdz
ð11Þ
1
is then ot E ¼ uð0Þu ju j
Z
0
e dz:
ð12Þ
1
Thus, total energy is injected or extracted at the surface (depending on the sign of uð0Þu ) and dissipated into heat throughout the water column.
3. Discretisation The equations for momentum and buoyancy will be discretised on a grid in time and space, see Fig. 1. In order to obtain a finite number of grid points, a lower bound z ¼ H for the vertical coordinate has to be chosen such that the solution near the surface is not affected. The onedimensional space is divided into M intervals hj ¼ zjþ1=2 zj1=2 with H ¼ z1=2 < z3=2 < < zM1=2 < zMþ1=2 ¼ 0. The discrete values for momentum, uj , and buoyancy, bj , are located in the centres of the intervals at zj ¼ 12ðzjþ1=2 þ zj1=2 Þ, the turbulent quantities for eddy viscosity and eddy diffusivity (for convenience both denoted by mjþ1=2 ), turbulent kinetic energy, kjþ1=2 , and turbulent dissipation rate, ejþ1=2 , are located on the interval interfaces at zjþ1=2 , see Fig. 1. The constant time increment is denoted by Dt. For a fixed time t, uj denotes the velocity at time t, u^j the velocity at time t þ Dt and uj the velocity at time t þ rDt with 0 6 r 6 1, see Fig. 1. After these preliminaries, the idealised momentum equation (1) can be readily discretised as follows:
350
H. Burchard / Ocean Modelling 4 (2002) 347–361
Fig. 1. Spatial (upper panel) and temporal (lower panel) organisation and indexing of the numerical grid. Here, a time stepping slightly more implicit than the Crank–Nicholson scheme with r ¼ 0:6 is shown.
uM uM1 u^M uM ss mM1=2 zM zM1 ¼ 0; Dt zMþ1=2 zM1=2 ujþ1 uj
uj uj1
u^j uj mjþ1=2 zjþ1 zj mj1=2 zj zj1 ¼ 0; Dt zjþ1=2 zj1=2
for j ¼ 2; . . . ; M 1
ð13Þ
u u
j 2 u^1 u1 m3=2 z2 z1 sb ¼0 Dt z3=2 z1=2
with the surface stress ss and the bed stress sb (which is set to zero here). For r > 0, this discrete equation for u^j is implicit, such that a linear system with a tri-diagonal matrix has to be solved. In analogy to the derivation of the mean kinetic energy equation (3), a discrete mean kinetic uj þ uj Þ. After some algebraic operenergy equation can be derived by multiplying (13) with 12ð^ ations, the following equation results for the discrete mean kinetic energy e^j ¼ 12 u^2j :
H. Burchard / Ocean Modelling 4 (2002) 347–361 ejþ1 ej
351
ej ej1
e^j ej mjþ1=2 zjþ1 zj mj1=2 zj zj1 1 ¼ mjþ1=2 2 Dt zjþ1=2 zj1=2 rð^ ujþ1 u^j Þð^ ujþ1 uj Þ þ ð1 rÞðujþ1 uj Þðujþ1 u^j Þ ðzjþ1=2 zj1=2 Þðzjþ1 zj Þ 1 mj1=2 2 rð^ uj u^j1 Þð^ uj uj1 Þ þ ð1 rÞðuj uj1 Þð^ uj uj1 Þ ðzjþ1=2 zj1=2 Þðzj zj1 Þ l u Pj1=2 ; ¼: Pjþ1=2
for j ¼ 2; . . . ; M 1:
ð14Þ
The equations for the surface and the bottom layer are omitted here. In the surface layer, the first term in the nominator on the left hand side has to be replaced by the surface energy flux, 1 ð^ u þ uM Þss and in the bottom layer, the second term in this nominator will be 12 ð^ u1 þ u1 Þsb . 2 M l u ¼ P1=2 ¼ 0. The terms on the right hand side are discrete sinks of discrete Furthermore, PMþ1=2 mean kinetic energy. They have to be source terms for the discrete turbulent kinetic energy equation: k
k
k
k
jþ3=2 jþ1=2 jþ1=2 j1=2 k^jþ1=2 kjþ1=2 mjþ1 zjþ3=2 zjþ1=2 mj zjþ1=2 zj1=2 ¼ Pjþ1=2 þ Bjþ1=2 ejþ1=2 ; Dt zjþ1 zj
for j ¼ 2; . . . ; M 1: ð15Þ
For the surface and the bottom layers, Eq. (15) would have to be modified such that fluxes of TKE at z ¼ zM and z ¼ z1 are given. For the shear dominated problem investigated here, these fluxes are set to zero in accordance with the law of the wall. Due to the staggering of the grid with the discrete values of momentum being located between the discrete values of turbulent kinetic energy, the discrete shear production terms must be of the following form in order to guarantee conservation of total energy: Pjþ1=2 ¼
l u Pjþ1=2 ðzjþ1=2 zj1=2 Þ þ Pjþ1=2 ðzjþ3=2 zjþ1=2 Þ
¼ mjþ1=2
zjþ1 zj ð ujþ1 uj Þð~ ujþ1 u~j Þ ðzjþ1 zj Þ2
;
for j ¼ 2; . . . ; M 1;
ð16Þ
where u~j ¼ 12 ð^ uj þ uj Þ. For deriving (16), the discrete Eq. (14) was first multiplied with ðzjþ1=2 zj1=2 Þ in order to obtain layer-integrated discrete energy fluxes. When inserting the discrete production terms Pju and Pjl into the discrete turbulent kinetic energy equation, they have to be divided by the thickness of the reference layers for turbulent kinetic energy. Such a weighted averaging is specifically important in three-dimensional models where the TKE equation and the momentum equations are usually located at different horizontal grid points and thus at different water depths. For terrain-following coordinates, this has the consequence that the layer thicknesses for the momentum equations and for the TKE equation can be substantially different, even when vertically equidistant r coordinates are used. However, using also the weighting in time
352
H. Burchard / Ocean Modelling 4 (2002) 347–361
included in (16) for three-dimensional models would require to additionally store the old velocity fields, which is usually not needed in two-level models. This formulation differs from the form which would be obtained by directly discretising the definition of the turbulent shear production given in (3). Such a direct ‘ad hoc’ discretisation would be for example: " #2 u^jþ1 u^j : ð17Þ Pjþ1=2 ¼ mjþ1=2 zjþ1 zj It can easily seen that (16) and (17) are identical only for the steady-state (^ uj ¼ uj ) case. Another ‘ad hoc’ discretisation could be the semi-implicit calculation for Pjþ1=2 : " #2 u~jþ1 u~j Pjþ1=2 ¼ mjþ1=2 ð18Þ zjþ1 zj which is identical to (16) for r ¼ 0:5. This discretisation is often applied when a leap-frog timestepping is used, see e.g. the Princeton Ocean Model (POM; Blumberg and Mellor, 1987). This semi-implicit choice for r is known to cause numerical instabilities, see e.g. Patankar (1980). Investigating this further here, is however beyond the scope of this paper. An equivalent formalism can be applied to the buoyancy equation (6). With the formalism given above, the discretisation is of the following form: b
b
b b
j j jþ1 j1 b^j bj mjþ1=2 zjþ1 zj mj1=2 zj zj1 ¼ 0: Dt zjþ1=2 zj1=2
ð19Þ
After multiplication with zj and some algebraic operations, the following equation for the discrete potential energy pj ¼ zj bj is obtained: b
b
b b
j j jþ1 j1 1 1 p^j pj 2 ðzjþ1 þ zj Þmjþ1=2 zjþ1 zj 2 ðzj þ zj1 Þmj1=2 zj zj1 þ Dt zjþ1=2 zj1=2 1 1 bjþ1 bj bj bj1 ¼ mjþ1=2 þ mj1=2 ¼: Bljþ1=2 Buj1=2 : 2 zjþ1=2 zj1=2 2 zjþ1=2 zj1=2
ð20Þ
The terms on the right hand side are the discrete expressions for the sinks and sources of potential energy due to vertical mixing of density. In order to achieve total energy conservation, these have to be used for the discretisation of the turbulent buoyancy production. The discrete buoyancy production for kjþ1=2 would thus be of the following form: bj bjþ1 ð21Þ Bjþ1=2 ¼ mjþ1=2 zjþ1 zj which again differs from a direct discretisation which would be formulated as: Bjþ1=2 ¼ mjþ1=2
b^jþ1 b^j : zjþ1 zj
ð22Þ
Formulations (21) and (22) are identical only for steady-state flow or fully implicit time stepping with r ¼ 1 and, consequently, bj ¼ b^j .
H. Burchard / Ocean Modelling 4 (2002) 347–361
353
Since the production terms involve the new velocity u^j and the new buoyancy b^j , the velocity and buoyancy equation are solved before the turbulent kinetic energy equation in each time step. Thus, the old eddy viscosity and diffusivity value have to be used in the velocity and buoyancy equation. Due to the completely independent derivation of the energy-conserving shear and buoyancy production terms, different values of r can of course be used for momentum and buoyancy. For temperature and salinity, which are usually calculated first as the two quantities defining buoyancy, the same r-values should be used. The physical model has to be completed by formulations for eddy viscosity mt , eddy diffusivity m0t and turbulent dissipation rate, e. Here, a model as simple as possible has been chosen in order to reduce the complexity of the calculation, but still retain basic physical properties. Eddy viscosity and diffusivity are calculated as mt ¼ cl
k2 ; e
m0t ¼
mt Pr
ð23Þ
with the empirical constant cl ¼ 0:09, see Rodi (1980) and the turbulent Prandtl number Ri Ri Pr ¼ Pr0 exp 0 1 þ 1 Pr Ri Ri
ð24Þ
2
with the gradient Richardson number Ri ¼ oz b=ðoz uÞ , the neutral turbulent Prandtl number Pr0 ¼ 0:74 and Ri1 ¼ 0:25, see Schumann and Gerz (1995). The dissipation rate is calculated as follows: e ¼ c3=4 l
k 3=2 L
ð25Þ
with the macro length scale L¼
1 1 þ k 1=2 jðz þ z0 Þ c N
!1 ;
ð26Þ
where the first term is a length limitation near the surface given by the law of the wall and the second term is a length limitation due to the effect of stable stratification. The surface roughness parameter is here set to z0 ¼ 0:1 m and the von Karman constant is given by j ¼ 0:4. The stability parameter c should depend on the stability functions cl and c0l and has been calibrated here to c ¼ 0:4 such that it reasonably reproduces the mixed layer depth evolution for the windentrainment experiment described below, see Eq. (30). This length scale parameterisation is similar to that of Axell and Liungman (2001) who used slightly different values for c. In order to guarantee positivity for the turbulent kinetic energy, the sink terms on the right hand side of the discrete k-equation have to be treated quasi-implicitly, see Patankar (1980). This leads to the following discrete form of the eddy diffusivity in the term Bjþ1=2 from (21): 8 k2 < cl jþ1=2 for bjþ1 < bj ; ejþ1=2 mjþ1=2 ¼ ð27Þ ^jþ1=2 kjþ1=2 k :c else: l
ejþ1=2
354
H. Burchard / Ocean Modelling 4 (2002) 347–361
An equivalent discretisation is applied when (22) is used for the buoyancy production term. Thus, the energy flux from the turbulent kinetic energy to the potential energy is not conservative for stable stratification. There is one other reason why the budget of the turbulent kinetic energy is not closed: In order to avoid too low values for the turbulent kinetic energy (which would lead to an underflow in numerical calculations), a lower limiting value is imposed, which is here: k > kmin ¼ 1010 J m2 :
ð28Þ
This has the consequence of a small artificial positive energy flux in areas with net turbulence dissipation. With (25), the dissipation term on the right hand side of the turbulent kinetic energy equation will be discretised according to Patankar (1980) as follows: ejþ1=2 ¼
c3=4 l
1=2 k^jþ1=2 kjþ1=2
Ljþ1=2
:
ð29Þ
This does however not violate the total energy conservation, since the turbulent dissipation rate is a sink of the total energy, see Eq. (12). The wind-mixing experiment used here is motivated by the laboratory experiment described by Kato and Phillips (1969). There, a mixed layer induced by a constant surface stress penetrates into a stably stratified fluid with density increasing linearly down from the surface. The water depth is assumed to be infinite. Price (1979) suggested a solution for the evolution of the mixed-layer depth Dm based on a constant Richardson number 1=2 1=2
Dm ðtÞ ¼ 1:05u N0
t
;
ð30Þ
where u is the surface friction velocity and N0 the constant initial Brunt–V€ais€al€a frequency. Following several authors, see e.g. Deleersnijder and Luyten (1994); Burchard et al. (1998) this laboratory experiment is transformed to ocean dimensions with u ¼ 102 m s1 and N0 ¼ 102 s1 . The empirical estimate (30) results in a mixed layer depth of Dm ¼ 34:5 m after 30 h of entrainment, which will be the run time of this numerical experiment. A total depth of H ¼ 50 m is chosen here which is sufficiently deep for simulating an infinite depth for this study.
4. Numerical performance study In this section, the numerical performance of the energy-conserving formulations for turbulent shear production and buoyancy production, (16) and (21) is investigated by means of the simple wind-mixing experiment given above and compared to the non-conservative schemes (17) and (22), respectively. The performance of the relatively simple turbulence closure model consisting of Eqs. (9), (23), (24), (25), and (26) is first compared to a more sophisticated closure model consisting of an additional transport equation for the dissipation rate e instead of using Eqs. (25) and (26). The stability functions in that model are based on an algebraic second-moment closure recently suggested by Canuto et al. (2001). The complete model has been described and tested in detail by Burchard and Bolding (2001). Fig. 2 shows profiles of eddy viscosity mt , turbulent kinetic energy k,
H. Burchard / Ocean Modelling 4 (2002) 347–361
355
Fig. 2. Profiles of eddy viscosity mt (left), turbulent kinetic energy k (middle), and macro length scale L (right) after 30 h of wind-entrainment. Calculations have been carried out with the one equation k model with algebraic length scale (––), see Eq. (26) and with the two-equation k–e model with algebraic second-moment closure by Canuto et al. (2001) (- - -).
and macro length scale L after 30 h of wind-entrainment calculated with both models. Despite the structural differences between the two turbulence closure models, the resulting turbulence profiles do not differ that much. This similarity to the well-tested, more complex model, gives some credibility to the simpler model with algebraic length scale. It is thus concluded that this model is sufficiently general in order to represent turbulence closure models for this numerical performance test. Three different resolutions in time and space have been used for this study: 1. A non-equidistant grid spacing with strong zooming towards the surface: tanh ðdl þ du Þ Mj dl þ tanhðdl Þ 1 ; j ¼ 0; . . . ; M zjþ1=2 ¼ H tanhðdl Þ þ tanhðdu Þ
ð31Þ
with the zooming coefficients du ¼ 3, dl ¼ 0 and the number of layers, M ¼ 200. This results in a near-surface layer depth of hM ¼ 0:008 m and a near bed layer depth of h1 ¼ 0:754 m. The time step chosen here is Dt ¼ 200 s. 2. An equidistant grid-spacing with M ¼ 200 (resulting in layer thicknesses of 0.25 m) and Dt ¼ 200 s. 3. An equidistant grid-spacing with M ¼ 200 (resulting in layer thicknesses of 0.25 m) and Dt ¼ 20 s. The level of implicitness has been set to r ¼ 0:6 for the momentum and the buoyancy equation, a value slightly larger than the semi-implicit value 0.5 ensuring second-order accuracy in time, but which is also known for causing numerical instabilities. For the turbulent kinetic energy equation, a fully implicit time-stepping with r ¼ 1 has been chosen in order to further stabilise the model. Fig. 3 shows profiles of eddy viscosity mt obtained with the new energy-conserving discretisation of turbulence production terms (16) and (21) and the old not energy-conserving discretisation (17) and (22), using the three spatial and temporal resolutions 1–3 given above.
356
H. Burchard / Ocean Modelling 4 (2002) 347–361
Fig. 3. Profiles of eddy viscosity mt obtained by applying the new ( ) and the old ( ) discretisations of turbulence production terms. The three spatial and temporal discretisations 1 (left), 2 (middle) and 3 (right) have been used.
It is clearly seen that the old discretisation is numerically unstable in the region of maximum eddy viscosity for non-equidistant grid-spacing with long time steps of Dt ¼ 200 s. In contrast to this, the new discretisation results in a stable and smooth solution which is not significantly deviating (by optical inspection) from the equidistant grid-spacing case with short of Dt ¼ 20 s, see Fig. 2 (left panel) and Fig. 3 (right panel). For equidistant grid-spacing and large time steps of Dt ¼ 200 s, the old method is numerically stable at the end of the run time, but significantly underestimates the mixed layer depth predicted obtained with high resolution. Only for equidistant grid spacing and small time steps of Dt ¼ 20 s, the old method results in a final viscosity profile optically not distinguishable from the profile obtained with the new method. It is also instructing to inspect some aspects of the temporal evolution of the numerical solutions. In Fig. 4, time series of the near surface velocity u at z ¼ 0:125 m are given. This coincides with the discrete velocities at M ¼ 200 for equidistant and M ¼ 187 for non-equidistant gridspacing given by (31) with du ¼ 3. In Fig. 4, the time axis is given in logarithmic scale in order to emphasise the instabilities occurring directly after the sudden onset of surface forcing over a laminar water column. These instabilities are visible for all model runs with long time steps of Dt ¼ 200 s. However, the near-surface velocity converges towards the result for Dt ¼ 20 s after 4 h for the non-equidistant case and after 1 hour for the equidistant case, when the new method is used. With the old method, some approximation of the high temporal resolution solution after about 6 h is visible, but even at the end of the run time, visible differences are present, see also Fig. 3. Finally, balances for the vertically integrated energy equations are plotted in Fig. 5 (new discretisation) and 6 (old discretisation) for non-equidistant grid-spacing and Dt ¼ 200 s. The total energy budget, Eq. (12), the mean kinetic energy budget, vertical integral of Eq. (3), the potential energy budget, vertical integral of Eq. (8), and the turbulent kinetic energy budget, vertical integral of Eq. (9) are shown. For the new discretisation, the total energy balance is basically an equilibrium between energy input due to surface stress, increase of mean kinetic energy due to mixing and turbulent dissipation, similarly to the mean kinetic energy budget where the energy input by wind is balanced by an increase of mean kinetic energy and shear production of turbulent kinetic energy. The turbulent kinetic energy budget is mainly a balance between shear production
H. Burchard / Ocean Modelling 4 (2002) 347–361
357
Fig. 4. Time series of near surface velocity u at z ¼ 0:125 m computed with the three spatial and temporal discretisations given above. Upper panel: new discretisation of turbulence production terms (16) and (21); lower panel: old discretisation of turbulence production terms (17) and (22). For time axis is given in logarithmic scale for clarity.
and dissipation, slightly modified by loss of turbulent kinetic energy to potential energy, for which the buoyancy production term fully balances the increase of potential energy. Some numerical instabilities are visible during the first three hours of the simulation for the shear production P, the dissipation rate e and the tendency of turbulent kinetic energy, ot k. The residuals for the mean kinetic and the potential energy are exactly zero (down to machine accuracy), but the turbulent kinetic energy budget and the total energy budget are not fully closed due to the two model features discussed above, the quasi-implicit treatment of sink terms due to Patankar (1980) and the lower limit set for the turbulent kinetic energy. The vertically integrated energy budgets for the old discretisation are shown in Fig. 6. It can be seen that the shear production is some orders of magnitude larger than the energy input due to surface wind stress. This is merely due to the non-conservative discretisation of the shear production term. This leads to an excessive production of turbulent kinetic energy, which is not fully compensated by dissipation, which depends on k and N, two variables which are lagging behind the shear production term. It seems that a positive feedback occurs in the mean kinetic energy budget during the first three hours of the simulation. Too large shear production causes too large turbulent kinetic energy, which in turn causes too large eddy viscosity. This again increases the
358
H. Burchard / Ocean Modelling 4 (2002) 347–361
Fig. 5. Results of the wind-entrainment experiment using the new, energy-conserving discretisations for the turbulence production terms, (16) and (21). Balance terms in the vertically integrated energy budget equations for total mechanic energy (upper left), mean kinetic energy (upper right), potential energy (lower left) and turbulent kinetic energy (lower right).
shear production. Also the potential energy budget is out of balance during the first five hours of the experiment due to the non-conservative discretisation of the buoyancy production term. It can be assumed that the unbalanced energy budgets cause the prevailing numerical instabilities seen in Figs. 3 and 4. However, a strict mathematical reasoning for this could not be found yet.
H. Burchard / Ocean Modelling 4 (2002) 347–361
359
Fig. 6. Results of the wind-entrainment experiment using the old, not energy-conserving discretisations for the turbulence production terms, (17) and (22). Balance terms in the vertically integrated energy budget equations for total mechanic energy (upper left), mean kinetic energy (upper right), potential energy (lower left) and turbulent kinetic energy (lower right).
5. Conclusions The new energy-conserving discretisations for the shear and buoyancy production terms of turbulent kinetic energy presented here differ from ‘ad hoc’ discretisations by the level of implicitness which is used for the vertical differences. While the procedure for the linear buoyancy
360
H. Burchard / Ocean Modelling 4 (2002) 347–361
production term seems to be straight-forward, the non-linear shear production term is based on a product of differences which are taken on different levels of implicitness. Although the differences between the ‘new’ and the ‘old’ formulations are small, their consequences for the overall performance of the numerical model are dramatic. After some initial instabilities due to the shock-like model initialisation, the ‘new’ discretisations compute even for strong zooming and large time steps numerically stable results for time series of near-surface velocity and also for profiles of eddy viscosity, a parameter known for its sensitivity to numerical instabilities, see Burchard and Deleersnijder (2001). In contrast to this, the initial instabilities do not vanish when the ‘old’ discretisation is used. Further experiments (not shown here) using twoequation k–e models with complex algebraic second-moment closure proved a similar stability behaviour than the simpler one-equation model discussed in this paper. It is assumed that also other turbulence models such as the k–kL two-equation model by Mellor and Yamada (1982) or the k–x model with x ¼ e=k, see Wilcox (1998) would profit from the new production term discretisations. It is satisfying that a new discretisation which ensures basic physical principles also for the numerical solution is indeed numerically more stable and thus also of practical relevance. An implementation of this new method into three-dimensional ocean models is straight-forward. It should however be noted that it is usually not the buoyancy b which is prognostically calculated in ocean models but the potential temperature and salinity from which the buoyancy is then diagnostically calculated by means of a non-linear equation of state. Thus, the energy flux between the potential energy and the turbulent kinetic energy is not completely closed, which is a consequence of the Boussinesq approximation usually made in ocean modelling.
Acknowledgements This study has been funded through a Habilitation grant by the Deutsche Forschungsgemeinschaft (German Research Foundation). Further support has been provided by the CARTUM (Comparative Analysis and Rationalisation of Second-Moment Turbulence Models) project, a concerted action of the MAST-III program of the European Commission (MAS3-CT98-0172). The author is indebted to Eric Deleersnijder (Universite Catholique de Louvain, Institut d’Astronomie et de Geophysique George Lema^ıtre, Belgium) for detailed discussions on this topic. Parts of the work presented in this paper have been carried out during a stay of Hans Burchard at this institute. The comments of two anonymous referees were greatly appreciated. For the numerical calculations in this study, of the General Ocean Turbulence Model (GOTM, http:// www.gotm.net) has been used.
References Arakawa, A., Lamb, V.R., 1977. Computational design of the basic dynamical processes of the UCLA general circulation model. Meth. Comput. Phys., 173–263. Axell, L., Liungman, O., 2001. A one-equation turbulence model for geophysical applications: comparison with date and the k–epsilon model. Environ. Fluid Mech. 1, 71–106.
H. Burchard / Ocean Modelling 4 (2002) 347–361
361
Blumberg, A.F., Mellor, G.L., 1987. A description of a coastal ocean circulation model. In: Heaps, N.S. (Ed.), Three dimensional Ocean Models. American Geophysical Union, Washington D.C, pp. 1–16. Burchard, H., Bolding, K., 2001. Comparative analysis of four second-moment turbulence closure models for the oceanic mixed layer. J. Phys. Oceanogr. 31, 1943–1968. Burchard, H., Deleersnijder, E., 2001. Stability of algebraic non-equilibrium second-order closure models. Ocean Modell. 3, 33–50. Burchard, H., Petersen, O., Rippeth, T.P., 1998. Comparing the performance of the Mellor–Yamada and the k–e twoequation turbulence models. J. Geophys. Res. 103, 10543–10554. Canuto, V.M., Howard, A., Cheng, Y., Dubovikov, M.S., 2001. Ocean turbulence I: one-point closure model. Momentum and heat vertical diffusivities. J. Phys. Oceanogr. 31, 1413–1426. Deleersnijder, E., Luyten, P., 1994. On the practical advantages of the quasi-equilibrium version of the Mellor and Yamada level 2.5 turbulence closure applied to marine modelling. Appl. Math. Modell. 18, 281–287. Haidvogel, D.B., Beckmann, A., 1999. Numerical Ocean Circulation Modelling. Series on Environmental Science and Management, vol. 2. Imperial College Press, London. Kato, H., Phillips, O.M., 1969. On the penetration of a turbulent layer into stratified fluid. J. Fluid Mech. 37, 643–655. Mellor, G.L., Yamada, T., 1982. Development of a turbulence closure model for geophysical fluid problems. Rev. Geophys. 20, 851–875. M€ uller, P., 1995. Ertel’s potential vorticity theorem in physical oceanography. Rev. Geophys. 33, 67–98. Patankar, S.V., 1980. Numerical Heat Transfer and Fluid flow. McGraw-Hill, New York. Price, J.F., 1979. On the scaling of stress-driven entrainment experiments. J. Fluid Mech. 90, 509–529. Rodi, W., 1980. Turbulence models and their application in hydraulics. Tech. Rep., Int. Assoc. for Hydraul. Res., Delft, The Netherlands. Schumann, U., Gerz, T., 1995. Turbulent mixing in stably stratified shear flows. J. Appl. Meteorol. 34, 33–48. Wilcox, D.C., 1998. Turbulence Modelling for CFD, 2nd ed. DCW Industries.