Laminar and weakly turbulent oceanic gravity currents performing inertial oscillations

The small scale dynamics of a weakly turbulent oceanic gravity current is determined. The gravity current considered is initially at rest and adjusts by performing inertial oscillations to a geostrophic mean flow. The dynamics is explored with a hierarchy of mathematical models. The most involved are the fully 3-D Navier-Stokes equations subject to the Boussinesq approximation. A 1-D and 0-D mathematical model of the same gravity current dynamics are systematically derived. Using this hierarchy and the numerical solutions of the mathematical models, the turbulent dynamics at the bottom and the interface is explored and their interaction investigated. Three different regimes of the small scale dynamics of the gravity current are identified, they are characterised by laminar flow, coherent roll vortices and turbulent dynamics with coherent streaks and bursts. The problem of the rectification of the turbulent fluxes, that is, how to average out the fluctuations and calculate their average influence on the flow, is considered. It is shown that two different regimes of friction are superposed, an Ekman friction applies to the average geostrophic flow and a linear friction, not influenced by rotation, to the inertial oscillations. The combination of the two makes the bulk friction non-local in time for the 0-D model. The implications of the results for parametrisations of the Ekman dynamics and the small scale turbulent fluxes in the planetary boundary layer are discussed.


Introduction
Oceanic gravity currents show a variability over a wide range of scales in space and time.To leading order they are governed by geostrophic equilibrium, that is a balance between the buoyancy and Coriolis force (Griffiths, 1986;Coleman et al., 1990;Price and Baringer, 1994;Killworth, 2001;Wahlin and Walin, 2001;Wirth, 2009).A gravity current in this equilibrium follows the lines of constant depth and does not descend.This equilibrium is perturbed by principally two processes: (i) at the meso-scale by its instability leading to wavelike disturbances and eddies; (ii) at even smaller scales the flow is three dimensional and turbulent, which leads to turbulent fluxes of mass and momentum.The quasi 2-D mesoscale dynamics is assumed to be well represented in today's high-resolution hydrostatic numerical models of the ocean dynamics and it is not investigated here.The appearence of meso-scale instability and variability is hindered by the homogeneous initial conditions in the direction parallel to the ocean floor and it is excluded in the numerical integrations by the small size of the domaine of integration.The subject of the paper is the second point the small scale turbulent dynamics which is fully 3-D, non-hydrostatic and involves scales smaller than a metre in all spatial directions.The influence of this small scale turbulent dynamics on the large scale has to be parametrised in today's and tomorrow's numerical models of the ocean dynamics, as they do not and will not explicitly resolve it.The small scale turbulent dynamics and its influence on the large scale dynamics is the subject of the present work.It is important to realise, that the small scale turbulence in gravity currents acts directly on the large-scale dynamics of the gravity current as it governs the friction, mixing and entrainment processes.
Published by Copernicus Publications on behalf of the European Geosciences Union.

A. Wirth: Laminar and weakly turbulent oceanic gravity currents
There is a substantial number of publications on the dynamics of gravity currents considering observations, laboratory experiments, analytical models and calculations and numerical simulations based on hydrostatic and non-hydrostatic mathematical models.For a comprehensive review on oceanic gravity currents I refer the reader to Griffiths (1986) and Price and Baringer (1994), concerning turbulent fluxes I direct the reader to the recent and comprehensive review by Legg et al. (2009) and the references therein.For numerical studies of turbulent fluxes in non-rotating gravity currents please see Özgökmen et al., 2006.The adjustment process of a gravity current to a topographic slope by performing inertial oscillations has been studied by Nof (1996) in the context of an oceanic turbidity current.When an oceanic gravity current moves along a topographic slope it does so at an average speed that is close to the geostrophic equilibrium.When the topographic slope changes the gravity current adjusts to the new slope by performing inertial oscillations, this is the nature of geostrophic adjustment known since the work of Rossby (1936).The process of inertial oscillations, the so-called fast dynamics, is often neglected although it is a dominant signal in the energy spectrum of the world oceans (Ferrari and Wunsch, 2009).To my understanding this neglect has principally two reasons: first, a substantial body of theoretical and numerical research were and are performed using quasi-geostrophic models, where inertial dynamics is filtered out by the mathematical model, ahead of the numerical integration.Second, in numerical studies based on primitive equation models snapshots are printed over time intervals substantially exceeding the inertial period, or even presenting mean quantities over time-intervals spanning many inertial periods.In such a way inertial dynamics present in the calculations is filtered out following the numerical integration.
A turbulent Ekman layer forms at the bottom of a gravity current.The turbulent Ekman layer below a stationary geostrophic flow has been studied numerically starting from the pioneering work of Coleman et al., 1990 andColeman, 1999.In their work the data were averaged over several inertial periods to average out the inertial oscillations which were present even in their case, with a flow that was initially adjusted.Inertial oscillations of oceanic gravity currents are also discussed in detail by Wang et al. (2003) using a 1-D 2-layer model.The here presented gravity current is initially at rest and adjusts to the sloping bottom by performing inertial oscillations.These oscillations, which are damped by friction, are also initiated each time a gravity current adjusts to a changing topographic slope or roughness and are thus likely to be a recurrent feature of oceanic gravity currents.If the dynamics is described by linear equations, these oscillations average out and leave no imprint on the slow dynamics (averaged over one or several inertial period).When nonlinearity becomes important these oscillations will influence the slow dynamics.It is shown here, that inertial oscillations have a strong influence on the turbulent fluxes.These turbu-lent fluxes lead to turbulent transport of mass and momentum which have to be parameterised in models of the ocean circulation.
The dynamics of the gravity current is investigated with a hierarchy of mathematical models of three, one and zero spatial dimensions.The purpose is twofold.First, the comparison of the results allows a deeper understanding of the processes and the evaluation of their impact on the slow dynamics of a gravity current.Second, in ocean general circulation models (OGCMs) the gravity current usually spans only part of a grid box in the vertical and its numerical (non-) resolution resembles the 0-D-model.In OGCMs with a very high vertical resolution at the bottom (see Laanaia et al., 2010) the dynamics is similar to the 1-D-model.

The physical problem considered
In the present work I consider the dynamics of a three dimensional gravity current on an inclined plane of constant slope α in a rotating frame, with a constant Coriolis parameter f = 4π/t 0 , where t 0 is the rotation period of the frame.Please see Fig. 1 for the configuration considered.The extension of the gravity current in the up slope and along-slope (along isobaths) direction is assumed to be very large as compared to its thickness.
Major oceanic gravity currents (strait of Gibraltar or the Denmark strait) extent several hundreds of kilometres in the along-flow direction, several tenths of kilometres in the cross-flow direction and are only around 100 metres thick Price and Baringer (1994).For the dynamics of gravity currents constrained in wide channels and the effect of channel geometry on the gravity current I refer the reader to the recent work of Umlauf and Arneborg (2009a, b), Cossu et al. (2010) and Wahlin (2002Wahlin ( , 2004)).The small scale dynamics of gravity currents considered here is for a gravity current on a flat inclined plane.Within the gravity current all variables are assumed to be statistically homogeneous in the directions parallel to the inclined plane.This means that time averaged quantities depend only on their distance to the ocean floor.The initial temperature anomaly of the gravity current, with respect to the surrounding water, is T 0 = −0.5K(< 0), H = 100 m is the thickness of the gravity current and h = 20 m measures the thickness of the interface between the gravity current and the surrounding water.This configuration is chosen to study the local dynamics of a turbulent gravity current neglecting the effect of large-scale gradients in the directions parallel to the plane, which are usually smaller than in the direction perpendicular to the plane.It has been shown by Wahlin and Walin (2001) and Wirth (2009) that the dynamics due to the large-scale gradients is well described by geostrophic balance subject to bottom friction, when the mesoscale dynamics is suppressed.In the present work I assume that large-scale gradients in the directions along the plane have a negligible influence on the n: A layer of cold (dense) water of thickness H superposed by warmer (less dense) water on an inclined plane of me with frequency Ω.The interface (blue line) has a thickness h.The temperature anomaly is indicated by the red by an angle α of the direction of gravity with respect to the z-direction.The gravity vector is in the x − z-plane.

Fig. 1. Initial condition:
A layer of cold (dense) water of thickness H superposed by warmer (less dense) water on an inclined plane of angle α in a rotating frame with frequency .The interface (blue line) has a thickness h.The temperature anomaly is indicated by the red line.Please note the tilt by an angle α of the direction of gravity with respect to the z-direction.The gravity vector is in the x-z-plane.local turbulent dynamics.The local turbulent dynamics govern, however, the large scale dynamics as it determines the friction laws and parameters and the dilution of the gravity current water (see Wirth, 2009).
The thermal expansion coefficient γ leads to a reduced gravity of g 0 = −gγ T 0 , where g is gravity.The problem considered here depends on 7 external parameters: (α, f, g 0 , h, H, ν, κ), where the last two are the kinematic viscosity and diffusivity, respectively.When dissipative effects are ignored the problem has a stationary solution with a geostrophic balance between the force of gravity and the Coriolis force leading to a geostrophic velocity of: The geostrophic speed v G is also called the Nof-speed (Nof, 1983).The "sin α" instead of the common "tan α" term appears as the rotation vector is perpendicular to the slope rather than aligned with gravity.For the small angle (α = 1 • ) considered here the difference is negligible (see Fig. 1).I will assume that the total depth of the fluid is much larger than H and has no dominant influence on the local turbulent dynamics.Several non-dimensional parameters can be formed which describe the dynamics of the problem.First is the Froude number F r = v G / g 0 H = g /H /f (sin α) comparing the geostrophic fluid speed v G to the speed of shallow water waves at the interface.The Froude number was found to be an important parameter in gravity currents and Killworth (2001) suggests that its value is always below ≈ 0.8, as gravity currents with larger Froude numbers are unstable and subject to turbulent fluxes that reduce the Froude number by entraining surrounding water.Please note that F r ∼ g 0 and F r ∼ H −1 , so that entrainment, which decreases g 0 and increases H , does decrease the Froude number efficiently.Another important parameter is the (gradient) Richardson number Ri = g 0 h/v 2 G giving the ratio between the stabilising effect of stratification to the destabilising effect of velocity shear.Please note that in my definition the Froude number is based on the thickness of the gravity current H and the Richardson number is based on the thickness of the interface h.Flows with a local Richardson number Ri < 0.25 are usually found to be unstable and become turbulent, whereas in flows with Ri > 1 turbulence is suppressed, please see Galperin et al. (2007) and Zilitinkevich et al. (2008) for a review and a critical discussion of the subject.It is important to note that in this definition the Froude number and the Richardson number are independent parameters and the ratio h/H = F r 2 Ri.The fundamental difference between the two is ignored in many research papers.When the effects of viscosity and dissipation are neglected there are three independent dimensionless parameters: (α, F r, Ri).The Rossby number, based on the geostrophic velocity and the layer thickness, is a function of these three non dimensional parameters, Ro = v G /(f H ) = F r 2 / sin α.The thus constructed Rossby number is large (as α 1), indicating that non-linearity is not dominated by rotation and that 3-D turbulence is likely to be important in the problem considered here.
A horizontal length scale is given by the Rossby radius L = g H /f = F rH / sin α which indicates the length scale for which there is a resonance between rotation and wave motion.In this work the mean motion is at infinite horizontal length scale and the dynamics studied at horizontal length scales which are smaller or comparable to the thickness H and thus much smaller than L. In the mathematical model employed (see Sect. 3) and in its numerical implementation (see Sect. 4), all variability on larger horizontal scales is suppressed by the periodicity L x and L y in the x and y-direction, respectively.This is beneficial to our goal of exploring the small scale turbulent fluxes.The Rossby radius L is much larger than the domain size and the appearance of mesoscale structures is thus artificially suppressed and L is not an important parameter in our experiment.The mesoscale dynamics is usually well represented in todays regional high-resolution hydrostatic ocean models.
All other independent non-dimensional variables involve dissipative variables, that is viscosity ν and/or diffusivity κ.They will in the following be called explicit viscosity/diffusivity as they appear explicitly in the governing equations and differ from the eddy viscosity/diffusivity ν eddy and κ eddy which are due to the resolved small scale dynamics.Their ratio is called the Prandtl number P r = ν/κ.The Reynolds number based on the velocity shear across the interface is Re h = v G h/ν.The thickness of the laminar Ekman layer δ = √ 2ν/f leads to an Ekman number Ek = (δ/(2H )) 2 , which gives the ratio of viscous force to the Coriolis force, in the case of a laminar boundary layer.The Reynolds number based on the layer thickness is Re H = H v G /ν and thus Ek = Ro/(4Re H ). A Reynolds number based on the laminar Ekman layer thickness is it compares the importance of the non-linear to the viscous term in the Ekman layer.
At the bottom of the gravity current a turbulent planetary boundary layer (PBL) forms.The dynamics of such PBLs is well studied in the case of a constant-in-space-and-time current above the PBL.The length scales of the turbulent boundary layer depend on the friction velocity u * = √ τ/ρ, where τ is the average friction force per unit surface area exerted by the fluid on the boundary and ρ is the density of the fluid.The ratio of the friction velocity and the fluid velocity above the boundary layer is the square root of the geostrophic drag coefficient The length scale z 0 is equal to the larger of the roughness of the ocean floor and ν/u * .The turbulent PBL is characterised by the surface Rossby number, which is related to the Reynolds number based on the Ekman layer thickness.Please note that Ro * (u * , c D ) are a result of the experiment and not an initial parameter.
The PBL can be decomposed into four layers (see Fig. 2).The first is the viscous sub-layer with a thickness of about 5z 0 and a horizontal velocity that varies linearly with the distance from the floor.In the buffer layer above, the dynamics transits to the turbulent log-layer which starts around 20z 0 .At a distance 0.1δ * Ek , with δ * Ek = u * /f , rotation becomes important and a turbulent Ekman layer with an extension to about 0.5δ * Ek forms (see Ferrero et al., 2005 andMcWilliams, 2006, for a concise introduction).Above this extends the quasi-geostrophic interior.The transport of momentum between the ocean floor and the viscous layer, as well as within the viscous layer is done by molecular friction.Above, turbulent transport takes over.The Ekman-layer influences the quasi-geostrophic interior by vertical advective transport (Ekman pumping), created through divergence of the horizontal Ekman transport (see Wirth, 2009).
The above considerations apply to the turbulent stationary PBL.Please note that in the problem considered here the dynamics above the PBL is oscillatory in time with a frequency which is equal to the Coriolis parameter f , which is also the adjustment time of the Ekman layer dynamics.In such configuration the PBL dynamics might be strongly altered.
Turbulent fluxes at the interface include mixing, entrainment and detrainment.They can be induced by: (i) local instability of the interface due to a low Richardson number (Kelvin-Helmholtz instability), or (ii) hydraulic jumps for Froude numbers exceeding unity, or (iii) by turbulence from the bottom boundary layer reaching the interface.Turbulent fluxes due to (i) are symmetric about the interface (in a Boussinesq fluid) whereas (ii) and (iii) can lead to asymmetric fluxes such as entrainment.The relative importance of (ii) versus (iii) can be expressed by: 0 This leads to critical angles α crit = arcsin( √ c G /2) which lies between one and two degrees, slopes typical for oceanic gravity currents.
The problem considered here depends on five independent dimensionless parameters: (α, Ri, F r, P r, Re δ ).

A hierarchy of mathematical models
The above introduced physical problem is studied with a hierarchy of mathematical models of three, one and zero dimensions.

Geometry of the initial configuration
The Cartesian coordinate system is chosen such that the inclined plane is given by z = 0.This implies that the x-y-plane is not normal to gravity.Normal to the x − y-plane is the zdirection, it forms a (small) angle α to the direction of gravity (see Fig. 1).
The gravity current is initially at rest: The initial temperature anomaly, with respect to the surrounding water, of the gravity current is given by: The gravity current is assumed to be periodic in the x-and ydirections with a periodicity L x and L y , respectively.Such a geometry allows a rigorous definition of averages in the x-y-plane and facilitates its numerical implementation.

The three dimensional model
The dynamics of and incompressible fluid is described by the three dimensional Navier-Stokes equations subject to the  Boussinesq approximation where u, v and w are the velocity components in the x, y and z-directions, respectively.P is the pressure, ν the explicit viscosity, g gravity, and ∇ 2 = ∂ xx + ∂ yy + ∂ zz is the Laplace operator.The boundary conditions are free slip at the upper surface and no-slip at the lower boundary: The equation of a scalar transported by a fluid is where T is temperature, κ is the explicit diffusivity of temperature.The boundary conditions for the temperature are no-flux at the lower and upper boundary: The boundary condition for all variables are periodic in x and y with the period L x and L y , respectively.The linear equation of state: allows to obtain the density anomaly from temperature using a constant thermal expansion coefficient γ .We further denote α c = gγ cos α and α s = gγ sin α

The one dimensional model
The problem is statistically homogeneous in the x and y directions, which allows to replace ensemble averages by averages in the x-y-plane.This leads to averaged quantities that depend only on the z-direction.We introduce the averages in the x-y-plane: When the Navier-Stokes equations are averaged over the xy-plane they read: Where all averages depend only on time and the zcoordinate.
The basis of our investigations are Eqs.( 16), ( 17) and ( 18).When the turbulent fluxes vanish, the model is linear.Wang et al., 2003 give an analytic solution of a two layer linear model when diffusion is neglected.Such a solution shows a persistent interfacial Ekman layer which is not a realistic feature, as can be verified in my numerical solutions of the 3-D model below and as it is discussed in Wirth (2011).The broadening of the interfacial dynamics is a key feature of gravity currents already noted by Ellison and Turner (1959).As I am not aware an analytical solution of this linear model with the initial and boundary conditions used here and a nonvanishing diffusivity, it is solved numerically.The numerical solution is close and converges towards a motion which is the sum of the geostrophic motion and inertial oscillation, the analytical solution for this motion is given in appendix A. An interesting feature is, that the inertial oscillation does not posses an Ekman spiral and no veering of the velocity vector at the bottom.This means that the friction force is always directed against the direction of motion.The nonoscillating part clearly shows a veering of the Ekman spiral at the bottom.
Although the one-dimensional model seems simpler than the full three-dimensional model, it suffers from the fact that it is not closed, that is, there are more unknowns than equations.This problem arises from the fact, that in general: Finding a closure for the above equations means expressing the vertical turbulent fluxes ( wu , wv and wT ) in terms of the calculated quantities ( u , v and T ).The most popular closure is introducing a eddy viscosity/diffusivity, a so called K-closure.This idea, which goes back to Prandtl (1925).In principal, the eddy viscosity is a fourth order tensor (∂ t u i + ... = ∂ l ν iklj ∂ k u j see Wirth et al., 1995), when averaged over the x-y-plane it reduces to a second order tensor.Please see Wirth 2010 for a discussion on anisotropic viscosity and its effects on the Ekman dynamics.The eddyviscosity closure in its anisotropic form is: If the eddy coefficients (ν u E , ν v E , ν uv E , κ E ), can be related to the 7 external parameters (see Sect. 2) and the variables averaged in the x-y-plane ( u , v , T ) the parameterisation problem is solved.As the small scale flow is very anisotropic, often formed by a roll structure approximately elongated in the flow direction, we can not expect ν u E (z) = ν v E (z) and ν uv E (z) = 0.This also means, that we can not determine the actual eddy-viscosity tensor from only observing average quantities as we have three unknowns in the only two momentum equations.
Variants of the 1-D model and the 0-D model in the next section are already discussed in Arneborg et al. (2007).

The zero dimensional model
The Eqs. ( 16), ( 17) and ( 18) can be averaged over the thickness of the gravity current H , these equations contain two difficulties.The first is that the thickness of the gravity current itself is subject to change due to entrainment, detrainment and mixing.The second are the turbulent fluxes at the interface, which appear in these equations.It is less involved to consider averages over the whole depth D of the fluid.This means integrating Eqs. ( 16), ( 17) and ( 18) over the whole depth D (see Fig. 1).The average along the zdirection is denoted by .= D 0 .dz/D.The double average ¯ . is an integral over the volume of the whole fluid, gravity current plus fluid above.The depth average of Eqs. ( 16), ( 17) and ( 18) becomes: Equation ( 25) expresses the conservation of heat.The volume averaged velocity vector depend only on time, it is a zero-(space)-dimensional problem.The problem lies in the determination of the velocity gradients at the bottom, which is written on the right hand side of Eqs. ( 23) and ( 24).The gradients for the 1-D and 3-D model can be seen in Fig. 5, which shows the velocity 0.5 m above the ground.
When the bottom friction is neglected, that is the right hand sides of Eqs. ( 23) and ( 24) are zero, the solutions are given by where the constants A and B are determined by the initial conditions, the motion is a geostrophically equilibrated constant velocity (g /f ) in the y direction which is superposed by inertial oscillations of amplitude √ A 2 + B 2 and frequency f .The dynamics of such a gravity current, which is initially at rest (A = −α s ¯ T /f , B = 0) is shown in Fig. 3.It is already presented in Nof (1996).
The bottom friction term can be parameterised by either a constant friction force (F x , F y ), a linear friction law (Rayleigh friction) r(u, v) where r is a friction time or by a quadratic drag law (turbulent friction) c G /H where c G is a drag coefficient (Wirth and Verron, 2008).Clearly a constant friction force, not depending on velocity, is not common, but its use is motivated by results presented in Sect. 5.For the constant and linear friction case an analytical solution exists: An important point is that the average transport has a downslope component, when friction opposes motion.It is r/ftimes the long-slope transport when the friction is linear.The dynamics for (−A = α s ¯ T = f = 1, B = F x = F y = 0 all non-dimensional) and different values of r/f are shown in Fig. 3.The downslope movement releases potential energy which transforms into kinetic energy and is drained by friction and dissipation.Friction also damps the amplitude of the inertial oscillations.Fig. 3 shows that for values of r/f exceeding ten percent the inertial oscillations are efficiently damped after only one inertial period.These models and solutions including a frictional process are, however, more mathematical curiosities and their physical meaning is questionable, as the velocity gradient at the boundary is unlikely to be linearly or quadratically related to the average velocity in the domain.Both representations of the friction also do not take into account the veering, turning, of the velocity vector in the PBL which leads to a friction force that is not aligned with the average velocity.Such friction terms added in Eqs. ( 23) and (24) will, if linear r ⊥ (−v, u), only change the the Coriolis parameter f = f + r ⊥ , as it is acting perpendicular to the average velocity like the Coriolis force.The veering and friction exerted by the floor is further considered in subsection 5.3.If we had a relation that connects volume averages of the velocities to the gradients at the ocean floor, the problem of finding a parametrisation would be solved.It is unlikely that such a relation can be found, there is not enough information included in the 0-D variables.Indeed, the equations do not even contain the thickness of the gravity current and its temporal evolution.

The 3-D numerical experiment and discussion of DNS approach
Strongly anisotropic grids are justified when only the large scale dynamics of oceanic gravity currents is considered explicitely as, at such scales, the dynamics is clearly anisotropic.At small scales, however, the turbulent structures are almost isotropic, a fact that has to be reflected in the numerical grid employed, which is x = y = 1 m and z = 0.5 m.The domain spans 256 m in the x-and ydirection and 223.5 m in the z−direction.The time step is 1 s.The viscosity is isotropic ν = 5 × 10 −3 m 2 s −1 and the Prandtl number is unity, that is diffusivity equals viscosity κ = ν.I checked that these values are sufficient to avoid a pile up of small scale energy caused by an insufficient viscous dissipation range, leading to a thermalized dynamics at small scales as explained by Frisch et al., 2008.The viscosity and diffusivity employed is more than three orders of magnitude larger than the molecular values of sea water, as the numerical grid resolution is too coarse to resolve the dynamics down to the dissipation scale (millimetre).Using large eddy simulation (LES) rather than the direct numerical simulation (DNS) employed here is clearly an option.In LES calculations, however, a turbulence model has to be used in the interior and a wall model at the boundary, somehow imposing the dynamics I like to study here, making the use of LES calculations controversial.A DNS that resolves the larger eddies, parameterising the effects of smaller eddies by constant eddy coefficients, as done here, might well be the best LES.Please see also Coleman andFerziger, 1990 andDimotakis, 2000 for a discussion of this point.The coherent structures (rolls and streaks) in the boundary layer studied here, have a typical size of tens of meters so that their viscous decay time t ν = l 2 /ν is on the order of ten hours, a little smaller than the inertial period, the dominant time scale in  the system.This shows, that the calculations are at the border of completely resolving these energy containing scales and the here presented results on the dynamics of the coherent energy containing structures can be extrapolated to higher Reynolds number flow.But, also note that the flow has not passed the mixing transition to a turbulent flow with an inertial range (see Dimotakis, 2000) as the Reynolds numbers are below 10 4 (see Table 2), which means that the route to dissipation of momentum and density gradients can not be extrapolated to higher Reynolds number flows.But I personally have doubts that this can be achieved by today's LES.It is found in laboratory experiments and numerical simulations of the turbulent wall layer that elongated coherent structures, streaks, are separated by around 100-times the inner scale z 0 of the boundary layer (see e.g.Robinson, 1991).The used domain size is only a little over twice this separation scale.Numerical experiments of larger domains with the same resolution, or finer resolution runs of the same domain, over a sufficiently long period (several inertial periods) are not attainable with my actual computer resources.A simulation of a complete gravity current measuring several tenths of kilometres in both horizontal directions is far beyond our actual and foreseeable future computer resources.
The parameter values of the experiment performed are given in Table 1.The definitions are given in Sect. 2.
These parameters lead to a reduced gravity of g 0 = 9.81 × 10 −4 ms −2 and a geostrophic velocity of v G = 1.71 × 10 −1 ms −1 .The laminar Ekman layer thickness is δ Ek = 10 m and the Reynolds number based on the Ekman layer thickness is Re δ = 340.The experiment starts from rest (u = v = w = 0 at t = 0) a small amplitude (.1T 0 ) white noise is added to the initial temperature anomaly.Snapshots of of all dynamical variables are printed every 5min.The total integration time is t max = 600 × 5 min= 50 h, which is almost 3 inertial periods (t inertial = 17.45 h).
The non-dimensional parameters are given in Table 2.
If we assume a geostrophic drag coefficient of c G = 10 −3 the friction velocity is u * ≈ 5.10 −3 ms −1 leading to a z 0 = 1m .
The numerical model used is HAROMOD (Wirth, 2004).HAROMOD is a pseudo spectral code, based on Fourier series in all the spatial dimensions, that solves the Navier-Stokes equations subject to the Boussinesq approximation, a no-slip boundary condition on the floor and a free-slip boundary condition at the rigid surface.The boundary conditions in the z-direction are imposed using a method based on the vertical boundary technique (Wirth, 2004).The boundary conditions in the x-and y-directions are periodic.The time stepping is a third-order low-storage Runge-Kutta scheme.
Please note that the large scale dynamics and the small scale turbulent dynamics have a large scale separation in space and time, which asks for substantial computer power in the direct numerical simulations.

Results
We have introduced a zero-dimensional, a one-dimensional and a three-dimensional mathematical model in Sect.3.Although the zero-dimensional and the one-dimensional model seem simpler than the full three-dimensional model they suffer from the fact that they are not closed, that is, they contain more unknowns than equations.The three-dimensional model presents a closed system and its results, when averaged over two or three dimensions and compared to the onedimensional and the zero-dimensional model, can be used to evaluate the closures.
The zero dimensional model allows for an analytical solution if we assume that the bottom friction is linear.An analytical solution that is close to the 1-D model is presented.For the one and three dimensional models, solutions are obtain numerically.

Qualitative description of 3-D results
In Fig. 6 the averages in the x-y-plane of the x-and ycomponent of the velocity vector are shown.They clearly show the inertial oscillations, as discussed in detail in Sect.3.4.They also give a first impression of the Ekman dynamics near the bottom boundary.Note, for example, that near the ocean floor the velocity is always down-slope, although above it changes sign periodically.
In Fig. 7 isosurfaces of the z-components of the velocity are shown.In a laminar gravity current dynamics this velocity component vanishes.They clearly show different nonlinear regimes characterised by laminar flow (not shown), followed by an instability leading to the well known roll structures in the Ekman layer (Fig. 7 left), which become unstable to a secondary instability with an about four times smaller wave length, Fig. 7 middle (see e.g.Dubos et al., 2008).Further instabilities lead to a turbulent boundary layer with bursts and streaks, (Fig. 7

Volume averages (zero dimensional results)
The variation of the volume integrated temperature anomaly during the total length of the experiment is less than 10 −4 times the initial value confirming the accuracy of the numerical advection scheme and the immersed boundary conditions.The time evolution of the volume averaged velocity components ¯ u and ¯ v are shown in Fig. 8. Their evolution is fitted by (see Eqs. (28, 28): The parameters that best fit the data from the simulations are shown in Table 3 and a comparison is found in Fig. 8.If friction were linear (Rayleigh friction) we would obtain B = −r/f A ≈ −1.5 × 10 −3 ms −1 , a result which is significantly different from value obtained by the best fit.This demonstrates that a simple linear Rayleigh friction based on the volume averaged velocity does not represent a good closure for the momentum flux at the bottom.Figure 8 also shows, that the frequency of the oscillation is indistinguishable from the Coriolis frequency f , demonstrating that there is no significant Ekman veering (r ⊥ ) associated to the oscillations.This is confirmed in Fig. 5, where the veering related to the mean motion but not to the inertial oscillations is depicted, and that the oscillations at different hight above the floor are in phase.See subsection 5.3 and appendix A for further explication concerning the absence of a veering for the inertial oscillations.This result, that the nature of the friction depends on the time-scale of the motion, leads to a friction which is non-local in time for a 0-D model.We continue by considering the volume averages of the kinetic and potential energy.The qualitative behaviour of the energy cycle in this experiment is clear: Potential energy is converted by the down-slope movement of gravity current water to kinetic energy, which is then irreversibly drained by diffusion and friction.This conversion is however far from being direct.As the gravity current water that is initially at rest slides down the incline, potential energy is converted to kinetic energy, the Coriolis force deviates the flow leading to inertial oscillations.The result is a periodic transfer between potential energy and kinetic energy.In the absence of dissipative effects the gravity current performs inertial oscillations converting potential to kinetic energy and back.These processes act at large scales, which are infinite in the x-ydirections, due to the periodic boundary conditions, and span the thickness of the gravity current in the vertical.Diffusion and friction, however, only act at the smallest scales.There is an infinite separation of scales in the horizontal, but not in the vertical direction.It is the vertical dynamics that allows for an energy cascade to the viscous/diffusive scales.The same conceptual picture might well apply to the ocean dynamics in general.The increase of potential energy due to diffusion in the vertical is negligible in the present experiment, that is, the mixing efficiency, which is the ration of energy used by diffusion to the energy injected, is close to zero.
The main quantities are the mean kinetic energies for unit mass where V 0 = 1 m 3 .The evolution of the potential energy due to the down-slope movement of the dense fluid is given by:  3 The time evolution of these three major energetic quantities and their sum is shown in Fig. 9. Molecular and turbulent diffusion and friction change the energy budget and have an overall effect of draining energy from the system.An other important point is the strong intermittency in time of the turbulent fluxes, which is clearly visible when considering the volume average of the vertical kinetic energy, shown in Fig. 4. Furthermore, the turbulent activity is strongest, when the mean velocity is small, showing that there is no connection between the average velocity (shear) and turbulent activity on (short) the inertial time scale.The kurtosis k 4 = ¯ w 4 / ¯ w 2 2 is around 5, larger than 3 (the value for a Gaussian distribution) which shows that the process is also intermittent in space.

Averages parallel to the floor (one dimensional results)
The previous subsection the governing equations for the volume averages were not closed as the space averaged friction force exerted by the bottom was missing.In this subsection I consider the results averaged in the x-y-plane as introduced in subsection 3.3, all this averages are then functions of the zcoordinate and time, only.The space averaged friction force exerted by the bottom can easily be determined by looking at the velocity averaged in the x-y-plane near the floor (u( z), v( z)) a small distance z (one grid point) above the floor, as the shear in the viscous sub-layer is given by: The velocities averaged in the x-y-plane at z = 0.5 m and 25 m above ground are shown in Fig. 5.A striking feature is that the oscillations near the ocean floor are highly reduced as compared to the geostrophic flow above.Near the bottom the water always flows down-ward, whereas the geostrophic part oscillates up and down.The instantaneous friction force exerted by the ocean floor is more related to the time average than to the instantaneous value of the geostrophic flow.In this mean sense, averaged over an inertial period, the direction of the friction force exerted by the floor is close to the 225 • to the geostrophic flow as predicted by laminar Ekman layer theory (see Fig. 8).Indeed, the friction force is almost constant with a decreasing inertial oscillation superposed of amplitude smaller than a quarter of the total friction force, whereas the velocities above are dominated by the inertial oscillations.For this reason the constant friction terms were introduced in Eqs. ( 28) and ( 29).An analytical solution for a stationary flow with an Ekman spiral, superposed by an inertial oscillation without Ekman dynamics is given in appendix A. This solution is very close to the problem discussed here.The frictional behaviour of the 3-D model is well reproduced by the 1-D model as can be verified in Fig. 5.
The averages in the x-y-plane of the x-and y-component of the velocity vector were already shown in Fig. 6.Our focus in this section will be on the turbulent fluxes so I start by considering the square of the z-component of the velocity vector (its average vanishes due to incompressibility).In Fig. 10 high values of w 2 in buffer layer with a positive skewness s 3 = w 3 / w 2 3/2 (skewness not shown) due to the dominance of strong upward ejections, called bursts, are observed.About two hours later a strong turbulent activity is observed at the interface, suggesting a propagation of turbulent activity from the bottom boundary layer to the interface.The skewness at the interface, however, has no definite sign.This indicates, that the turbulence is symmetric about the interface and that asymmetric transport processes across the interface, such as entrainment, are excluded.Entrainment is usually observed in gravity currents, when the Froude number exceeds unity, here it is below (see Table 2).Next, I consider the vertical advective transport of temperature (not shown).The values are mostly negative, indication that cold water is stirred upward.Near the interface (z = 200 × 0.5 m) at around t = 220 × 5 min large negative values are followed by positive values of similar magnitude.A closer inspection (not shown) reveals, that this is due to an increasing wave amplitude on the interface followed by a decline.This represents a stirring of water that is not followed by an irreversible mixing but by an "unstirring".In the present simulation this is most likely due to a flow that has not passed the mixing transition, as the Reynolds number is below 10 4 .
The diffusivity in the z-direction is given by: At the interface the large gradient leads to a small eddy diffusivity.This is somehow artificial, because when vertical stirring and mixing is considered we are more interested in the increase of potential energy (proportional to wT ) rather than mixing coefficients.Large values of the eddy viscosity due to small gradients, that is in almost homogeneous areas indicate a large mixing potential in areas where there is nothing to mix.
The other important quantity transported by the small scale turbulent fluxes is momentum.This transport is commonly modelled by an eddy viscosity.In its most general form the eddy viscosity is a fourth order tensor as explained in subsection 3.3.I start by considering its absolute value: It is considerable only in areas where the z-gradient of the x-and y-component of the velocity are small (not shown).All the criticism on the eddy diffusivity can be applied to the eddy viscosity.For the eddy viscosity we can, furthermore, consider the question about its isotropy, that is, if the transported momentum is aligned with the mean shear.A hypothesis is at the basis of the concept of eddy viscosity (see Wirth, 2010, for a detailed discussion).This question which is best considered by looking at which gives the cosine between the two vectors.Fig. 11 reveals, that there are large areas where this hypothesis is not met.
At the end of this subsection the appearance of Ekman layers at the floor and also at the interface is investigated.This problem is best considered by looking at the quantity which measures the veering or turning of the velocity vector averaged in the x-y-plane.The growth of the Ekman dynamics at the ocean floor is clearly visible in Fig. 12, with negative values indicating the clock-wise turning of the velocity vector with the distance from the ocean floor.A conspicuous feature is the absence of an Ekman dynamics at the interface.
When the viscous and diffusive effects are equal, the shear layer and the temperature gradient broaden at the same rate, and a geostrophic equilibrium is assured throughout the interface, leading to an absence of Ekman dynamics.This is discussed in detail in Wirth (2011).An interfacial Ekman layer has been observed by Umlauf and Arneborg (2009a,b) in a continuously forced gravity current in a canyon at one current section.This broadening reduces the shear at the interface.The magnitude of the velocity gradient at the interface is at least five fold smaller than the one at the bottom.This shows, that in the present regime friction is mostly due to the bottom friction.

Discussion and Conclusions
The spatial resolution of our model is around a thousand times coarser than the dissipation scale in the ocean.The explicit viscosity/diffusivity has to be increased by roughly the same factor as compared to the molecular values.In the presented calculations the turbulent fluxes are therefore largely dominated by those due to viscosity and diffusivity.The low Reynolds number implies that the mixing parameters and coherent structures identified in this study cannot easily be generalized to high-Reynolds number flows.This, however, does not prevent us from studying the turbulent regimes and fluxes Inertial oscillations occur whenever a gravity current adjusts to the topographic slope.I showed that the nature of the bottom friction is of a different nature for the geostrophic mean flow and the inertial oscillations.I demonstrated, that inertial oscillations have a determining influence on the turbulent fluxes, suggesting that parameterisations based on stationary flows might only capture part of the truth.The turbulence induced by inertial oscillations go through three different stages: laminar, roll structures and turbulence with upward bursts and streaks.The turbulent activity is maximal at the end of an inertial oscillation, when the average velocity is smallest.It is only during this stage that significant tur-bulent fluxes are observed.The data is used to demonstrate the anisotropy of the eddy viscosity tensor (see Wirth et al., 1995 andWirth, 2010) and the absence of an interfacial Ekman layer (Wirth, 2011).I furthermore show that when an attempt is made to parameterise the bottom friction in a model which does not resolve the bottom boundary layer the friction law depends on the time scale of the dynamics and is thus non-local in time.This shows the necessity of resolving the Ekman dynamics in ocean models, as already pointed out in Laanaia et al., 2010.When the vertical resolution at the bottom is fine enough the problem reduces to determining the anisotropic eddy viscosity tensor.Today's parameterisations which do not account for coherent structures and thus the anisotropy in the turbulent transport, can only be a first step towards solving the problem.
A deeper understanding and a better representation of bottom friction is key to the understanding and modelling of the ocean dynamics, not only when gravity currents are considered.The bottom friction also acts non-locally on the interior ocean dynamics through the Ekman pumping induced by the divergence in the horizontal Ekman transport.The surface stress of the winds transmitted through Ekman pumping to the interior ocean drives the ocean dynamics.The bottom stress transmitted through Ekman pumping to the interior ocean drains it (see Wirth, 2010, HDR).

Fig. 12 .
Fig.12.Sketch of the composition (sub-layers) of the turbulent planetary boundary layer.The thick horizontal line is the ocean floor, the distance from the ocean floor is given by the Z-coordinate.The extension of the different sublayers are given in the middle column the extension of the layers in terms of the physical parameters is given in the left column.The important processes in the intra and interaction of layers is given in the right column.

Fig. 2 .
Fig.2.Sketch of the composition (sub-layers) of the turbulent planetary boundary layer.The thick horizontal line is the ocean floor; the distance from the ocean floor is given by the Z-coordinate.The extension of the different sublayers are given in the middle column; the extension of the layers in terms of the physical parameters is given in the left column.The important processes in the intra and interaction of layers is given in the right column.

Fig. 4 .
Fig. 4. Left figure shows E 3 which is highly intermittent in time.Right figure gives the kurtosis k 4 with values that show an intermittence in space.

Fig. 5 .
Fig. 5. Time evolution u-component (black) and v-component (red) at 0.5 m (solid) and 25 m (dashed) above ground for the 1-D (linear) model (thin lines) and the averages of the 3-D (turbulent) model (thick lines).For better comparison the velocities at 0.5 m are multiplied by a factor 20. (the thickness of the laminar Ekman layer is 10 m).

Fig. 6 .
Fig. 6.Inertial oscillations dominate the up-slope velocity u (left) and cross-slope velocity v (right) in ms −1 .The interface is at z = 100 m.

Fig. 9 .
Fig.9.Time evolution of E 1 (dotted line), E 2 (dashed line), potential energy change E p (dashed dotted line) and the sum of the three (full line)

Fig. 11 .
Fig. 11.Cosine c of the angle between the shear in the z-direction of the velocity vector in the x-y-plane and the flux in the z-direction of the momentum in the x-y-plane (as defined in Eq. (38)).

Table 1 .
Parameters of numerical experiment.

Table 2 .
Non-dimensional parameters in experiments.