A 19-h time series of dissipation, stratification, and horizontal velocities has been obtained for a dense gravity current flowing into the Arkona Basin in the western Baltic Sea. The observations are compared with one-dimensional, quasi-steady theory, in which the gravity component in the flow direction is balanced by bottom friction, while that in the cross-flow direction is balanced by the Coriolis force. The observations deviate from the theory in that the bottom shear stress is more than 2 times as large as that required to balance the gravity. Several reasons for this discrepancy are discussed. A 1D turbulence model is also compared with the observations. Profiles of velocity, stratification, and dissipation rates generally show similar variations with depth as the observations, although the observed dissipation rates are somewhat larger than the modeled and the modeled transverse velocities are much larger than the observed. Subsequently, the model is used to investigate the variation of the entrainment parameter for a large range of Ekman and Froude numbers. Within the modeled parameter space, the entrainment parameter can be collapsed onto a curve that is an increasing function of both the Froude and the Ekman numbers. There is one puzzling result of the observations that differs from the model results and earlier observations: the observed entrainment rate increases dramatically during the observation period, where the Froude number decreases slightly. Some reasons for this increase are discussed.
Deep-water masses in ocean basins are often formed at local sources, where relatively narrow, dense gravity currents slide down the sloping ocean bottom to the deepest parts of the basin or intrude at their level of equilibrium density. Examples of such currents and field studies of these currents are the Denmark Strait overflow (Girton and Sanford 2003), the Faroe Bank Channel outflow (Duncan et al. 2003), the Mediterranean Sea outflow (Baringer and Price 1997), and the Red Sea outflow (Peters and Johns 2005). On the way downward, the currents entrain more or less of the surrounding ocean water into them, increasing the volume transport but reducing the density difference. Today, the entrainment process is still not well understood. An early entrainment model is that by Ellison and Turner (Turner 1986, hereinafter the ET model) based on laboratory experiments (Ellison and Turner 1959). The ET model predicts increasing entrainment rates for internal Froude numbers Fr that are larger than 1.12, and zero entrainment for smaller Fr numbers, where the internal Froude number is defined on the bulk layer velocity, thickness, and reduced gravity. The main problems with the ET model are (i) that many natural flows have Fr numbers smaller than unity and still entrain surrounding waters, and (ii) the ET model is based on laboratory experiment slopes between 12° and 90°, which are very seldom seen in ocean basins. In more recent papers, laboratory results for weaker slopes (Baines 2001) and rotation (Cenedese et al. 2004) have been presented. Alternative entrainment models for smaller Fr and realistic slopes have been proposed by Pedersen (1980) and Stigebrandt (1987). These models have shown to work well for small-scale deep-water renewals in fjords (Liungman et al. 2001; Arneborg et al. 2004a), and the latter model has been applied with success in Baltic Sea long-term circulation models (Stigebrandt 1987). The entrainment relations by ET, Cenedese et al. (2004) and Pedersen/Stigebrandt are compared with observational estimates in the Denmark Strait overflow, the Mediterranean outflow, and the Gullmar Fjord inflow in Fig. 1. It is seen that there is a quite large scatter of the observations, and that none of the entrainment relations do describe this scatter. However, all observational results are based on salinity and/or volume flux budgets, which means that both the entrainment rates and the Froude number are averages over large areas with the possibility of large variations within the area. In addition, the definition of the Froude number varies from case to case. To obtain a tighter parameterization of entrainment rates it is probably necessary to perform turbulence measurements in gravity currents. Johnson et al. (1994) presented dissipation measurements from the Mediterranean outflow, but mainly used them for an investigation of bottom and interfacial stresses. Thorpe scale analysis for estimating the diapycnal buoyancy flux has been used by Peters and Johns (2005) for the Red Sea overflow and by Fer et al. (2004) for the Storfjorden overflow. One problem with those measurements was that only one profile was obtained from each station, which causes an extreme data scatter that makes it difficult to separate spatial variations from local random fluctuations. Repeated casts at several stations in the Storfjorden overflow with a shear microstructure profiler have recently been reported by Fer (2006); however, no attempt has been made to compute entrainment rates.
In February 2004, a cruise within the Quantification of Water Mass Transformation Processes in the Arkona Sea (Quantas) project aimed at observing mixing in saline inflows to the Arkona Basin in the Baltic Sea (Fig. 2) obtained detailed observations of the temporal and spatial variations of such an inflow (Burchard et al. 2005b; Sellschopp et al. 2006). The main focus in the present work is on a 19-h time series of dissipation rates of turbulent kinetic energy, which was obtained during this cruise at a position north of the Krieger Shoal (Fig. 2). Our main objective is to investigate if a state-of-the-art turbulence closure model can describe the observed properties of the current, with main focus on the entrainment, and in given case use the model to extend the results to other parts of the parameter space.
The paper is organized as follows. In section 2, the theoretical framework is presented, including definitions of bulk layer parameters and a plume dynamics model. The third section is a presentation of data, with estimates of main dynamics and entrainment parameters. This is followed by a numerical turbulence modeling section, where the observations are placed into a broader parameter space. In section 5 the results are summarized and discussed.
a. Definition of layer variables and entrainment rate
We consider a layer of fluid situated on a sloping bottom, with slope, s = tan(α), and below a thick, passive layer with homogeneous density ρ0 (Fig. 3). The coordinate system is placed with the z axis pointing upward normal to the bottom, the y axis in the downhill direction, and the x axis along the depth contours. We define the layer depth D, the reduced gravity acceleration G′, and the layer velocities U and V by the relations
which are seen to reproduce the box-profile density and velocity values. With this definition the layer depth becomes 2 times the distance between the bottom and the center of gravity of the excess mass within the gravity current.
We use the following definition of the entrainment velocity, based on the parameters above,
In the stationary three-dimensional case the entrainment velocity is the import of ambient water into the gravity current per unit area and time. In the one-dimensional case it is the rate of layer thickening due to entrainment. The nondimensional entrainment parameter E is defined as
where Us = (U2 + V2)1/2 is the current speed.
In the present work we will concentrate on the one-dimensional case, where the conservation of mass equation can be written as
in which Kρ is a turbulent eddy diffusivity and ρ′ = ρ − ρ0. This is subject to the boundary conditions
In (7) it has been assumed that the density is a linear function of salinity and temperature, and that the turbulent diffusivities for salt and heat are equal. Integration of (7) with respect to z and multiplication with g/ρ0 yields
The relation between the entrainment rate and the diapycnal diffusion can be found by multiplying (7) with zg and integrating with respect to z. By use of the boundary conditions and (2) this can be written as
The right-hand side of (11) is the work against buoyancy forces by diffusion, while the left-hand side is the increase in potential energy.
In our coordinate system, the one-dimensional Boussinesq approximation momentum equations can be written as
where α is the bottom slope angle, Km is a turbulent eddy viscosity, and f ′ is the projection of the Coriolis parameter on the z direction [i.e., f ′ = f cos(α)]. The equations are subject to the boundary conditions
Integration with respect to z, yields
which by use of the boundary conditions (14) can be reduced to
where τbx and τby are the bottom shear stresses in the x and y directions. If we, finally, parameterize the bottom friction in terms of a drag coefficient Cd and introduce the entrainment parameter as defined in (6), (17) and (18) can also be written as
We do not expect the bottom drag coefficient Cd to be a constant, but we expect it to have values in the vicinity of 0.002–0.003. The main approximation in the bottom friction parameterization therefore is that it is directed oppositely to the main flow direction. In classical Ekman theory the bottom friction is turned 45° from the free flow direction, but in more realistic turbulence models this angle is much smaller. It is seen that the entrainment of fluid with zero momentum into the current acts as a drag term with drag coefficient E.
If we restrict ourselves to subcritical flows, the entrainment parameter, E, is much smaller than the bottom drag coefficient (Fig. 1) and can be neglected in (19) and (20). For time scales longer than the inertial period, we can also assume quasigeostrophic conditions, which means that we can neglect the accelerations, so that (19) and (20) reduce to
Last, we introduce the new coordinate (s, n), rotated by the angle β with respect to the (x, y) system, such that the s direction coincides with the flow direction (see Fig. 4). In this system, the quasigeostrophic, one-dimensional, layer momentum equations for subcritical flow can be written as
where we have used the relations U sinβ − V cosβ = 0 and U cosβ + V sinβ = Us, and where the slope angles, αs and αn, in the s and n directions are related to the total slope angle α and the angle between the x and the s axis β through
In nondimensional form (23) can be written as
is the Froude number. One therefore sees that for a subcritical flow, the Froude number is proportional to the square root of the slope in the flow direction. One also sees that the flow becomes critical when the slope in the flow direction is equal to Cd. For near-critical flows the entrainment parameter may, however, become comparable with the bottom drag coefficient, and the acceleration terms in (19) and (20) may also not be negligible, which means that (26) breaks down.
where K is the Ekman number,
which is basically a relation between the thickness of the Ekman layer and the thickness of the gravity current. For small Ekman numbers the Ekman layer is much thinner than the gravity current thickness, and the current is a geostrophic current that follows the depth contours. For large Ekman numbers the Ekman layer is larger than the gravity current thickness, and the current tends to flow straight down the hill. Then (23) and (24) can be combined to yield the following relation between the Froude number, the Ekman number, and the bottom slope:
So, the Froude number of the flow depends on the bottom slope divided by the bottom drag coefficient, and the Ekman number. The relation is shown in Fig. 5.
Alternatively, if the velocity is not known, it may be convenient to use the geostrophic Ekman number,
which by use of (26) and tanαs ≈ sinαs gives the Froude number,
which is also shown in Fig. 5. It is seen that the difference between using the real or the geostrophic Ekman number for predicting the Froude number is small. For KG < 1 the Froude number is very close to that predicted from geostrophic theory based on the total slope
which is independent of the bottom drag. For large geostrophic Ekman numbers the Froude number becomes independent of the Ekman number and approaches (sinα/Cd)1/2.
Last, we want to obtain an energy equation for later use in the entrainment parameterization. Multiplication of the 1D momentum equations (12) and (13) with their respective velocities, depth integration, and summation of the x and y components yields
where P is the production of turbulent kinetic energy
S = [(∂u/∂z)2 + (∂υ/∂z)2]1/2 is the vertical shear, and the coefficients γk and γb are defined by
Both γk and γb are equal to 1 for box profiles of velocity and density and close to 1 for realistic profiles.
In the quasi-steady limit the right-hand side of (36) is negligible. If we further assume weak slopes and subcritical flow, the entrainment parameter is much smaller than the bottom drag coefficient, meaning that (23) can be used to eliminate G′D from (36), which reduces to
Equation (40) expresses that the turbulence production is dominated by that caused by the bottom shear stress. Note that this does not imply that all turbulence used for the entrainment process is generated at the bottom. The most important turbulence for the entrainment is probably generated by shear instabilities within the interface.
Most of the produced turbulent kinetic energy P is lost to dissipation ɛ, while the rest is used to work against buoyancy by increasing the potential energy due to mixing. The work against buoyancy forces was found in (11) and can be written as
By introducing the bulk flux Richardson number Rf as the fraction of P that is used to work against buoyancy due to mixing,
or, using (26),
In the Pedersen (1980) and Stigebrandt (1987) model, (44) is used with a constant bulk flux Richardson number (≈0.035) to calculate the entrainment parameter. Our approach is different, in that we want to investigate how the bulk flux Richardson number depends on various parameters. This means that the only approximations behind (44) are the quasi-steady approximation and the neglected production of turbulence due to entrainment relative to the bottom turbulence production. The value of Rf may be an increasing function of K since the bottom boundary turbulence becomes fully constrained to homogeneous fluid when the layer is much thicker than the Ekman depth. On the other hand Rf may also be an increasing function of Fr as indicated by laboratory experiments on high-Fr number mixed layer turbulence.
We will investigate these functional relationships with the help of turbulence modeling in section 4, but first we will look into the observations from Arkona Basin in order to obtain a validation base for the model.
a. Location and instrumentation
The observations were conducted as part of a basin wide survey of the Arkona Sea to investigate the inflow and the subsequent spreading and mixing of salty North Sea water on its way into the Bornholm Basin. From 6 to 7 February 2004, during the most intensive period of the inflow, a 38-m deep station was established north of Kriegers Shoal (Fig. 2).
The local topography near the measuring site corresponds to a narrow, constricted channel with a downchannel slope of approximately αs = 0.0005 upslope of the constriction (Fig. 2). In a channel the flow is forced to flow in the channel direction. In that case, the 1D theory developed in the previous section is strictly not applicable since there will be cross channel thickness variations of the gravity current. But if we assume that 1D theory holds even for this case, with the slope, α, being the total slope of the interface rather than the bottom slope, and αs being the interface slope in the channel direction, then (26) can be used to predict the Froude number of the flow. Further complications may arise owing to the three-dimensional nature of the constriction near Kriegers Shoal. However, we will see in the following that the 1D approach, being the simplest applicable theory at our measuring site, can already explain many of the observed features. Some shortcomings of 1D theory have to be expected and will be discussed in detail below.
The experimental setup was designed to acquire currents, thermal microstructure, and turbulence simultaneously over the whole water column. During more than 19 h, currents were measured with a vessel-mounted acoustical Doppler current profiler (ADCP) and for the high-resolution measurements of thermohaline stratification, as well as microstructure shear and temperature, a loosely tethered free falling microstructure profiler (MST profiler) was used.
The transducer head of the 600-kHz version of RDI Workhorse Monitor ADCP was mounted in the sea chest at the level of the hull at 3-m depth. Every 1.4 s a current velocity profile was recorded with a vertical resolution of 1 m. Expecting that strong shears and high velocities would dominate the vertical current structure north of Kriegers Shoal, the water profiling mode 1 for “dynamic sea state” was selected. During the whole period at anchor, the ADCP operated in bottom tracking mode to remove ship movement around the anchor from the measured current profile data. A closer comparison between gyro heading and direction of the absolute current confirmed that this correction by means of the bottom-tracked ship motion worked perfectly. With a standard beam angle of 20° and a distance of 35 m between the transducer head at 3 m and the bottom at 38 m, the current measurements in the lowest 6% (=2.1 m) of the measured vertical profile can be contaminated due to bottom echoes from side lobes of the downward looking ADCP. Although the bottom around the anchorage is very flat, a certain variability of the depth due to movement of the ship around the anchor has to be considered regarding the thickness of the contaminated layer at the bottom. To eliminate disturbing effects near the transducer due to ringing effects and a possible contamination near the bottom, only current data in the depth range between 5 and 35 m are accepted so that the velocity structure in the surface layer as well as in the bottom boundary layer cannot be resolved properly.
The MST profiler (Prandke et al. 2000) is an early version of the MSS profiler (Lueck et al. 2002; Lass et al. 2003; Arneborg et al. 2004b) and is equipped with the same type of sensors. In the present setup it was equipped with a pair of airfoil shear probes and a fast microthermistor, which measured vertical profiles of horizontal velocity fluctuations and thermal microstructures at the same time. The airfoil shear probes are calibrated by the manufacturer for measuring the dissipation rate in a range between 10−2 to 10−11 W kg−1 with a resolution of 10−11 W kg−1 and a time constant of 4 ms. The time constant of the fast microthermistor is 10 ms and its resolution is 0.002°C. In addition, standard precision CTD sensors yielded the vertical density stratification. Disturbing vibrations induced by the cable and the housing were recorded by means of a horizontal acceleration sensor. The sampling frequency was 1024 Hz. The sinking velocity of the profiler was adjusted by additional weights to 0.78 m s−1. Even if sampling frequency and sinking velocity allow a vertical resolution regarding the horizontal velocity fluctuations of less than 1 mm, the resolution of the vertical velocity shear is limited by the size of the airfoil tip to roughly 5 mm. Further details of the microstructure sensors including the measuring principle, construction and calibration regarding the air foil shear probes, as well as the temperature microstructure sensor are explained in Prandke (2005).
b. Data analysis
While the ship rode at anchor north of Kriegers Shoal, a dataset of 496 microstructure and turbulence profiles was collected with a deployment rate of one profile every two to three minutes. The dissipation rate ɛ (W kg−1) of turbulent kinetic energy was calculated from the microstructure airfoil data using Taylor’s hypothesis of frozen turbulence and the usual assumption of isotropic turbulence in the relevant wavenumber band to obtain dissipation from the microstructure shear. For isotropic turbulence, the dissipation rate is related to the variance of vertical shear of the horizontal velocity through
where ν (m2 s−1) is the kinematic viscosity estimated from the actual vertical temperature profile. Prior to the data analysis, the raw microstructure profiles were carefully validated following the suggestions of Prandke et al. (2000). Subsequent to the conversion of raw data into physical values, spikes and values outside more than 2.7 times the standard deviation from the local mean were identified and eliminated by linear interpolation.
The pressure signal used for calculating the sinking velocity of the profiler can be possibly contaminated by surface wave induced fluctuations as well as by the tumbling profiler passing a strong vertical current shear layer. To remove these unwanted disturbances in the pressure signal, the pressure channel was low-pass filtered by means of a Butterworth filter.
The horizontal velocity shear and the profiling speed were calculated from the raw data of the shear sensor, the elapsed time, the low-pass filtered pressure signal and a mean density of 1010 kg m−3 according to Prandke et al. (2000). It has to be pointed out that uncertainties in the sinking velocity enter the computation of the physical shear with a power of 2 and the computation of the dissipation rate (45) with a power of 4.
The vertical shear variance is calculated in 0.5-dbar bins by spectral integration in the wavenumber interval from 2 to 30 cpm, which excludes high-frequency noise caused by, for example, probe vibrations. The lower limit is set by the probe length. Finally, the shear variance is corrected to compensate the unresolved part of the disturbed shear spectrum outside the chosen wavenumber interval using the universal Nasmyth spectrum (Prandke et al. 2000, Lass et al. 2003). For dissipation rates in the order of 10−5 W kg−1 the correction factor is about a factor of 3.
The hourly averaged dissipation rates, velocities in the main flow direction, and isopycnals are shown in Fig. 6. The gravity current is seen to be a relatively steady feature with a density surplus of around 7 kg m−3, a thickness of about 12 m, and velocities up to 0.6 m s−1. The maximum velocities are observed in a relatively sharp peak at the lower edge of the interface. Above the gravity current the dissipation rates are highly variable with a background level of 10−8 W kg−1 and patches of up to 10−5–10−4 W kg−1. Within the gravity current the dissipation rates are generally higher, ranging from 10−6 to 10−4 W kg−1. Generally, there is a minimum in dissipation rates about 11 m above the bottom, which is also the depth of maximum velocity. This is not surprising since the turbulence production is small at this level. Below the dissipation rate minimum, within the relatively homogeneous bottom layer, the dissipation rates increase toward the bottom and are relatively constant in time. Above the minimum, in the halocline, the dissipation rates are much more variable in time, with maximum dissipation rates ranging from 10−6 to 10−4 W kg−1.
The bulk layer variables defined in section 2 are shown in Fig. 7. The missing velocities close to the bottom are taken to be 0.83 times the lowest bin velocity, assuming a logarithmic velocity profile with roughness length z0 = 0.01 m (see below). The layer depth is relatively steady at about D = 12 m during most of the period, followed by a decrease toward 11 m in the end. The reduced gravitational acceleration increases more steadily from G′ = 0.065 m s−2 in the beginning to 0.083 m s−2 in the end. The eastward velocity is relatively steady around 0.48–0.50 m s−1 during most of the period and decreases toward 0.42 m s−1 in the end. The southward velocity is relatively steady around 0.12 to 0.14 m s−1. The main flow direction based on the layer velocities is 105°, that is, 15° south of eastward. Taken together, these results lead to a Froude number that varies around 0.55 during most of the period and decreases toward 0.45 in the end.
Profiles of density, velocities, gradient Richardson number
and dissipation rates, averaged over the whole 19-h period and over the first and last five hours are shown in Fig. 8, where N is the buoyancy frequency and S is the vertical shear of the horizontal velocities. The angle brackets in (46) denote ensemble averaging. The density shows a clear two-layer structure with a gradient layer between 25- and 28-m depth. The velocity has a strong shear layer in the same depth interval with a velocity maximum at about 29-m depth and decreasing velocities below that depth, as mentioned above. The transverse velocity has a somewhat unexpected shape with rightward currents closest to the bottom and leftward farther up, when looking in the flow direction. One would expect the opposite from Ekman theory. Probably, the difference is caused by transverse pressure gradients set up by transverse density gradients within the current. In fact, such density gradients are clearly seen in CTD cross sections north of Kriegers Shoal (Sellschopp et al. 2006). The gradient Richardson number is generally larger than one above the gravity current, except closest to the surface. At the top of the shear layer it decreases below 1 and obtains values near 0.1 at the observations closest to the bottom. The dissipation rates increase dramatically at the top of the shear layer when moving downward. They obtain a maximum in the center and decrease slightly at the bottom of the shear layer before increasing again toward the bottom. The local minimum at the bottom of the shear layer is located at the depth of maximum velocity, as mentioned above. The local maximum within the shear layer is much stronger during the end of the period than in the beginning, which is of large importance for the entrainment, as will be discussed below.
The upper quartile of the dissipation rates (Fig. 8d) is generally lower than the average value, which shows that a small number of high values are responsible for the average. This tendency is particularly large above the gravity current and within the shear layer. Closer to the bottom the average sinks below the upper quartile, indicating that turbulence bursts in the bottom boundary layer turbulence are not as extreme as they are in the shear layer and above the gravity current.
The gradient Richardson number is generally larger than the 0.25 needed to obtain shear instability. In the shear layer it is relatively constant around 0.7 although the dissipation rates indicate strong turbulence especially toward the end of the period. One possible explanation for this is that it is an effect of the finite (∼1 m) depth resolution for the velocities, and maybe even of the time averaging that may remove large short-term velocity gradients.
d. Bottom shear stress and momentum balance
The bottom shear stress can be obtained either by fitting a logarithmic profile to the velocity observations or from the dissipation rate measurements by assuming a logarithmic velocity profile, local balance between production and dissipation of turbulent kinetic energy close to the bottom, and a constant shear stress. With the latter assumptions the bottom shear stress can be written as
where κ = 0.4 is the von Kármán constant and z is the distance from the bottom. We will use the dissipation rate method to determine the bottom shear stress since the velocity measurements are not sufficiently close to the bottom to obtain a good logarithmic velocity profile fit.
The friction velocity u* = (τb/ρ)1/2 = Cd1/2Us calculated with the dissipation rate method, is shown in Fig. 9 as a function of the distance from the bottom. One sees that the friction velocity is reasonably constant around 0.03 m s−1 within 8 m from the bottom, with only a slight tendency of an increase toward the bottom. There is no difference between the estimates based on the whole period and those based on the first and last five hours. The quartiles show that there is a risk of 25% of obtaining an estimate that is more than 33% smaller than the average values (or 56% too small for the shear stress) if one uses a single dissipation profile to calculate the friction velocity.
The error bounds on dissipation rate calculations is often stated to be about a factor of 2, which seems to be reasonable even in our case where only one-third of the estimated dissipation spectrum is located within the measured wavenumber band, as mentioned above. This means that a reasonable estimate of the friction velocity based on these results is u* = 0.03 ± 0.007 m s−1. With a layer speed of Us = 0.49 m s−1 this gives a drag coefficient of Cd = 0.0037 ± 0.0015. This estimate is rather high relative to the typical values Cd = 0.002–0.003 often cited in the literature, although these values are included within the error band. As we will see below, it is also higher than the value obtained from large-scale dynamics based on the 1D theory in section 2. The large drag coefficient makes it appear likely that some of the assumptions involved in (47) are not fully satisfied. We will discuss some possible error sources in section 5.
If the assumptions behind (26) are correct, namely (i) that the flow is in a local quasi-steady balance between bottom friction and the gravitational forcing in the flow direction, (ii) that the interface slope in the flow direction is equal to the bottom slope in that direction, and (iii) that the bottom drag is much larger than the entrainment drag, then we can determine the bottom drag coefficient from the observed Froude number and from the bottom slope. The local bottom slope is not that easy to determine since the bottom levels out from a sinαs = 0.0005 slope upstream of the station to a more or less horizontal bottom downstream of the station. An upper estimate of the drag coefficient can be found by using a bottom slope, sinαs = 0.0005 and a Froude number, Fr = 0.5, which gives Cd = 0.002. This is a reasonable value, but only about half that found from the dissipation rate method above, and outside the relatively large error bounds used there. A drag coefficient Cd = 0.0037 and the observed mean Froude number Fr = 0.54 would require a bottom slope of sinαs = 0.0011, which is much steeper than the actual bottom slope. We are therefore in a position, where we without dissipation rate measurements would accept the 1D theory, where the gravitational component in the flow direction is balanced by bottom friction, whereas our dissipation measurements do not support this conclusion. In section 5 we will quantify some possible deviations from the 1D theory and discuss some error sources in the bottom shear stress estimate.
The entrainment rate can be found from the turbulent diapycnal diffusion of density, using (11). However, we need to first determine the turbulent diffusivity from either the dissipation rate data or from the temperature microstructure data. Since the noise levels of the temperature microstructure data are too large to use the Osborn and Cox method, we will concentrate on the dissipation data. An often-used parameterization of the diapycnal diffusivity of density is the Osborn (1980) model:
with Γ being the mixing efficiency with a maximum value of 0.2. This corresponds to a flux Richardson number of Rf = Γ (Γ +1)−1 = 0.17. Recent analyses of direct numerical simulations of shear generated turbulence (e.g., Shih et al. 2005) indicate that (48) is valid for intermediate values of the turbulence intensity parameter, 7 < ɛ/νN 2 < 100, while the diffusivity becomes a decreasing function of the turbulence intensity parameter for higher turbulence intensities. This also means that the mixing efficiency decreases toward zero for high turbulence intensities, which is natural since the fraction of work against buoyancy to dissipation must decrease toward zero in the limit of infinite dissipation or homogeneous fluid. There are also theoretical suggestions about the value Γ = 0.2 being too high for patchy turbulence in the interior water column (Arneborg 2002) since the mixed fluid loses potential energy in the collapse following the mixing event. Nevertheless, the interface of a gravity current is a typical case of turbulence generated by a constant shear, so, if the direct numerical simulations and laboratory experiments on mixing in shear-generated turbulence are relevant anywhere, then it is there. This means that we expect (48) to be valid within the interface of the gravity current with Γ = 0.2, if the turbulence intensity parameter is within the mentioned range 7–100, which is seen to be the case at least in the first part of our measurements (Fig. 10).
Below the interface there are at least two reasons why (48) cannot be expected to be valid. The main reason is that bottom-generated turbulence is characterized by a boundary layer length scale rather than an overturning length scale. Second, (48) breaks down also because the turbulence intensity parameter increases toward infinity when approaching the bottom since the dissipation increases while the stratification decreases toward zero. Instead, a reasonable assumption is that the relatively homogeneous bottom layer remains homogeneous, which means that the rate of density change is independent of depth within the bottom layer, which in turn means, using (7) and (8), that the buoyancy flux decreases linearly toward the bottom. Our approach for estimating the entrainment is therefore to use (48) for the diffusivity within the interface (−28 m < z < −25 m), where the turbulence intensity parameter is generally below 100, and a linear decrease of the buoyancy flux below this, within the relatively homogeneous bottom layer. The profiles of density diffusivity and buoyancy flux calculated this way are shown in Fig. 10 for the whole 19-h period and for the first and last five hours.
One sees that the diapycnal diffusivity within the interface increases dramatically from the first five hours to the last five hours of the 19-h period and that the buoyancy flux increases accordingly. The entrainment rates can be found from the integrated buoyancy flux through (11), and the results are given in Table 1. The entrainment rate also increases dramatically with a factor of 15 from the first 5 h to the last 5 h. The increased entrainment is opposite to what one would expect for a decreasing Froude number. The Ekman number is constant (Table 1) and can therefore not explain the change. We have no explanation for this unexpected behavior but will discuss some possible mechanisms in section 5.
4. One-dimensional turbulence modeling
a. Model description
The numerical model used here is based on the one-dimensional Reynolds-averaged equations for the along-slope and cross-slope velocity components, (12) and (13), and for density, (7), where the turbulent fluxes have been expressed in terms of turbulent diffusivities. The turbulent diffusivities of momentum Km and mass Kρ are calculated by means of the Prandtl–Kolmogorov relation, which in our case reads
where k is the turbulent kinetic energy per unit mass and ɛ is the dissipation rate per unit mass, such that ω denotes an inverse turbulence time scale; cμ and c′μ are nondimensional stability functions, which are discussed below. The one-dimensional transport equations for k and ω read
where σk and σω are turbulent Schmidt numbers used for the gradient-transport models and Cω1, Cω2, and Cω3 are model parameters, discussed in detail by Wilcox (1988) and Umlauf et al. (2003). The stability functions result from the algebraic second moment closure (SMC), discussed by Umlauf and Burchard (2005), and depend on the nondimensional shear number αS = S2k2/ɛ2 and the nondimensional buoyancy number αN = N 2k2/ɛ2. For our explicit algebraic SMC, we adopted the coefficients suggested by Cheng et al. (2002).
Equation (51) is simply the transport equation of turbulent kinetic energy with the well-known shear and buoyancy production terms and the dissipation rate, on the right-hand side, and a downgradient model for the turbulent transport terms, on the left-hand side. The physics contained in (52) can be understood by using (49) and (50) in order to replace the diffusivities and the dissipation rate on the right-hand side. Neglecting the transport term, (52) can then be seen to describe the relaxation of the inverse time scale, ω, toward the weighted average of S and N, which are the inverse time scales set by the shear and the buoyancy, respectively. For the given form of the stability functions, it can further be shown that for weak stratification ω ∝ S, whereas for strong stratification, ω ∝ N. This corresponds to the expected behavior in stratified shear flows (see Baumert 2005; Shih et al. 2000).
The model has been compared with direct numerical simulation (DNS) data on stratified shear flows. In agreement with DNS by Shih et al. (2000), stationary turbulence is predicted for a so-called stationary Richardson number of Ri = 0.25. For all values of Ri up to this value, the key parameters of stratified turbulence, most notably the turbulent Froude number and the mixing efficiency, were shown to be in good agreement with DNS data from various authors (see Umlauf 2005). For Ri = 0.25, the predicted turbulent Froude number is close to 1 and the mixing efficiency is close to Γ = 0.2, as expected from DNS and laboratory data. For Ri > 0.25 the model predicts exponentially decaying turbulence if no turbulence sources from turbulent transport are available. Near a solid boundary the model is in agreement with boundary layer similarity theory.
For the surface and bottom boundary conditions, flux conditions derived from the classical boundary layer theory have been used (see Umlauf and Burchard 2003). The spatial discretization is based on a staggered grid, where the mean flow properties are computed in the grid cell centers and the turbulent quantities at the interfaces (see Burchard et al. 2005a).
b. Validation against observations
To compare with the observations, the model was run for 72 h on a slope tanα = 0.0013 with initial box profile salinity and velocity distributions with D = 10 m, G′ = 0.089 m s−2 (ΔS = 10, ΔT = 0), and a geostrophically balanced velocity. The bottom roughness is set to z0 = 0.01 m, which is quite rough, but necessary in order to obtain a drag coefficient in the vicinity of the observed Cd = 0.0037. The total slope is chosen so that it is equal to the combination of the observed transverse interface tilt (tanαn = 0.0008) and the streamwise bottom slope (tanαs = 0.0011) needed to explain the observed bottom shear stress. The other initial parameters are chosen so that DG′ is similar to the observations but leaving space for the thickening and reduced gravity reduction during the initial adjustment toward quasi-steady conditions.
After 15 h the integrated layer quantities corresponded well with the observed values as seen in Table 2. Here, and in the following intercomparison, the observed quantities are ensemble-averaged quantities for the whole 19-h period observational campaign.
The profiles of density anomaly, velocities, gradient Richardson number and dissipation rates are compared in Fig. 11. The modeled salinity profile is very similar to the observed. The streamwise velocities correspond reasonably, although the observed velocity maximum is situated closer to the bottom than the modeled, the observed interface shear is smaller than the modeled, and the observed velocity falls off quicker than the modeled below the maximum. The modeled transverse velocity profile is much stronger than the observed, although there is some correspondence in structure within the interface. Closer to the bottom, the modeled and observed transverse velocities are oppositely directed with the modeled being directed toward the left when looking in the flow direction, as expected from Ekman veering. We have no explanation for the differences between the observed and modeled transverse currents, except that they may be the result of 2D or 3D effects that are not included in the 1D modeling results. We will look further into this in future work.
The observed and modeled gradient Richardson numbers show similar features with a local minimum within the interface, a local maximum close to the velocity maximum (where the shear is small), and a decrease toward zero below that. The modeled gradient Richardson number has the expected critical value Ri = 0.25 in the interface, while the observed value is closer to one. The large observed value may be an effect of too low vertical resolution of the velocity profile, but it may also be an effect of intermittency in the real data that is absent in the model results. However, the resolution of the dataset is too coarse to show such intermittency.
The observed and modeled dissipation rates show a similar pattern with a local maximum within the interface, a local minimum close to point of minimum Richardson number where the production is small, and a sharp increase toward the bottom. However, the observed dissipation rates are almost a factor of 3 larger than the modeled within the interface, but the difference decreases toward the bottom. An almost perfect correspondence could be obtained by excluding the last hours of the dissipation time series from the average, which illustrates the uncertainty of the dissipation rate profile. We are therefore relatively satisfied with the model performance.
Although the observed dissipation rates within the interface are so much larger than the modeled, the modeled entrainment rate is only 15% smaller than that estimated from the dissipation measurements (Table 2). The reason for this is that most of the integrated buoyancy flux takes place below the interface, where the estimate from the observations is very similar to that from the model.
The similarities between the observations and modeling results make us confident that it is interesting to extend the model results to other parameter values. However, we must keep in mind that assumptions behind the 1D theory (section 3) and therefore also the 1D model differ from reality in at least a couple of aspects, namely (i) that a more than 2 times as steep slope as the real bottom slope is needed in order to obtain the observed bottom shear stress and (ii) that the observations show very large variations on time scales of several hours without any strong variations in the velocity and density fields. These features will be discussed further in section 5.
c. Multiparameter runs
Now we will present model results for a number of cases with initial conditions listed in Table 3, which have been chosen to cover an as large area of the Froude and Ekman number space as possible. All of the runs were initialized with box profile salinity and velocity distributions, with the initial velocities being in geostrophic balance. The bottom roughness was set to z0 = 0.01 m. The water depth was 500 m for runs 1–5 and 100 m for runs 6–12. The model was run for 10 days with a time step of 10 s and a stretched grid with finest resolution closest to the bottom and 200 grid cells.
The model results are plotted in Fig. 12 in terms of Froude number versus Ekman number with 6-h intervals between the points. The corresponding theoretical solutions (30) are plotted as continuous curves. Each of the runs start in the high-Fr and high-K end of the curve, which corresponds to thin and fast currents, and evolve downward to the left as the gravity currents become thicker and slower due to entrainment. They all adjust toward the theoretical solution within some inertial periods. The initial differences between numerical and theoretical model results are caused by inertial oscillations, which are not covered by the theory. Runs 1, 2, 3, and 6 are supercritical parts or all of the time but still follow the theoretical curve, which must mean that the modeled entrainment drag is negligible relative to the bottom drag even for these cases.
The entrainment parameter is shown as a function of the Froude number in Fig. 13. It is seen that the entrainment rates generally cluster around the Pedersen (1980) and Stigebrandt (1987) curve. For a given run, that is, a given bottom slope, the Froude number dependency (curve slope) is similar to that found by Cenedese et al. (2004). The clustering around the Pedersen and Stigebrandt curve means that the bulk flux Richardson number generally is in the vicinity of 0.035. For a single run it is, however, an increasing function of the Froude number, as shown in Fig. 14, with an overall minimum value of about 0.01 and a maximum between 0.05 and 0.06. It is worth remembering that these small bulk flux Richardson numbers occur even though the local mixing efficiency within the interface is Γ = 0.2, corresponding to a local flux Richardson number of 0.17. The reason for the much smaller bulk value is that the major dissipation happens close to the bottom, whereas the main buoyancy flux happens in the interface, where the dissipation is relatively small. Simulations with smaller roughness length z0 = 0.001 m (not shown) generally gave smaller entrainment rates and bulk flux Richardson numbers below 0.03.
A large part of the scatter in Fig. 14 is caused by a dependency of the Ekman number, so in order to obtain a better collapse of data we fitted a function of type
to the data. In log–log scale this becomes a linear function in a, b and log(A). The best fit was obtained for the function
which is shown in Fig. 15. The six outliers on the low side of the dashed line are all obtained within the first inertial period, where the current is still adjusting toward the quasi-steady state. The otherwise nice fit in Fig. 15 shows that the flux Richardson number is an increasing function of both the Ekman number and the Froude number. An increase with Froude number could be expected from the Ellison and Turner (1959) and the Cenedese et al. (2004) results. An increase with Ekman number is also to be expected intuitively since a small Ekman number means that the most energetic turbulence is restricted to a thin boundary layer close to the bottom where it cannot mix fluid through the interface, whereas a large Ekman number means that the bottom-generated turbulence extends out to the interface so that the fluid that is mixed through the interface can also be mixed all the way to the bottom. This result is, however, not supported by our observations, where a dramatic increase of the entrainment rate was observed with a decreasing Froude number and an almost constant Ekman number. This will be discussed further in section 5.
5. Summary and discussion
A 19-h time series of dissipation, stratification and horizontal velocities was obtained for a dense gravity current flowing into the Arkona Basin in the western Baltic Sea. The observations were compared with one-dimensional, quasi-stationary theory, where the gravitational component in the flow direction is predicted to be balanced by bottom friction, while the Coriolis force is balanced by a cross-flow pressure gradient. The cross-flow geostrophic balance was seen to work well, whereas the streamwise balance between gravity and bottom friction was more doubtful. The bottom drag coefficient estimated from large-scale dynamics was reasonable, whereas that estimated from dissipation rates was almost twice as large. Possible reasons for this discrepancy will be discussed below. Another feature of the observations was a more than factor of 10 increase of the entrainment rate whereas the Froude number was slightly decreasing and the Ekman number practically constant. Such a variation is in contradiction to earlier observations, where entrainment rates tend to increase with increasing Froude numbers. On the other hand, this is the first long time series of dissipation rates from one position within a gravity current.
The observations were also compared with 1D turbulence model results. When using a slope that was larger than the real slope and equal to that predicted from 1D theory for the observed bottom shear stress, the correspondence was good for stratification and streamwise velocities. The modeled dissipation rates were somewhat too small relative to the observed, but the overall structure compared well, and the entrainment rates also compared well. The Ekman veering predicted by the model was not seen in the observations, which has probably to do with secondary transverse pressure gradients that are not included in the 1D model. The observed gradient Richardson number in the interface was larger than the modeled, which was the typical 0.25. The difference may have been caused by a too low vertical resolution of the observed velocities.
Subsequently, the model was run for a number of different cases. The modeled entrainment rates could be collapsed into one increasing function of Froude number and Ekman number, and were seen to correspond well to the low Froude number part of earlier reported entrainment rates. The bulk flux Richardson number were also seen to be an increasing function of Froude number and Ekman number and were generally situated in the interval 0.01 < Rf < 0.06. The reason for these small values is that most of the dissipation takes place within a very small distance from the bottom where the fluid is more or less homogeneous.
a. Error sources in the bottom shear stress estimate
The main result of the 1D theory, which was also seen in the model results, is that the gravitational component in the flow direction is balanced by the bottom shear stress. This important feature was not supported by the observations, where the bottom shear stress estimated from dissipation rates was more than a factor of 2 larger than the gravitational component estimated from the bottom slope. One possible reason for this discrepancy is that the bottom shear stress estimate is too large. What speaks in favor of this explanation is that the bottom drag coefficient is much larger than what is usually reported from bottom boundary layers. One problem with the dissipation rate method may be enhanced dissipation levels caused by turbulence production from convective overturns in the bottom boundary layer (Moum et al. 2004; Lorke et al. 2005). In a gravity current an upslope density gradient exists due to entrainment, and the large velocity shear near the bottom therefore causes dense water to ride over less dense water, causing boundary layer convection and an additional turbulence source near the bottom. According to (47), this process will lead to enhanced apparent values of the bottom shear.
b. Error sources in the 1D theory
Alternatively, the observed bottom shear stress is correct and the 1D balance is wrong. Some possible reasons for this discrepancy are (taking the missing terms in the momentum equation one by one) (i) that the velocity time derivative is important, (ii) that the fluid is decelerating in the flow direction, (iii) that the interface slope is steeper than the bottom slope, and (iv) that there are lateral shear terms accelerating the fluid.
(i) The gravity current is decelerating from about Us = 0.51 to 0.44 m s−1 during the 19-h period. The bottom shear stress needed to accomplish this deceleration is about 0.1 N m−2, one order of magnitude smaller than the estimated bottom shear stress τb = 0.9 N m−2.
(ii) and (iii) The channel is widening downstream of the station (Fig. 2), which for a subcritical current could lead to a deceleration and thickening in the flow direction, which would balance each other in the absence of bottom friction. However, in a bottom-friction dominated current it may rather be the friction that balances the deceleration, leaving the thickness constant or even thinning. Assuming the thickness to be constant, the advective term caused by channel widening is of order Us2W−1∂W/∂s, where W is the gravity current width. The bottom shear stress needed to balance such a term is ρDUs2W−1∂W/∂s. In our case, a contribution of about one-half of the observed bottom shear stress would require a twofold increase of the gravity current width within about 3 km.
(iv) If we happened to hit a local spot with large bottom roughness, the surrounding fluid would tend to accelerate the local fluid with lateral shear stresses (horizontal eddies). To investigate this, we need dissipation measurements at more than one point. Such measurements have actually been done in 2005 and will be presented in a future publication.
It is possible that the missing balance between the gravitation component and bottom friction is caused by a combination of the aforementioned points. A more important question is whether it makes sense to compare a 1D numerical turbulence model forced by gravity with observations where other types of forcing are present. What speaks in favor of this is that all of the additional forces and accelerations mentioned above are relatively constant over the depth of the gravity current, just as the gravity forcing.
c. Entrainment rates
Another puzzling result of the observation is the more than factor of 10 increase in dissipation and entrainment rates during the 19-h period in combination with a slightly decreasing Froude number. Such a variation is not predicted by the model and is in contradiction with previous entrainment studies (Fig. 1). One explanation for this result is that interfacial waves generated at some other position by some unknown process reach our location at the end of the series and cause the increased dissipation. Another explanation is the change in transverse structure seen in Fig. 8. The transverse shear at the level of maximum velocity is seen to be much larger during the last five hours than during the first five hours, and consequently the Richardson number maximum at this level is absent during the last five hours, enabling the production of turbulence and transport of buoyancy through that level. We will investigate this process further in future work with 2D modeling of a gravity current flow in a channel. We cannot explain this result at present, but it needs attention in the future since it has a large impact on our hope to investigate gravity currents with reasonably short dissipation time series. It also shows that previous dissipation data in gravity currents with only one profile at each station should be taken with some caution.
The entrainment results obtained by the model (Fig. 13) are smaller than many of the earlier observations presented in Fig. 1. One reason for this may be that our model results are not valid at the high Froude numbers where these high entrainment rates are found. Another reason may be that many of these observations are based on budgets rather than point measurements. For example, the entrainment at gravity current fronts seems not to be negligible relative to the entrainment in the main current (e.g., Lass and Mohrholz 2003), and probably the entrainment at the sides of a current may also be larger than that away from the sides. The latter effect will be investigated in a continuation of the present work.
The work for this study was carried out in the framework of the international QuantAS Consortium (Quantification of water mass transformation processes in the Arkona Sea), which is partially funded by the QuantAS-Off project by the German Federal Ministry of Environment, Nature Conservation and Nuclear Safety (BMU). Lars Arneborg was supported by the Swedish Research Council. We are grateful for the support and flexibility of the captain and crew on board the R/V Helmsand. Last, we deeply acknowledge the professionalism and enthusiasm of, and the pleasant time together with, the other cruise members: Jürgen Sellschopp, Michaela Knoll, Frank Gerdes, Jens Benecke, Horst Schock, Dietmar Rüss, Günter Plüschke, Hans-Ulrich Lass, Volker Mohrholz, and Holger Prandke, without whom the present work would not have been possible. The numerical studies have been performed with a modified version of the General Ocean Turbulence Model (see www.gotm.net). The map is plotted using the m_map toolbox developed by Rich Pawlowicz. We also thank two anonymous reviewers for useful comments on an earlier version of this manuscript.
Corresponding author address: Lars Arneborg, Göteborg University, Box 460, 405 30 Göteborg, Sweden. Email: firstname.lastname@example.org