On the sensitivity of internal waves and mixing at a sill to horizontal resolution and
sub-grid scale mixing
Jarle Berntsen
a,∗
aDepartment of Mathematics, University of Bergen, Johannes Brunsgate 12, N-5008 Bergen, Norway
Jiuxing Xing
bAlan M. Davies
bbProudman Oceanographic Laboratory, Joseph Proudman Building, 6 Brownlow Street, Liverpool L3 5DA, UK
Abstract
A non-hydrostatic terrain following model in cross sectional form is applied to study the generation and propagation of tidally forced internal waves near a sill in an idealized loch. On inflow internal waves are generated behind the sill, and they propagate from the sill. There is a transfer of energy from the barotropic tide to internal waves and then to irreversible mixing. The range of length scales involved goes from the scale of the forced tide to scales associated with wave breaking.
This wide range of scales makes it very difficult to resolve all relevant length scales with a numerical model using a uniform finite difference grid. The sensitivity of the numerical results to the grid resolution and the parameterization of sub-grid scale mixing is investigated. Initial calculations are performed with constant values of viscosity and diffusivity, and subsequently they are repeated using large eddy simulations and including Richardson number dependant vertical mixing.
In the studies with sub-critical flow over the sill, the results are robust to changes in the grid size. With stronger forcing, the flow over the sill becomes super-critical during maximum inflow, and the frequency and the nature of the internal waves generated in the lee of the sill becomes grid size dependant. For the super-critical case, the solution is sensitive to the sub-grid closure scheme. Calculations show that with constant and small values of vertical viscosity and diffusivity, strong internal waves are generated on inflow in the lee of the sill. When using Richardson num- ber dependant vertical viscosity and diffusivity, more of the tidal energy goes into vertical irreversible mixing, and less into internal waves.
Key words: Non-hydrostatic, ocean modelling, tides, sills, internal waves, grid resolution, mixing
∗ Corresponding author
Email addresses: [email protected] Tlf.:47-55-584854 Fax.:
47-55-589672(Jarle Berntsen), [email protected](Jiuxing Xing),[email protected] (Alan M. Davies).
1 Introduction
The role of turbulent mixing on the ocean circulation has been discussed in many recent papers, see for instance Munk and Wunsch (1998); Samelson (1998); Spall (2001); Wunsch and Ferrari (2004). How energy due to the major forcing mecha- nisms namely wind and tide is converted into small scale mixing and how this feeds back into the large scale circulation is a major problem in oceanography. It has been suggested that internal waves interacting with topography play an important role in this mixing process, see Polzin et al. (1997); Kunze and Llewellyn Smith (2004). Many studies focus on the processes along the continental slopes, see Wun- sch (1970); Armi (1978); Legg and Adcroft (2003) and the discussion in Thorpe (2005). The process by which large scale energy of wind and tidal energy is con- verted into mixing involves the cascade of energy through a spectrum of time and space scales. This makes it very challenging to measure or model all details of internal wave generation, propagation, breaking, and mixing. Many notable studies of these processes are summarised in Thorpe (2005) and Vlasenko et al.
(2005).
Internal waves are also found in fjords. Fjords are smaller and shallower than the open ocean, and many of the important physical phenomena of the deep ocean are also found in fjord systems. Therefore, fjords may be used as laboratories to investigate baroclinic processes, and in particular processes near sills where en- ergy is transferred by non-linear effects from large to small scales and eventually mixing, see for instance Stigebrandt (1976); Stigebrandt and Aure (1989); Pars- mar and Stigebrandt (1997); Stigebrandt (1999); Cushman-Roisin et al. (1989);
Tverberg et al. (1991); Farmer and Freeland (1983); Farmer and Armi (1999);
Cummins et al. (2003); Stacey (1985); Stacey and Gratton (2001); Klymak and Gregg (2001, 2003, 2004); Skogseth et al. (2005, 2007); Fer et al. (2003); Fer (2006); Inall et al. (2004, 2005); Eliassen et al. (2001); Vlasenko et al. (2002).
To gain some insight into how barotropic tidal energy is converted into internal waves and mixing Xing and Davies (2006a) (hereafter XD06a) applied an idealised cross sectional model of Loch Etive (a region of recent intensive measurements, see Inall et al. (2004, 2005)) to study the dynamics near the sill. Subsequently they studied the influence of stratification and topography upon internal wave spectra (Xing and Davies (2006b) hereafter XD06b), and the influence of stratification and tidal forcing upon in particular how stratification in the region near the upper part of the sill influenced mixing (Davies and Xing (2007)). The importance of non-hydrostatic pressure effects is demonstrated in Xing and Davies (2007).
The results reported in XD06a,b were produced with the MITgcm, see Marshall et al. (1997a,b); Adcroft et al. (1999), using a horizontal grid size of 10 m near
the sill and gradually coarser resolution away from it. A grid size of 10 m is very small if we relate it to grid sizes often used in numerical ocean modelling, see for instance Haidvogel and Beckmann (1999). However, both observations and numerical model results for Knight Inlet show processes around this inlet that may not be well resolved with a grid size of 10 m, see for instance Farmer and Armi (1999); Cummins et al. (2003). The topography at Knight Inlet is different from the topography considered for the idealised Loch Etive in XD06a.
The corresponding forcing and stratification are also different. However, the flow near the sill described in XD06a is very strong, and one may expect small scale processes also in this loch. Therefore, we want to investigate the sensitivity of numerical results for the model area in XD06a to the grid resolution and the sub-grid scale parameterisation. Using the XD06a model domain it will be very computationally demanding to refine the grid further. The approach chosen in the present paper is, therefore, to apply a sequence of equidistant grid sizes ranging from 100 m to 12.5 m to examine convergence. The model applied in the present studies is a non-hydrostatic cross sectionalσ-coordinate model, see Berntsen et al.
(2006). By using aσ-coordinate model, we hope to achieve a better representation of the processes in the bottom boundary layer, see for instance Haidvogel and Beckmann (1999), than XD06a,b who used a z-coordinate model with shaved cells in the bed region.
The experiments in XD06a,b were performed with constant values of viscosity and diffusivity, both horizontally and vertically. Since there may be unresolved small scale processes in such studies of stratified flow over sills, one may expect the numerical results to be sensitive to the choice of sub-grid scale closure. This sensitivity is also investigated in the present study.
In section 2 the numerical experiments are described. Some results from experi- ments where the maximum velocities over the sill are less than the internal wave speed are given in section 3. In section 4 corresponding results for the case with super-critical flow over the sill are given. The sensitivity of the results to the parameterisation of sub-grid scale effects is discussed in section 5. A summary and a general discussion is given in section 6.
2 The numerical experiments
2.1 The model and experimental setup
Theσ-coordinate ocean model applied in the present studies is a two dimensional, (x, z), version of the model described in Berntsen (2000) where x and z are the horizontal and vertical Cartesian coordinates respectively. The model is available from www.math.uib.no/BOM/. The variables are discretized on a C-grid. In the vertical, the standardσ-transformation,σ= H+ηz−η, whereηis the surface elevation, and H the bottom depth, is applied. For advection of momentum and density a Total Variance Diminishing (TVD)-scheme with a superbee limiter described in Yang and Przekwas (1992) is applied in the present studies. The standard second order Princeton Ocean Model (POM) method is applied to estimate the internal pressure gradients (Blumberg and Mellor (1987); Mellor (1996)). The model is mode split with a method similar to the splitting described in Berntsen et al.
(1981) and Kowalik and Murty (1993).
Using a splitting technique, the non-hydrostatic pressure may be computed by in- serting the expressions for the non-hydrostatic velocity corrections into the equa- tion of continuity (Kanarska and Maderich (2003); Heggelund et al. (2004)). Due to theσ-transformation, additional terms appear in the pressure equations, which complicates the computations. An alternative method, that has been adopted in the present study, is to view the non-hydrostatic pressure directly as P(x, σ, t), wheretis time, or a pressure due to convergence or divergence in theσ-coordinate system (Berntsen and Furnes; 2005).
The time steps are performed with a predictor-corrector method where the leapfrog method is used as predictor and the θ-method, where θ is the degree of implicit- ness, is used as corrector. For stability θ must be in the range 0.5≤ θ ≤1, and in the present experimentsθ = 1.0 (i.e. fully implicit) (Haidvogel and Beckmann;
1999; Casulli; 1999).
The model has recently been applied to study lock release gravity currents and the propagation of solitary waves up an incline (Berntsen et al. (2006)) and the results are related to measurements from laboratory experiments and to corre- sponding results from the MITgcm (Marshall et al.; 1997a,b; Adcroft et al.; 1999).
In Berntsen et al. (2006) more details about the estimation of the non-hydrostatic pressure are given.
The set up of the numerical experiments is very similar to that described in XD06a. The topography and the initial stratification for the first set of exper-
iments are given in Fig. 1. The initial buoyancy frequency N is 0.01 s−1 in all experiments. The depth profileH(x) in meters is specified according to
H(x) =
−50 + 1+(x/W35sill)∗∗4 , x <0
−100 + 1+(x/W85
sill)∗∗4 , x >0
assuming x = 0 m at the top of the sill and Wsill is the half width of the sill. In the first set of experiments Wsill is chosen to be 1000 m, and in the second set it is set to 500 m. The Coriolis frequency f = 1.2×10−4 s−1.
Figure 1
There is no flow through the sea bed or the closed right end of the loch. Initially the water elevation is −η0, and there is no flow. The M2 tide is forced into the loch by specifying a velocity in the x-direction u as
u=utide−c(η−ηtide)/HI,∀z
at the left boundary. In the above equation HI is the depth at the inflow point, c=√
gHI,η is the water elevation just inside the left open boundary,ηtide is the forced tidal water elevation, and utide is the corresponding forced tidal velocity.
The tidal forcing is computed from a standing wave solution given in Gill (1982) p. 113
ηtide =−η0cos(kLx) cos(ωM2t), and
utide = cη0
HI
sin(kLx) sin(ωM2t),
whereωM2 = 0.000140526 rad s−1, k =ωM2/c, and Lx the length of the channel (Lx = 21500 m). The value of η0 is 0.3 m in the first set of experiments giving sub-critical flow over the sill at maximum inflow. In the second set of experiments η0= 3.8 m giving super-critical flow over the sill at maximum inflow as in XD06a.
The standing wave solution is strictly for a channel with constant depth. However, by forcing with the above open boundary condition, the tidal signal propagate with very little disturbances into and out of the domain. For the temperature a Neumann boundary condition is applied at the open boundary.
In the vertical 100 equidistant σ-layers are used in all experiments. Horizontally the grids applied are equidistant and the experiments are run with horizontal grid sizes ∆x equal to 100 m, 50 m, 25 m, and 12.5 m. In the experiments with ∆x = 12.5 m the internal time step used is 0.56 s, and 30 external time steps are used for each internal step. The time steps are increased proportionally with ∆x in the experiments with coarser grid resolutions.
The experiments are performed with constant values of viscosity and diffusivity as
in XD06a, that is: horizontal viscosity Ah = 10−1m2s−1, vertical viscosity AV = 10−3m2s−1, horizontal diffusivityKh= 10−7m2s−1, and vertical diffusivityKV = 10−7m2s−1. With the present grid sizes all relevant length scales are not resolved and the effects to be parameterised by eddy diffusivity/viscosity will depend on the grid resolution. Therefore, the experiments are also performed with values of eddy viscosity and eddy diffusivity that depend on the grid size and that vary in time and space. For the horizontal viscosity a Smagorinsky (1963) type formulation is applied, and the viscosities are computed from
Ah =CM∆x2 [(∂u
∂x)2 +1 2 (∂v
∂x)2]12, (1) where CM is a non-dimensionless viscosity parameter andv is the velocity com- ponent in the cross loch direction. Statement (1) is a 2D version of the formula- tion applied in the Princeton Ocean Model (Blumberg and Mellor (1987); Mellor (1996)). Values ofCM in the range 0.05 to 0.20 are suggested in Mellor (1996) and in Haidvogel and Beckmann (1999). In the present studiesCM is chosen to be 0.2.
To compute the vertical eddy viscosities and diffusivities, the Mellor and Yamada (1982) 2 1/2 level model with the modifications due to Galperin et al. (1988) are used. The minimum allowed value of AV is set to 10−5m2s−1 and the minimum allowed value ofKV is 10−7m2s−1. The experiments with the Smagorinsky model and the Mellor and Yamada model are denoted as large eddy simulations (LES).
2.2 Model diagnostics
In the numerical experiments the kinetic energy Ekin is computed in each time step according to
Ekin(t) =
Z 15000m
−5000m
Z η(x,t)
H(x) 0.5ρ(x, z, t)(u(x, z, t)2+v(x, z, t)2+w(x, z, t)2)dzdx , (2) where ρ is the density computed from the temperature and a linear equation of state and u, v, and w are the velocity components in x, y, and z direction respectively. Defining
(u(x, t), v(x, t)) = 1 D(x, t)
Z η(x,t)
H(x) (u(x, z, t), v(x, z, t)) dz , (3) where D(x, t) = H(x) + η(x, t), the barotropic energy Ekin2D is computed from
Ekin2D(t) =
Z 15000mZ η(x,t)
0.5ρ(x, z, t)(u(x, t)2+v(x, t)2)dzdx . (4)
The baroclinic kinetic energy Ekin−bclin is computed from
Ekin−bclin(t) = Ekin(t) − Ekin2D(t). (5) At each time step the maximum velocity componentUmaxinx-direction is recorded according to
Umax(t) =max∀x,zu(x, z, t). (6) To investigate the sensitivity of the stability of the water masses in the area around the sill to the parameters involved in the various calculations, the area of the water masses with Richardson numbers (Ri) less than the critical value (ARicrit) is computed at each time step, namely,
ARicrit =
Z 2000m
−1500m
Z η(x,t)
H(x) δRi(x, z)dz dx , (7) where
δRi(x, z) =
1 , ifRi(x, z) < 0.25 0 , otherwise,
and
Ri = − g ρ0
∂ρ
∂z
∂u
∂z
2
+∂v∂z2 . (8)
The area (7) will be time dependent. The time averaged value of ARicrit
ARicrit = 1 2T
Z 2T
0 ARicritdt (9)
where T is a M2 tidal period, is used to indicate the potential for mixing in different calculations.
The time and space integrated values of the energy losses due to bottom friction Eτ are computed according to
Eτ =
Z 2T 0
Z 15000
−5000mu~b(x, t)·τ~b(x, t)dx dt ,
where u~b(x, t) = (ub(x, t), vb(x, t) is the velocity vector in the nearest bottom grid cell, and ub(x, t) and vb(x, t) are the velocity components in the x-direction and the cross channel direction respectively. The bottom stress vector τ~b(x, t) is specified by
~
τb(x, t) =ρ0CD
qu2b +vb2u~b(x, t) (10)
where ρ0 is the reference density. The drag coefficient CD is given by CD = max[0.0025, κ2
(ln(zb/z0))2] (11) and zb is the distance of the nearest grid point to the bottom. The von Karman constantκis 0.4 and the bottom roughness parameter is chosen to bez0 = 0.01m, see Blumberg and Mellor (1987).
Similarly integrated values of the losses due to horizontal viscosityEH and vertical viscosity EV are computed
EH =
Z 2T 0
Z 15000
−5000m
Z η(x,t) H(x) ρAh
∂u
∂x
!2
+ ∂v
∂x
!2
+ ∂w
∂x
!2
dz dx dt ,
and
EV =
Z 2T 0
Z 15000
−5000m
Z η(x,t) H(x) ρAV
∂u
∂z
!2
+ ∂v
∂z
!2
+ ∂w
∂z
!2
dz dx dt .
Errors in the internal pressure gradient estimation may create artificial flows in σ-coordinate models, see for instance Haney (1991); Mellor et al. (1998);
Shchepetkin and McWilliams (2003); Thiem and Berntsen (2006). However, in the present experiments with relatively small grid sizes this is not a serious prob- lem. For instance in experiments with ∆x = 12.5 where the flow is driven only by the artificial internal pressure gradients, the kinetic energy Ekin remains less than 1×10−4J m−1, and Umax less than 1×10−4m s−1.
3 Results from the sub-critical experiments
In the first set of calculations the Froude number, F r = N Hu , at the sill crest, whereu is the across sill velocity andH is the sill depth, is strictly less than one throughout the tidal cycle. The half width of the sill is 1000 m. In general the results for this case are independent of the horizontal grid size for the range of grids examined here. This may for instance be seen from the temperature fields at maximum inflow in the second tidal cycle, see Figs. 2 and 3.
In Fig. 4 the values of the various diagnostics are given. The results for this low Froude number case are in general less sensitive to the grid size than to the choice of viscosity and diffusivity for the range of grid sizes applied. When using the present values of constant viscosity and diffusivity, the results are fairly robust to the grid size. However, from Fig. 4a and 4b we notice that the solution becomes less energetic as the grid size is reduced. This is due to slightly larger viscous losses and losses due to bottom friction when ∆x is reduced. In a mathematical sense we can produce a converged solution by reducing ∆x for this low Froude number case keeping the diffusivity constant and the viscosity large and constant.
There is a different message from the large eddy simulations. For this case, there is a clear increase in energy as the grid size is reduced. This is due to the grid dependent viscosities (1) allowing more small scale processes to enter the solution as ∆x is reduced. This may for instance be seen in the time series of vertical velocities for the case with ∆x= 12.5 m, see Fig. 5. It is also seen in the increase in the estimates of the areas of the water masses with Ri less than the critical value, see Fig. 4e. These low Ri water masses are found near the sea bed over the sill, but also higher up in the water column near the sill, see Fig. 6. For this case with grid dependent viscosity, the turbulent part of the water column will probably grow in extent if the grid size was reduced even further. In essence the non-linear nature of the problem (even at a sub-critical Froude number) in this region is such that there is a rapid cascade of energy from long waves to short waves and associated mixing. This cascade is balanced by the viscous terms, which in the case of constant values may lead to a rapidly converged solution, as ∆x is reduced (see in particular Fig. 4e). However, in the LES calculations, reducing
∆x decreases the viscous term and hence the balance in the energy cascade to short waves is affected. In a mathematical sense the problem changes with ∆x, and this is reflected in the rate of convergence.
Fig. 2 Fig. 3
Fig. 4 Fig. 5 Fig. 6
4 Results from the super-critical experiments
In the second set of calculations the Froude number is larger than one at maximum inflow and outflow, and the sill width is reduced to 500 m. This has the effect of substantially increasing the non-linear nature of the problem compared to previous set of experiments. The calculations show that for this case the numerical solutions are much more sensitive to the grid size, see Figs. 7 and 8. Note in particular the increase in amplitude of the internal waves at depth as ∆xdecreases from 25 m to 12.5 m in the constant viscosity case (Fig. 7). However, in the LES simulations these waves are smoothed by the viscosity and diffusivity terms (Fig.
8). In essence the solutions qualitatively change nature when the sub-grid closure is changed. In the constant viscosity and diffusivity experiments, internal waves appear in the lee of the sill as reported in XD06a,b. The period of these waves depend on the grid size, see the time series of vertical velocities in the lee of the sill given in Fig. 9a. For the highest resolution, ∆x = 12.5 m, the period is approximately 10-12 min consistent with the results produced with the MITgcm and ∆x= 10 m near the sill, see XD06b. The vertical velocities at (x, z) = (500 m,- 20 m) are generally much smaller in in the LES calculations, see Fig. 9b, and the time intervals during inflow with vertical excursions are much shorter. The vertical motions are also less harmonic, see Fig. 9b. However, in the time series for ∆x = 12.5 m a few oscillations with periods of approximately 12 minutes are still seen, even if they are much weaker than in the corresponding results produced with constant viscosity and diffusivity.
Since the internal waves produced in the LES calculations generally are much weaker, the solutions away from the sill, see Figs. 10 and 11, are also affected by the choice of sub-grid closure. With constant viscosity and diffusivity the strong internal waves produced in the lee propagate away from the sill. In the corre- sponding results from the large eddy simulations this signal is much weaker. The internal waves leaving the sill may subsequently break at topographic features away from the sill. In the constant viscosity/diffusivity calculations there may thus be a substantial turbulent dissipation away from the sill. In the LES calcu- lations more mixing will take place near the sill. The choice of sub-grid closure scheme thus affects the spatial distribution of the mixing, and the background stratification and the residual flow will be affected in long time integrations.
The propagation speed of the internal waves leaving the sill is approximately 0.3 ms−1 in the constant viscosity case. This internal wave speed has been repro- duced also with the MITgcm with an equidistant horizontal grid size of 12.5 m.
With a coarser resolution away from the sill, the wave propagation in the coarse grid part of the domain will not be captured, and this may also affect the spatial
distribution of mixing in the calculations.
By studying the diagnostics for the super-critical case given in Fig. 12 we note that they vary both with grid size and the sub-grid closure. However, contrary to the findings for the weakly forced case, the results from the LES calculations seem to be less sensitive to the grid size, see for instance the maximum velocities reported in Fig. 12d. The reason for this, can in part be understood in terms of the much larger short wave components in the constant viscosity/diffusivity solutions compared to the LES solutions. As in any numerical calculation, if the solution is dominated by short waves, a finer grid must be used than when the solution is composed of mainly long waves. Thus, for the constant viscosity case, the energy losses due to horizontal viscosity increase very much when the grid size is reduced, see Fig. 12g. These losses are due to the term ∂x∂ Ah∂u
∂x in the governing equations, and since Ah is constant, it must mean that ∂u∂x changes substantially with the grid size since the solution contains large short wave components. For the LES calculations it is noticeable that the corresponding losses are approximately equal for ∆x = 25 m and ∆x = 12.5 m since the short wave components of the solution are small.
Focusing on the areas of the water masses with Ri less than the critical value, we notice by comparing Figs. 4e and 12e that these areas are more than 10 times larger in this more strongly forced case. Furthermore, in the higher Froude number case the areas seem to increase as the grid size is reduced. However, the areas for the LES calculations using ∆x = 25 m and ∆x = 12.5 m are very similar, and almost equal to the value for ∆x = 12.5 m using constant viscosity and diffusivity. Even if the total areas are almost the same, the spatial distributions of these regions with low Ri are different, see Figs. 13 and 14. Comparing Figs.
10 and 13 we notice a good correspondence between the temperature fields and the Ri number distributions.
Fig. 7 Fig. 8 Fig. 9 Fig. 10 Fig. 11 Fig. 12 Fig. 13 Fig. 14
5 Sensitivity to the sub-grid closures
In the previous sections the results were found to be very sensitive to the sub-grid closure. In order to investigate this sensitivity further, four additional experiments for the super-critical flow case with ∆x= 12.5 m were performed. The horizontal diffusivity Kh is equal to 10−7m2s−1 in all studies. The values, or methods for computing, Ah, AV, and KV in all six experiments with ∆x = 12.5 m are given below
Experiment 1: Ah = 10−1m2s−1, andAV = 10−3m2s−1, and KV = 10−7m2s−1 as in XD06a,b.
Experiment 2:Ahis computed using (1) withCM = 0.2,AV andKV are estimated using using the Mellor-Yamada scheme as described in Section 2.
Experiment 3: Ah is computed using (1) with CM = 0.2, AV = 10−3m2s−1, and KV = 10−7m2s−1.
Experiment 4:Ah is computed using (1) with CM = 0.05,AV = 10−3m2s−1, and KV = 10−7m2s−1.
Experiment 5: Ah is computed using (1) with CM = 0.05, and AV and KV are estimated using using the Mellor-Yamada scheme.
Experiment 6:Ah = 10−1m2s−1, and AV andKV are computed with the method given in Pacanowski and Philander (1981) (hereafter PP81). The equations are
AV = A0
(1 +αRi)k +Ab , KV = AV
(1 +αRi)+Kb (12)
where we have chosen A0 = 0.015 m2s−1, α = 5, k = 2, Ab = 10−4m2s−1, Kb = 10−5m2s−1, and Ri is computed from (8).
The temperature fields after 5/4 T for experiments 3 to 6 are given in Fig. 15, and the diagnostics are given in Fig. 16, including also the diagnostics from the previous studies with ∆x = 12.5 m. By comparing Fig. 15 with Figs. 7d and 8d, we notice that the internal waves in the lee of the sill appear clearly in all the results produced using both AV and KV constant. If the vertical viscosity and diffusivity are computed with the Mellor-Yamada scheme or (12), the internal waves are far less evident.
The results are also sensitive to the horizontal viscosityAh, see Fig. 16. However, the solutions do not change qualitatively when replacing Ah = 10−1m2s−1 with the values computed using (1), or by adjusting the coefficient CM. A noticeable
result from Fig. 16 is that the energy losses due to horizontal viscosity, EH, are substantially reduced when replacing the constant values of AV and KV with values given by the Mellor-Yamada scheme or PP81. The explanation is that when using the Mellor-Yamada scheme or PP81 (equation 12), more of the tidal energy goes into irreversible mixing, and the internal wave generation is reduced.
The horizontal gradients ∂u∂x consequently become smaller, and the viscosity terms
∂
∂xAh∂u∂x are accordingly reduced.
From Fig. 16e we notice that the areas with less than critical Richardson numbers are fairly consistent for the first five experiment. The area is considerably reduced when using PP81. Generally larger values of AV and KV occur when using PP81 and more mixing may be the explanation. This is also consistent with larger energy losses due to bottom friction and vertical viscosity in Experiment 6, see Figs. 16f and 16h.
Fig. 15 Fig. 16
6 Discussion
A cross sectionalσ-coordinate model has been applied to study the dynamics near a sill in an idealised loch. The buoyancy frequency is assumed to be constant in all studies. For studies with sub-critical flow over a wide sill, the results are very robust to changes in the horizontal grid size in the range from 100 m to 12.5 m.
With stronger tidal forcing, a narrower sill, and super-critical flow over the sill during inflow, we get flow separation and generation of high amplitude internal waves in the lee of the sill. The frequency and the amplitude of these waves changes noticeably as the horizontal grid size is reduced from 100 m to 12.5 m.
With a grid size of 12.5 m and constant values of viscosity and diffusivity, the results are very similar to those presented XD06a,b.
Calculations show that the results presented are sensitive to the choice of sub- grid closure. In the low Froude number studies, the results are more sensitive to the sub-grid closure scheme than to changes in the grid size, even if the re- sults are generally consistent for the sub-critical case. In this case the short wave components are less significant than in the high Froude number calculations. For the high Froude number case, the solution may change qualitatively nature when changing the sub-grid closure scheme. With constant and low values of vertical diffusivity as in XD06a, strong internal waves are generated on inflow in the lee of the sill. These waves subsequently propagate away from the sill with almost constant wave speed. In this case processes at the scale of the grid resolution and sub-grid scale processes are more important than in the low Froude number cal- culations. By using a Richardson number dependant vertical diffusivity (Mellor and Yamada (1982) or Pacanowski and Philander (1981)) more of the tidal energy goes into vertical irreversible mixing on inflow, and far less into the generation of internal waves propagating away from the sill.
There is a large literature on sub-grid scale closures in ocean modelling, see Haid- vogel and Beckmann (1999); Kantha and Clayson (2000); Burchard (2002). One could, however, hope that as the grid sizes are gradually reduced, there will be less unresolved important processes, and that the sensitivity to the sub-grid closure scheme therefore should be reduced. The results from the present studies indicate that with a grid size of 12.5 m, there are still important unresolved processes, par- ticularly in studies with super-critical flow. This may not be surprising, consider- ing the results in Farmer and Armi (1999); Cummins et al. (2003). However, many studies of internal wave generation, propagation, breaking, and mixing have been performed with grid sizes in the range from 1 m to 100 m. Especially in the gener- ation and breaking phases there may be important unresolved processes in such studies, see for instance Klymak and Gregg (2001); Vlasenko et al. (2002); Legg
and Adcroft (2003); Bourgault and Kelley (2003); Wang (2006); Lamb (2007).
One may consider to apply numerical fine scale simulations, in essence direct nu- merical simulations, to get the ’true’ solution and to estimate diffusivities and viscosities, and finally to assess or improve sub-grid scale parameterisations, see Smyth et al. (2005); Berntsen et al. (2006). However, in direct numerical simula- tions the grid sizes applied are typically smaller than 0.01 m allowing only studies for tank scale problems. In physical oceanographic problems such as those con- sidered here, inaccuracies in, for example, topographic variations and changes in bed types and forms on scales of the order of 10 m, which are rarely known, will have an appreciable effect and may negate any advantages of direct numerical simulations even if it was computationally possible. Several complementary ap- proaches are therefore necessary to improve our understanding of the processes occurring in the lee of the sill during super-critical inflow. One approach may still be to use numerical models and to focus more on the area near the sill, allowing a horizontal grid size closer to 1 m as in Cummins et al. (2003). As this grid size is significantly larger than that used in direct numerical simulations, there will still be unresolved processes. However, shear interface instabilities, and the effects of these instabilities, may still be captured. Another approach may be to seek high resolution measurements. A potential problem with this approach is that it may be difficult to find a loch or a fjord of sufficient limited extent that a synoptic data set of density measurements can be made particularly in the sill region for model investigations, and subsequent measurements for model validation. In the Arctic there may, however, be fjords with almost constant buoyancy frequency at least during parts of the year, see Fer and Widell (2007) which would at least provide initial conditions, provide detailed topography was available. A third approach would be to undertake laboratory experiments of the kind reported in Sveen et al.
(2002); Guo and Davies (2003); Guo et al. (2005), but for conditions similar to those used in the present study. There are strengths and weaknesses with all three approaches. However, a combined use of all three may lead to deeper insight into the processes in the lee of the sill under conditions of strong inflow, and not least into the effects of these processes on the larger scale fields.
Three dimensional effects are ignored in the present study. In a real loch or fjord, vortices will be created behind headlands and they will influence the wave generation and mixing, see Klymak and Gregg (2001); McCabe et al. (2006);
Zhao et al. (2006). However, with limiting computer resources, three dimensional studies will require the use of larger grid sizes. The choice made here is therefore to apply a cross sectional model with as high a resolution as possible to allow a better representation of the waves in the lee of the sill. On a much finer scale there will also be three dimensional effects associated with shear instabilities, see D¨ornbrack (1998); Smyth et al. (2005). However, these three dimensional effects
may not be resolved with the grid sizes used here.
Acknowledgements. This work was done during the first author’s sabbatical year at Proudman Oceanographic Laboratory. He thanks the people at the lab- oratory for their great hospitality, and the Faculty of Mathematics and Natural Sciences, University of Bergen, for support during the sabbatical year.
Figure captions
Fig. 1. Initial temperature and topography for half sill width equal to 1000 m.
Fig. 2. The temperature fields for the sub-critical case after t = 5/4 T (where T is the M2 tidal period) for the experiments with constant viscosity and diffusivity.
The results are for ∆x = 100 m, 50 m, 25 m, and 12.5 m from top to bottom.
Fig. 3. The temperature fields for the sub-critical case after t = 5/4 T for the large eddy simulations. The results are for ∆x = 100 m, 50 m, 25 m, and 12.5 m from top to bottom.
Fig. 4. Sensitivity of a)Ekin(5/4T), b)Ekin(2T), c)Ekin−bclin(2T), d)Umax(5/4T), e)ARicrit, f)Eτ, g)EH, and h)EV to the grid resolution for the sub-critical case.
The results for the constant viscosity and diffusivity experiments are given with solid lines and squares. The results for the large eddy simulations are given with dashed lines and circles.
Fig. 5. Vertical velocities for the sub-critical case atx = 1000 m during a part of the second tidal cycle for the large eddy simulations with ∆x = 12.5 m. (Time is scaled with T.)
Fig. 6. Richardson number distribution for the sub-critical case att= 13/8 T for the large eddy simulations with ∆x = 12.5 m.
Fig. 7. The temperature fields for the super-critical case after t = 5/4 T for the experiments with constant viscosity and diffusivity. The results are for ∆x = 100 m, 50 m, 25 m, and 12.5 m from top to bottom.
Fig. 8. The temperature fields for the super-critical case after t = 5/4 T for the large eddy simulations. The results are for ∆x = 100 m, 50 m, 25 m, and 12.5 m from top to bottom.
Fig. 9. Time series of vertical velocities in m s−1 at (x, z) = (500 m, -20 m) for the super-critical case. The results for ∆x= 100 m are given with solid black lines, the results for ∆x= 50 m with red lines, the results for ∆x= 25 m with blue lines, and the results for ∆x= 12.5 m with dashed black lines. The time is given in minutes.
The results from the constant viscosity and diffusivity experiments are given in the top figure, and the corresponding results from the large eddy simulations are given below. Notice the change in scale on the axis for the vertical velocities.
Fig. 10. The temperature fields for the super-critical case after t = 3/2 T for a) the experiment with constant viscosity and diffusivity and b) the large eddy simulation. The results are for ∆x = 12.5 m.
Fig. 11. The temperature fields after t = 2 T for a) the experiment with constant viscosity and diffusivity and b) the large eddy simulation. The results are for ∆x
= 12.5 m.
Fig. 12. Sensitivity of a)Ekin(5/4T), b)Ekin(2T), c)Ekin−bclin(2T), d)Umax(5/4T), e) ARicrit, f) Eτ, g) EH, and h) EV to the grid resolution for the super-critical case. The results for the constant viscosity and diffusivity experiments are given with solid lines and squares. The results for the large eddy simulations are given with dashed lines and circles.
Fig. 13. The Richardson number distribution for the super-critical case after t
= 3/2 T for a) the experiment with constant viscosity and diffusivity and b) the large eddy simulation. The results are for ∆x = 12.5 m.
Fig. 14. The Richardson number distribution for the super-critical case after t = 2 T for a) the experiment with constant viscosity and diffusivity and b) the large eddy simulation. The results are for ∆x = 12.5 m.
Fig. 15. The temperature fields for the super-critical case after t = 5/4 T for a) Experiment 3, b) Experiment 4, c) Experiment 5, and d) Experiment 6.
Fig. 16. Sensitivity of a)Ekin(5/4T), b)Ekin(2T), c)Ekin−bclin(2T), d)Umax(5/4T), e)ARicrit, f) Eτ, g) EH, and h) EV for the super-critical case to the sub-grid clo- sure scheme for the experiments 1 to 6.
−6000 −4000 −2000 0 2000 4000 6000 8000 10000 12000 14000
−100
−90
−80
−70
−60
−50
−40
−30
−20
−10 0
distance [m]
depth [m]
8.25 8.25
8.5 8.5
8.75 8.75
9 9
9.25 9.25
9.5 9.5
9.75 9.75
10 10
10.25 10.25
10.5 10.5
10.75 10.75 10.75
11 11 11
11.25 11.25 11.25
11.5 11.5 11.5
11.75 11.75 11.75
12 12 12
12.25 12.25 12.25
12.5 12.5 12.5
12.75 12.75 12.75
13 13 13
T (ci=0.25 C)
(a)
Fig. 1. Initial temperature and topography for half sill width equal to 1000 m.
11 12 12.5
T (ci=0.25 C)
distance [m]
depth [m]
0 200 400 600 800 1000 1200 1400 1600 1800 2000
−100
−90
−80
−70
−60
−50
−40
−30
−20
−10 0
(a) ∆x = 100 m
10 11.5 12.5
T (ci=0.25 C)
distance [m]
depth [m]
0 200 400 600 800 1000 1200 1400 1600 1800 2000
−100
−90
−80
−70
−60
−50
−40
−30
−20
−10 0
(b) ∆x = 50 m
11.5 13
T (ci=0.25 C)
distance [m]
depth [m]
0 200 400 600 800 1000 1200 1400 1600 1800 2000
−100
−90
−80
−70
−60
−50
−40
−30
−20
−10 0
(c) ∆x = 25 m
11 12
12.5
T (ci=0.25 C)
distance [m]
depth [m]
0 200 400 600 800 1000 1200 1400 1600 1800 2000
−100
−90
−80
−70
−60
−50
−40
−30
−20
−10 0
(d) ∆x= 12.5 m
Fig. 2. The temperature fields for the sub-critical case after t = 5/4 T (where T is the 23
10 11 11.5 12
12.5
T (ci=0.25 C)
distance [m]
depth [m]
0 200 400 600 800 1000 1200 1400 1600 1800 2000
−100
−90
−80
−70
−60
−50
−40
−30
−20
−10 0
(a) ∆x = 100 m
11
T (ci=0.25 C)
distance [m]
depth [m]
0 200 400 600 800 1000 1200 1400 1600 1800 2000
−100
−90
−80
−70
−60
−50
−40
−30
−20
−10 0
(b) ∆x = 50 m
T (ci=0.25 C)
distance [m]
depth [m]
0 200 400 600 800 1000 1200 1400 1600 1800 2000
−100
−90
−80
−70
−60
−50
−40
−30
−20
−10 0
(c) ∆x = 25 m
9.5 11.5
12 12.5
T (ci=0.25 C)
distance [m]
depth [m]
0 200 400 600 800 1000 1200 1400 1600 1800 2000
−100
−90
−80
−70
−60
−50
−40
−30
−20
−10 0
(d) ∆x= 12.5 m 24
10 20 30 40 50 60 70 80 90 100 110 2.7
2.75 2.8 2.85 2.9 2.95 3 3.05 3.1 3.15
3.2x 105
DX
Kinetic energy [J/m]
(a) Ekin(5/4T)
10 20 30 40 50 60 70 80 90 100 110
0.8 0.85 0.9 0.95 1 1.05 1.1 1.15 1.2 1.25
1.3x 105
DX
Kinetic energy [J/m]
(b) Ekin(2T)
10 20 30 40 50 60 70 80 90 100 110
7 7.5 8 8.5 9 9.5 10 10.5 11 11.5
12x 104
DX
Kinetic energy [J/m]
(c) Ekin−bclin(2T)
10 20 30 40 50 60 70 80 90 100 110
0.06 0.065 0.07 0.075 0.08 0.085 0.09
DX
U [m/s]
(d) Umax(5/4T)
10 20 30 40 50 60 70 80 90 100 110
500 1000 1500 2000 2500 3000
DX
Area [m**2]
(e) ARicrit
10 20 30 40 50 60 70 80 90 100 110
1.5 2 2.5 3 3.5 4 4.5 5 5.5
6x 104
DX
Energy loss [J/m]
(f) Eτ
10 20 30 40 50 60 70 80 90 100 110
1000 2000 3000 4000 5000 6000 7000
DX
Energy loss [J/m]
(g)EH
10 20 30 40 50 60 70 80 90 100 110
1 1.5 2 2.5 3 3.5 4 4.5
5x 104
DX
Energy loss [J/m]
(h) EV
Fig. 4. Sensitivity of a)Ekin(5/4T), b)Ekin(2T), c)Ekin−bclin(2T), d)Umax(5/4T), e) ARicrit, f) Eτ, g) EH, and h) EV to the grid resolution for the sub-critical case. The results for the constant viscosity and diffusivity experiments are given with solid lines
25
1.1 1.15 1.2 1.25 1.3 1.35 1.4 1.45 1.5 1.55 1.6
−70
−60
−50
−40
−30
−20
−10 0
Time
Depth [m]
0
0
0
0 0
0 0
0
0
0
0 0
0
0 0
0 0
0 0
0
0 0
0
00
0.2 0.2
0.2
0.2 0.2
0.2
0.2 0.2
0.2
0.2 0.2
0.2 0.2
0.2 0.4
0.4
0.4 0.4
0.6
W (ci=0.2 mm/s)
−0
−0
−0
−0
−0
−0
−0
−0
−0 −0
−0
−0
−0
−0
−0
−0
−0
−0
−0
−0
−0
−0
−0
−0
−0
−0.2 −0
−0.2
−0.2
−0.2
−0.2
−0.2
(a)
Fig. 5. Vertical velocities for the sub-critical case at x = 1000 m during a part of the second tidal cycle for the large eddy simulations with ∆x = 12.5 m. (Time is scaled with T.)
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
−1000 −800 −600 −400 −200 0 200 400 600 800 1000
−30
−25
−20
−15
−10
−5 0
Ri (ci=0.05 )
distance [m]
depth [m]
(a)
Fig. 6. Richardson number distribution for the sub-critical case at t = 13/8 T for the large eddy simulations with ∆x= 12.5 m.
9
10 9.5 10
10.5
11 11
11.5
11.5
11.5 12
12
12
T (ci=0.25 C)
distance [m]
depth [m]
0 200 400 600 800 1000 1200 1400 1600 1800 2000
−100
−90
−80
−70
−60
−50
−40
−30
−20
−10 0
(a) ∆x = 100 m
9
9.5 10
10.5
10.5 11
11
11.5 11.5
12
12
12
12.5
T (ci=0.25 C)
distance [m]
depth [m]
0 200 400 600 800 1000 1200 1400 1600 1800 2000
−100
−90
−80
−70
−60
−50
−40
−30
−20
−10 0
(b) ∆x = 50 m
8.75 9.25 9 9.759.5
10
10.25 10.510.75
11 11.25 11
11.25 11.5
11.5 11.75
12
12 12
12.25 12.5
T (ci=0.25 C)
distance [m]
depth [m]
0 200 400 600 800 1000 1200 1400 1600 1800 2000
−100
−90
−80
−70
−60
−50
−40
−30
−20
−10 0
(c) ∆x = 25 m
9 8.5
9
9.5 9.5
9.5
10
10
10
10
10.511 1110.5 1110.5
11 11.5
11.5
11.5 11.5
11.5 12 12 12 12
1212 12
12
12 12 12
12
12 12.5 13
T (ci=0.25 C)
distance [m]
depth [m]
0 200 400 600 800 1000 1200 1400 1600 1800 2000
−100
−90
−80
−70
−60
−50
−40
−30
−20
−10 0
(d) ∆x= 12.5 m
Fig. 7. The temperature fields for the super-critical case after t = 5/4 T for the exper- 27
8.5 9
10 10.5 11
11.5
12 12.5
T (ci=0.25 C)
distance [m]
depth [m]
0 200 400 600 800 1000 1200 1400 1600 1800 2000
−100
−90
−80
−70
−60
−50
−40
−30
−20
−10 0
(a) ∆x = 100 m
9
9.5 10
10.5 11
11.5 12
12.5 13
T (ci=0.25 C)
distance [m]
depth [m]
0 200 400 600 800 1000 1200 1400 1600 1800 2000
−100
−90
−80
−70
−60
−50
−40
−30
−20
−10 0
(b) ∆x = 50 m
8.5 8.75 9
9.25 9.5
10 11.259.75 10.2510.7510.5 11 11.5 11.75
11.75
12.25 12.5
12.75
13
T (ci=0.25 C)
distance [m]
depth [m]
0 200 400 600 800 1000 1200 1400 1600 1800 2000
−100
−90
−80
−70
−60
−50
−40
−30
−20
−10 0
(c) ∆x = 25 m
8.5 9.59 10
10.5 11 11.5
12
12 12.5
12.5
13
13
T (ci=0.25 C)
distance [m]
depth [m]
0 200 400 600 800 1000 1200 1400 1600 1800 2000
−100
−90
−80
−70
−60
−50
−40
−30
−20
−10 0
(d) ∆x= 12.5 m 28
850 900 950 1000 1050 1100
−0.2
−0.15
−0.1
−0.05 0 0.05 0.1 0.15 0.2
t
W [m/s]
(a)
850 900 950 1000 1050 1100
−0.05
−0.04
−0.03
−0.02
−0.01 0 0.01 0.02 0.03 0.04 0.05
t
W [m/s]
(b)
Fig. 9. Time series of vertical velocities in m s−1 at (x, z) = (500 m, -20 m) for the super-critical case. The results for ∆x = 100 m are given with solid black lines, the results for ∆x = 50 m with red lines, the results for ∆x = 25 m with blue lines, and the results for ∆x = 12.5 m with dashed black lines. The time is given in minutes.
The results from the constant viscosity and diffusivity experiments are given in the top figure, and the corresponding results from the large eddy simulations are given below.
Notice the change in scale on the axis for the vertical velocities.