Internal pressure gradient errors in σ-coordinate ocean models in high resolution fjord studies

Ocean Modelling 92 (2015) 42–55

Contents lists available at ScienceDirect

Ocean Modelling journal homepage: www.elsevier.com/locate/ocemod

Internal pressure gradient errors in σ -coordinate ocean models in high resolution fjord studies Jarle Berntsen a,∗, Øyvind Thiem b, Helge Avlesen b a b

Department of Mathematics, University of Bergen, Allegt. 41, N-5008 Bergen, Norway UNI Computing, Uni Research, Thormøhlensgt. 55, N-5008 Bergen, Norway

a r t i c l e

i n f o

Article history: Received 30 January 2015 Revised 19 May 2015 Accepted 27 May 2015 Available online 9 June 2015 Keywords: Sigma-coordinates Internal pressure gradients High resolution Steepness fjords

a b s t r a c t Terrain following ocean models are today applied in coastal areas and fjords where the topography may be very steep. Recent advances in high performance computing facilitate model studies with very high spatial resolution. In general, numerical discretization errors tend to zero with the grid size. However, in fjords and near the coast the slopes may be very steep, and the internal pressure gradient errors associated with σ models may be significant even in high resolution studies. The internal pressure gradient errors are due to errors when estimating the density gradients in σ -models, and these errors are investigated for two idealized test cases and for the Hardanger fjord in Norway. The methods considered are the standard second order method and a recently proposed method that is balanced such that the density gradients are zero for the case ρ = ρ (z ) where ρ is the density and z is the vertical coordinate. The results show that by using the balanced method, the errors may be reduced considerably also for slope parameters larger than the maximum suggested value of 0.2. For the Hardanger fjord case initialized with ρ = ρ (z ), the errors in the results produced with the balanced method are orders of magnitude smaller than the corresponding errors in the results produced with the second order method. © 2015 Elsevier Ltd. All rights reserved.

1. Introduction Sigma coordinate ocean models are applied in numerous marine science studies. Their ability to represent the processes near the free surface and near the bottom makes these models attractive even if there are errors associated with the estimation of the internal pressure gradients (IPG). This is discussed in textbooks, see Haidvogel and Beckmann (1999), Kantha and Clayson (2000), Griffies (2004). The IPG errors are addressed in many papers including Haney (1991), Beckmann and Haidvogel (1993), Mellor et al. (1994), McCalpin (1994), Chu and Fan (1997), Mellor et al. (1998), Ezer et al. (2002), Shchepetkin and McWilliams (2003), Martinho and Batteen (2006), Sikiric´ et al. (2009), Berntsen and Oey (2010). For a more complete reference list see Berntsen (2011). Let x and z be the horizontal and vertical coordinates in an orthogonal Cartesian system, ρ the density and H the depth. Assuming that the free surface is at z = 0 m, the sigma transformation is given by σ ≡ z/H. The x-components of the density gradients in the σ -

∗

Corresponding author. Tel.: +47 55 584854; fax: +47 55 589672. E-mail addresses: [email protected], [email protected] (J. Berntsen), [email protected] (Ø. Thiem), [email protected] (H. Avlesen). http://dx.doi.org/10.1016/j.ocemod.2015.05.009 1463-5003/© 2015 Elsevier Ltd. All rights reserved.

coordinate system are given by

 ∂ρ  ∂ρ σ ∂ H ∂ρ − = . ∂ x z ∂ x H ∂ x ∂σ

(1)

After discretization using second order differences and averaging, numerical errors are introduced and for the case ρ = ρ (z ) the den-



sity gradients may be estimated to be non-zero even if ∂ρ ∂ x  is zero z analytically. The two terms on the right hand side of statement (1) should for this case cancel, but when using standard discretization techniques they may not. A sketch of the variables involved when approximating the density gradients on a C-grid is given in Fig. 1. We let x = (xi−1 + xi )/2, σ = (σk−1 + σk )/2, H = (Hi−1 + Hi )/2, ρ xi,k = (ρi−1,k + ρi,k )/2, ρ σi,k = (ρi,k−1 + ρi,k )/2, δ x = xi − xi−1 , δσ = σk−1 − σk , and let δ x and δ σ be the discrete central difference operators in x- and σ -direction, respectively. The density gradient approximated at (x, σ ) = (x, σ ) with the second order POM-method may then be given by

 ∂ρ  δx ρ σ σ δx H δσ ρ x − ( x, σ ) ∼ . ∂ x z δx H δ x δσ

(2)

A modified perfectly balanced method is suggested in Berntsen (2011). The idea is to calibrate or adjust the approximation (2) such that the computed value becomes zero in all points where it is needed

J. Berntsen et al. / Ocean Modelling 92 (2015) 42–55

Fig. 1. The locations of the variables involved when approximating (1). The positions of the u-points are marked on the lateral edges of the grid cells. The density gradient is approximated at (x, σ ) = (x, σ ).

for the case ρ = ρ (z ). This is done by replacing σ in (2) with σ SOPT such that the two terms on the right hand side of (2) cancel. In Berntsen (2011) numerical evidence is given that σ SOPT tend to σ quadratically as δ x and δσ tend to zero. The testing in Berntsen (2011) is performed with regional and large scale applications in mind and include tests for the seamount case, see Beckmann and Haidvogel (1993), Mellor et al. (1998), and the Nordic Seas. In this paper, the focus is on higher resolution smaller scale studies over domains with steep topography such as fjords. Following Mellor et al. (1994) the error of the approximation (2) is to leading order

H E= 4

  2  2 ∂ ρ δx H 2 δx H 2 (δσ ) − σ . δx ∂ z2 H

(3)

Repeating the error expansion around (x, σ ) = (x, σ ) for the SOPTmethod, the corresponding error becomes

ESOPT =

 ∂ρ δx H  σ − σ SOPT δx ∂z   2  2 ∂ ρ H δx H H δ x (δσ )2 − σ σ SOPT + 4 δx H ∂ z2 +

δx H  2

σ − σ SOPT

 ∂ 2ρ . ∂ x∂ z

(4)

In (4) only the three lowest order terms are included. Terms one and three on the right hand side of (4) are new compared to statement (3). However, they both tend to zero as σ SOPT tends to σ . Term two is a modified version of the right hand side of (3) with σ 2 replaced by σ σ SOPT . It may be noted that previously another method for balancing the IPG errors to zero for ρ = ρ (z ) was presented by Stelling and Van Kester (1994). Their method was based on interpolation to zlevels before estimating the gradients. For each velocity point, two estimates of the density gradient, bounding the exact density gradient, are computed, and the gradient that is smallest in magnitude is used as the approximation to the real density gradient. If the two estimates have different signs, the estimate of the density gradient is set to zero. Stelling and Van Kester (1994) method has inspired the present work. The method may, however, lead to systematic underestimation of density gradients and geostrophic transports, see Slørdal (1997). Over the years, different methods for estimating the internal density and internal pressure gradients in terrain following models have been tested both on idealized test cases and for more realistic model

43

domains. The Gaussian seamount case introduced in Beckmann and Haidvogel (1993) has become a standard case. The methods are also tested in a range of more realistic regional and large scale model studies, see for instance Mellor et al. (1994), Barnier et al. (1998), Berntsen (2002), Shchepetkin and McWilliams (2003), Berntsen and Thiem (2007), Adcroft et al. (2008), Ciappa (2008), Berntsen and Oey (2010), Berntsen (2011). In these studies the horizontal grid size is in the range from 1 to 200 km. The terrain following models are also applied in a range of fine scale studies for fjords, throughflows, and inlets. In three-dimensional studies the horizontal grid size is typically in the range from 50 m to 1 km (Ali et al., 2011; 2013; Asplin et al., 1999; Eliassen et al., 2001; Lynge et al., 2010; Myksvoll et al., 2014; Sutherland et al., 2011). To allow detailed studies of vertical motion, two-dimensional versions of terrain following models are often applied and in such studies the horizontal grid size may be close to 1 m (Berntsen et al., 2008; 2009; Cummins, 2000; Lamb, 2004). Fjords are formed by glaciers and tend to have a cross-fjord U-shape. The slopes may become very steep, and there may be underwater vertical cliffs where ∂∂Hx is infinity. Numerical errors should normally tend to zero with the grid size. However, with very steep slopes, the internal pressure errors may still be very large when using high resolution terrain following models, see Eq. (3). Haney (1991) stated that the hydrostatic consistency condition

   σ δx H   <1 r= δσ δ x 

(5)

should be satisfied. It is later shown that σ -models may produce accurate results also for values of r larger than 1 (Ezer et al. (2002); Mellor et al. (1994)). In Mellor et al. (1998) it is argued that this concept is not particularly meaningful. There is still some discussion on the usefullness of the hydrostatic concistency concept and Sikiric´ et al. (2009) for instance regard r to be a ’very useful measure’. In the present paper further evidence to support the view that r is not a useful parameter will be given. The slope parameter

   δx H    sp =  2H 

(6)

is, however, a more relevant concept. For a real slope α and a grid size δ x, δ x H will be αδ x showing that for a given slope sp will tend to zero with the grid size. This slope parameter is used in Mellor et al. (1994), without the factor 2 in the denominator, and later on in Mellor et al. (1998) and Martinho and Batteen (2006). Numerical evidence is given that the acceptable values are less than 0.2. In model studies for fjords we may, however, have one very shallow grid cell near land and the next cell in the deep fjord giving δ x H ∼ H and sp ∼ 0.5 indicating that it may be difficult to satisfy the sp less than 0.2 criterion in fjord studies with σ -models without substantial topography smoothing. An overview of the physical oceanography of fjords is given in Farmer and Freeland (1983). The fjords are formed by glacial carving, creating steep side walls. The vertical salinity, temperature, and density profiles are affected by fresh water supply at the surface, heat exhanges, circulation, and mixing, see for instance Pickard (1961), Hansen and Rattray Jr. (1966), Stigebrandt (1981). There will be a brackish layer near the surface creating a vertical profile of density that in many situations may be approximated with an exponential density profile. Mixing near the surface may create two-layer systems that may be approximated with tanh density profiles. The methods for estimating the density gradients are tested using two idealized two-column test cases based on both types of density profiles in combination with steep topography. The main findings from the two idealized cases are robust to the choice of vertical density profile, and in tests based on the topography from the Hardanger fjord an exponential density profile is used. The effects of topograhy smoothing are also considered in all test cases.

44

J. Berntsen et al. / Ocean Modelling 92 (2015) 42–55

ative values may accordingly cancel. The integral (12) may be computed exactly when inserting (11) and Tg is given by



ρβ g H0 Lz  β Lx /Lz e − e−β Lx /Lz ρ0 f β      α L2 β Lx β Lx − 2z eβ Lx /Lz − 1 + e−β Lx /Lz +1 Lz Lz β 2  L − z eβ Lx /Lz − e−β Lx /Lz β

 (β +α )L /L  L2z −H0 /Lz −(β +α )Lx /Lz x z + e e −e . (β + α )

Tg = −

Fig. 2. The configuration of test case 1 with a density isoline indicated in red.

2. Two-column test case with exponential density profile To investigate the σ -errors a simple two-dimensional test case is constructed. The bottom is tilted with an angle θ and a slope α . The depth H = H0 − α x. An exponential density profile tilted the other way with an angle φ and a tilt β is given by:

ρ (x, z ) = ρ0 − ρ e(z+β x)/Lz

(7)

where Lz is the vertical scale. The configuration is outlined in Fig. 2, and the test case will be referred to as test case 1. For this case the density gradients are given by

∂ρ

ρβ (z+β x)/Lz (x, z ) = − e . ∂x Lz

(8)

Following Slørdal (1997), this may be translated into force per unit mass

g ∂p = ρ0 ∂ x ρ0 1



0 z

∂ρ dz , ∂x

(9)

and after division by the Coriolis parameter f equivalent geostrophic velocities are obtained

1 f ρ0

vg (x, z ) =

∂p . ∂x

(10)

In all experiments, g = 9.81 m/s2 , f = 1.3 × 10−4 s−1 , and the reference density ρ0 = 1020.0 kg/m3 . For the density gradient (8), vg is given by

vg (x, z ) = −



ρβ g β x/Lz e − e(z+β x)/Lz . f ρ0

(11)

The exact solutions (8) and (11) will be used to check the errors in the numerical approximations. It may be noted that as expected for the profile (7), the exact solutions do not depend on the bottom slope α . However, the numerical approximations will depend on α . To reduce the σ -errors, it is suggested to smooth H. This means to reduce α and the hydrostatic consistency parameter r. By doing this one may certainly improve the quality of the approximations to vg (x, z) at a point with undisturbed depth. However, the changes in the topography will also affect the geostrophic transports. To illustrate this, the geostrophic velocities are integrated across a two-dimensional cross section



Tg =

Lx −Lx



0

−H (x )

vg (x, z )dzdx .

(12)

It may be noted that the geostrophical velocities vg may change sign over the two-dimensional cross section and positive and neg-

(13)

The solution (13) is exact for Lx > H0 /α and Tg is a function of α . Statement (13) will be applied to quantify the errors in corresponding numerical approximations, with and without topography smoothing. Statement (8) is used to compute ∂ρ ∂ x (0, z ) and the following numerical methods are used to produce corresponding numerical approximations (i) the standard POM-method, (ii) the optimal σ -method (SOPT), and (iii) the z-coordinate method using central differences. The POM-method (Blumberg and Mellor (1987)) and the optimal σ -method are based on the C-grid (Mesinger and Arakawa (1976)) and described in Berntsen (2011). The densities are given at lateral distances x/2 from velocity points given at x = 0 when computing estimates to the density gradients. The z-coordinate method considered is also based on the C-grid and central differences are taken along iso-surfaces of z to compute ∂ρ /∂ x. It may be noted that usually the use of z-level methods to compute density gradients in σ -models involves vertical interpolation, see Ciappa (2006); Fortunato and Baptista (1996); Mellor et al. (2002); Stelling and Van Kester (1994). However, in the present idealized tests with known density distributions ρ (x, z) there is no need for vertical interpolation before estimating ∂ρ ∂ x . The errors when using the z-coordinate method are strictly due to errors involved when approximating ∂ρ /∂ x. With a tilting topography and a velocity point at x = 0, there will be a step in the numerical grid at x = 0 when applying this method. For velocity points at the wall of the step, underneath the shallowest cell, ∂ρ /∂ x is set to zero, consistent with the no-flow boundary condition at a closed wall. In the first set of experiments, H0 = 250 m, Lx = 50 m, and x = 50 m which is an achievable grid size for high resolution studies. The density surface has a tilt β = 1/100. For all methods there are 200 equidistant layers. Furthermore, ρ = 1.0 kg/m3 , and the vertical scale Lz = 50 m. The exact solution and the approximations are computed for z ∈ [−250m, 0m] and for α ∈ [0, 9]. For α = 9, the depth at x = − x/2 is 475 m and 25 m at x = x/2. The errors for the POM- and SOPT-methods are based on numerical approximations to Eq. (1) and the exact value given in the right hand side of Eq. (8). The errors for the z-coordinate method is based on central difference approximations to the left hand side of (8) and the exact value. The exact values of the density gradients are small and the relative errors are given in the right column of the figure. The errors in Fig. 3 are given using log10 of the absolute values due to the wide range of magnitudes. The errors in the upper part of the water column are positive above the minima lines seen in Figs. 3a–d and negative below them. For the z-coordinate method the errors are positive above the lines of discontinuity in Figs. 3e and f. For cases when the z-coordinate method has density points on both sides of a velocity point, this method is as expected clearly the best with a constant relative error of 10−5.38 , see Fig. 3f, regardless of α . When the velocity points are on the step, the relative errors become 1, which explains the lines of discontinuity in Figs. 3e and f. For such situations the POM-method produces estimates of the density gradients with relative errors larger than 1, see the large area of positive values of log10 of the relative errors in Fig. 3b.

J. Berntsen et al. / Ocean Modelling 92 (2015) 42–55

45

Fig. 3. Errors in the approximations to ∂ρ /∂ x, see the right hand side of statement (8), at x = 0 as functions of z and α . In these experiments β = 0.01, x = 50 m, and z = 1.25 m. In the top panel log10 of the errors in the approximations from the POM-method ( ∂ρ ∂ x POM ) are given, in the middle panel of the errors in the approximations from the SOPT-method ∂ρ ( ∂ρ ∂ x SOPT ), and in the lower panel the errors in the approximations from the z-coordinate model ( ∂ x Z ) are given. Absolute errors are given to the left and relative errors to the right.

The SOPT-method is balanced such that the errors are zero for

β = 0. For β = 0.01, the errors are also considerably smaller than cor-

0

2 3 4 5 6 7 8 9 10

−100

1

−50

11 12 13 14 15 16 17 18 19 20 21 22 2 243 25

Z −200

−250 0

2 3 4 5 6 7 8 9 10

−150

1

responding errors for the POM-method, see Figs. 3c and d, and the relative errors are smaller than 1 even for the larger values of α . The values of the hydrostatic consistency parameter r for α ≤ 1 are given in Fig. 4. The value of r = 358.2 for α = 9 at the bottommost velocity point. For this case with a high vertical resolution, r becomes larger than the critical value 1 for most values of z and α considered. By doubling the number of equidistant layers to 400, the r values will double, see statement (5). However, it has been tested that the errors given in Fig. 3 will not be affected, supporting that hydrostatic consistency is not a useful concept, see the discussion in Berntsen and Oey (2010). The results from Figs. 3a to d show nevertheless that the internal pressure gradient errors are increasing with larger values of α for both σ -methods. Due to the steps in the topograhy this may also be the case for the z-coordinate method, see Figs. 3e and f. It will accordingly be beneficial to smooth H if the only goal is to reduce the internal pressure gradient errors. The slope parameter sp , see Eq. (6), is a more useful concept and the relationship between sp

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

1

Fig. 4. The values of the hydrostatic consistency parameter r for α ∈ [0, 1].

46

J. Berntsen et al. / Ocean Modelling 92 (2015) 42–55

and α is given by s p = α x/(2H0 ) and for the present values of x and H0 , s p = 0.1α . Acceptable values are suggested to be less than 0.2 (Martinho and Batteen (2006)). From Fig. 3b we notice that the relative errors when using the POM-method are, for most values of z, larger than 1 for α > 2 and this is consistent with Martinho and Batteen (2006). It may be noted that for z < -150 m the relative errors when using the POM-method are larger than 1 already for α ∼ 1 or sp ∼ 0.1. However, for z < -150 m the real density gradients and the absolute values of the errors are small. The results given in Fig. 3d

suggest that the 0.2 threshold may be increased when applying the SOPT-method. Grid convergence for the SOPT-method is demonstrated in Berntsen (2011) for parameters relevant to regional scale studies. The exact solution for Tg , see Eq. (13), has been used to investigate convergence of the methods considered here. The results for gentle slopes (small values of α ) confirm the main results from Berntsen (2011) on second order grid convergence. Here the results from a convergence study with a very steep topography, α = 9, are given, see Fig. 5. For

(a)

(b)

(c)

(d)

(e)

(f)

Fig. 5. The relative errors of the approximations to the geostrophic transport Tg given by Eq. (13) for α = 9 for the three methods considered. The errors are given as function of x in the left panel. In these experiments z is constant and equal to 0.009765625 m. The errors are given as functions of z in the right panel. In these experiments x is constant and equal to 1.5625 m. Results for β = 0.01 are given in the top panel, results for β = 0.1 are given in the middle panel and results for β = 1.0 are given in the lower panel. The black dashed lines scale with x2 and z2 , respectively.

J. Berntsen et al. / Ocean Modelling 92 (2015) 42–55

47

Fig. 6. Relative errors in the approximations to Tg (α ), see Eq. (13), as functions of the real slope α and the reduced slope α  . The relative errors in Tg (α  ) are given in a), and the errors for the numerical approximations Tg−POM (α  ), Tg−SOPT (α  ), and Tg−z (α  ) produced with the POM-method, the SOPT-method and the z-coordinate method are given in (b)–(d), respectively. In these experiments β = 0.01, x = 12.5 m, and z = 1 m.

this case, Lx is set to 25 m due to the requirement Lx > H0 /α for the solution (13). In one set of experiments, the vertical grid size z is held constant and equal to 0.009765625 m, and the horizontal grid size is varied from x = 25 m to 0.78125 m by repeated division by 2. The grid size x is chosen such that there will be 2p sub-intervals of width

x over the domain [−Lx , Lx ] for p = 1...6. The partial derivatives are computed using central differences and the integrals are computed with the repeated trapezoidal rule over sub-domains of width x and height z . In another set of experiments, the horizontal grid size is held constant and equal to x = 1.5625 m and the vertical grid size is varied from 10 m to 0.078125 m by repeated division by 2. During downwelling or upwelling situations the density interface tilt β may become large, and to explore the internal pressure gradient errors for such cases, the experiment is repeated for β = 0.01, 0.1, and 1.0. The results confirm that also for this steep topography case, the convergence rate is second order for the three methods considered. It may be noted that for a fixed horizontal grid size x, the errors may increase when the vertical grid size z is further decreased, see Figs. 5b, d, and f. The reason is that for small values of z, the errors may be dominated by the errors due to the choice of x. Having the error expansions (3) and (4) in mind, we may notice that as z or δσ is reduced a fortuitous cancellation of errors may occur for specific values of δσ . With further reduction of z or δσ this cancellation does not occur and the errors may again increase. The slope parameter sp will be reduced with x and the threshold value 0.2 is achieved for xthresh = 11.11 m ( xthresh = 0.2 ∗ 2H0 /α ). From Fig. 5a we notice that the relative errors in the approximations to Tg produced with the POM-method are larger than 1 for x larger than approximately 10 m, which is consistent with Martinho and Batteen (2006). For larger values of β , see Figs. 5c and e, the geostrophic transports increase and the relative errors in the approximations from the POM-method are reduced. For the other two meth-

ods, the relative errors are smaller than 1 for the grid sizes considered. By multiplying β with 10, the exact value of Tg becomes approximately 10 times larger. The errors for the SOPT-method and the z-coordinate method also become approximately 10 times larger, making the relative errors for these two methods robust to the change in β , see Fig. 5. The absolute values of the errors for the POM-method are robust to the change in β . This means that the relative errors in the results produced by this method are reduced by approximately a factor 10. This also indicates that to leading order the errors for the POM-method are not dependent on ∂ρ /∂ x, but rather on ∂ρ /∂ z which is consistent with Eq. (3). Statement (13) shows that the geostrophic transports are functions of the bottom slope α . If the real slope is α and the topography is smoothed such that the modified slope becomes α  , there will be an associated error. In Fig. 6a the relative errors |(Tg (α ) − Tg (α  ))/Tg (α )| are given for α = [0, 9], α  = [0, α ], and β = 0.01. For this case, the relative errors are smaller than 1% except for large values of α and much smaller values of α  . The numerical approximations to Tg (α  ) produced with method M will be denoted Tg−M (α  ) and are computed for all three methods considered with x = 12.5 m and z = 1.0 m. Grid sizes of this order of magnitude may for instance be used in local area fjord studies of dispersal from fish farms, see Ali et al. (2013). The relative errors |(Tg (α ) − Tg−M (α  ))/Tg (α )| are given in Figs. 6b to d. When considering |Tg (α ) − Tg−M (α  )| to be a measure of the error in Tg−M (α  ), we should have in mind that there are two sources of errors: one is due to the slope reduction from α to α  and the other is due to the internal pressure gradient errors associated with the numerical approximation. When using the POM-method, the errors are large even for α  = α and the errors are dominated by the internal pressure gradient errors. The errors may be reduced by topography smoothing for all

48

J. Berntsen et al. / Ocean Modelling 92 (2015) 42–55

values of α . However, as the reduced slope α  tends to zero for a given α , the errors may again increase. It may accordingly be possible to identify an α  ∈ (0, α ) that gives the minimal error when using the POM-method. The SOPT-method is calibrated such that the errors are zero for β = 0. The errors for the SOPT-method are very similar to the errors in Tg (α  ), see Figs. 6a and c, suggesting that for a small tilt β = 0.01, the internal pressure gradient errors associated with this method are small. For the larger values of α it may, however, still be beneficial to smooth the topography slightly. The numerical errors in the results from the z-coordinate method are due to the steps in the topography. However, the total errors are still smaller than 1 % for most values of α and α  considered. The errors may be reduced by reducing α  because the steps in the numerical grid will be reduced accordingly. When comparing the plots in Fig. 6, it becomes clear that when considering to smooth the topography to reduce the IPG errors one should have the properties of the chosen IPG method in mind. In particular, the results show that it is almost always beneficial to apply topography smoothing when using the POM-method. In contrast, it is

not beneficial to apply topography smoothing when using the SOPTmethod except in areas with very steep topography. 3. Two-column test case with two-layer tanh density profile In fjord or coastal studies, two layer systems with a less dense surface layer and a deeper and denser layer with an interface at depth zI may be represented with a tanh density profile

ρ (x, z ) = ρ0 − ρ tanh((z − zI + β x )/Lz ) .

(14)

The configuration will be as depicted in Fig. 2 and the density interface may have a tilt β . In the discussion, the present case with tanh density profile will be referred to as test case 2. The density gradients are given by

∂ρ

ρβ 1 (x, z ) = − . ∂x Lz cosh2 ((z − zI + β x )/Lz )

(15)

For this case, the solution (15) is used to quantify the errors in the numerical approximations produced with the three methods

Fig. 7. Errors in the approximations to ∂ρ /∂ x, see statement (15), at x = 0 as functions of z and α . In these experiments β = 0.01, x = 50 m, and z = 1.25 m. In the top panel log10 of the errors for the POM-method are given, in the middle panel the values for the SOPT-method, and in the lower panel values for the z-coordinate model are given. Absolute errors are given to the left and relative errors to the right.

J. Berntsen et al. / Ocean Modelling 92 (2015) 42–55

considered. As for test case 1, H0 = 250 m and x = 50 m. The density surface has a tilt β = 1/100. There are 200 equidistant vertical layers. Furthermore, ρ0 = 1020.0 kg/m3 , ρ = 1.0 kg/m3 , the vertical scale Lz = 50 m, and zI = −50 m. The exact solution and the approximations are computed at x = 0 m for z ∈ [−250 m, 0 m] and for α ∈ [0, 9]. The errors are given in Fig. 7. When comparing Fig. 3 for test case 1 and Fig. 7, the main impression is that the errors for each method considered as function of the slope and vertical coordinate are similar. There are, however,

49

minima in the errors for test case 2 associated with the density interface at z = −50 m that are not found in Fig. 3. These minima are due to changes in the sign of the errors across the density interface. For the POM-method and the SOPT-method, see Figs. 7a to d, the errors are negative above the near surface density interface and positive below it. The errors are again negative below the deeper minima lines. For the z-coordinate method, the errors are positive near the surface, negative in an intermediate layer down to approximately z = −80 m and again positive down to the line of discontinuity.

(a)

(b)

(c)

(d)

(e)

(f)

Fig. 8. The relative errors of the approximations to the geostrophic transport Tg given by Eq. (17) for α = 9 for the three methods considered. The errors are given as function of x in the left panel. In these experiments z is constant and equal to 0.009765625 m. The errors are given as functions of z in the right panel. In these experiments x is constant and equal to 1.5625 m. Results for β = 0.01 are given in the top panel, results for β = 0.1 are given in the middle panel and results for β = 1.0 are given in the lower panel. The black dashed lines scale with x2 and z2 , respectively.

50

J. Berntsen et al. / Ocean Modelling 92 (2015) 42–55

Fig. 9. Relative errors in the approximations to Tg (α ), see Eq. (17), as functions of the real slope α and the reduced slope α  . The relative errors in Tg (α  ) are given in a), and the errors for the numerical approximations Tg−POM (α  ), Tg−SOPT (α  ), and Tg−z (α  ) produced with the POM-method, the SOPT-method and the z-coordinate method are given in (b),–(d), respectively. In these experiments β = 0.01, x = 12.5 m, and z = 1 m.

The results for the POM-method, see Figs. 7a and b, are also for this case supporting that the slope parameter sp should be less than 0.2. For the SOPT-method, the relative errors are greater than unity below the interface for large values of α . We may, however, have in mind that well below the interface the density gradients are small. Following the same procedure as for test case 1, a solution for the equivalent geostrophic velocities vg is derived

vg (x, z ) = −

ρβ g (tanh((−zI + β x )/Lz )−tanh((z−zI + β x )/Lz )) . f ρ0 (16)

Eq. (16) is inserted in Eq. (12) to give an expression for the exact geostrophic transport Tg



ρβ g H0 Lz Tg = − (ln(cosh((−zI + β Lx )/Lz )) ρ0 f β

− ln(cosh((−zI − β Lx )/Lz ))) α L x Lz − (ln(cosh((−zI +β Lx )/Lz ))+ln(cosh((−zI −β Lx )/Lz )))

β

+

α Lz β



− Lz



Lx

−Lx Lx −Lx

ln(cosh((−zI + β x )/Lz ))dx

[ln(cosh((−zI + β x )/Lz ))

− ln(cosh((−H0 − zI + (α + β )x )/Lz ))]dx .

(17)

The one-dimensional integrals from −Lx to Lx in the equation above are estimated numerically to machine precision using the repeated trapezoidal rule. The convergence studies performed for test case 1 and α = 1 and α = 9 are repeated for test case 2 and the results are given in Fig. 8.

The results are consistent with the main findings from test case 1. The convergence rate is second order. As for the test case with an exponential density profile, the errors may increase when the vertical grid size z is further decreased, see Figs. 8b, d, and f. The explanation is the same as for the previous test case and it is given in the discussion of Fig. 5. Cancellation of errors may also occur in the integration of Tg , see Eq. (12), and this may explain that the errors for the POM-method decrease when increasing x from 12.5 m to 25 m, see the left panel of Fig. 8. The errors in the results from the SOPT-method and the z-coordinate method are generally smaller the errors in the results from the POM-method. For large values of the density interface tilt β , the benefits from using the SOPT-method or the z-coordinate method rather than the POM-method are gradually reduced. The explanation for the SOPT-method is that the benefits by using a method that is perfectly balanced for β = 0 are gradually reduced as β increases. The explanation for the z-coordinate method is that the errors due to interactions between the tilting density interface and steps in the topography increase with β . The exercise on topography smoothing performed for test case 1 is repeated for the present test case and the results are given in Fig. 9. Also for this case the exact values of Tg are robust to the slope, see Fig. 9a. The errors in the estimates of Tg produced by the POMmethod may be substantially reduced by topography smoothing and as for test case 1 there is an optimal level of topography smoothing. The errors in the estimates of Tg produced by the SOPT-method and the z-cordinate method are generally small, and except for large values of α there is less need for topography smoothing. 4. The Hardanger fjord case The Bergen Ocean Model (BOM), see Avlesen et al. (2001), Thiem et al. (2006), Berntsen et al. (2009), is set up for the Hardanger fjord, see Fig. 10. The BOM is a σ -coordinate ocean model based on the

J. Berntsen et al. / Ocean Modelling 92 (2015) 42–55

51

0

Depth [m]

−100 −200

N=0

−300

N=1 N=5

−400

N = 20

−500 −600 −700 −800 −900 0

Fig. 10. The inner part of the Hardanger fjord. The colors indicate depth in m.

1000

1500

2000

2500

Distance [m] Fig. 12. A cross section of the Hardanger fjord topography near Utne, see the line across the fjord in Fig. 10, for topographies that are Shapiro filtered N times, see the legend.

1

0.8

ing of the flow. The effects on the bottom profile over a section across the fjord near Utne, see the black line in Fig. 10, are demonstrated in Fig. 12. Repeated use of the filter basically moves mass from the shallow parts of the fjord to the deeper parts to reduce sp . This means that in studies with focus on processes near land, where for instance the fish farms are located, repeated use of a Shapiro filter may be questionable. From Fig. 12 the U-shape of a typical fjord cross section is also demonstrated, and even after 20 times use of the filter the Ushape is still very clear. There are 1123 × 1136 horizontal grid cells and 80 equidistant vertical σ -layers. The initial density profile is given by

N=0 N=1

0.6 sp

N=5

F

N = 20 0.4

0.2

0 0

500

ρ (z ) = 1025.0 kg/m3 − 5.0kg/m3 ez/25 m . 0.05

0.1 Sp

0.15

0.2

Fig. 11. The probability density function Fsp for bottom topographies that are Shapiro filtered N times, see the legend.

same governing equations as the POM. The differences in numerical techniques are explained in Berntsen (2011). The discretization is based on the C-grid, and a version with the SOPT-method included may be downloaded from http://www.math.uib.no/BOM/. The model has been applied in numerical studies of dispersal from fish farms (Thiem (2013)). In the present paper, the focus is on possible artifacts due to IPG errors in numerical studies of flow and processes in fjords. The horizontal grid size is 50 m, and the bottom topography is interpolated on to the computational grid using data from the Norwegian hydrographic service. At each velocity point in the C-grid, the slope parameter sp is computed. The probability that sp in a velocity point is larger than a given value x is denoted Fsp (x) (Fsp (x ) = Prob(s p ≥ x )). The probability density function Fsp (x ) is given in Fig. 11. This analysis of the interpolated data shows that the sp value is larger than the threshold value 0.2 in approximately 5% of the velocity points. The three point Shapiro filter (Shapiro (1975)) f i = 0.25 × ( fi−1 + 2 fi + fi+1 ) is applied in each grid cell both in x- and y-direction to reduce the sp values. After filtering the topography once, the sp values are larger than the threshold value 0.2 in approximately 3% of the velocity points. By applying the filter five times, most values of sp become smaller than 0.2, and after twenty iterations with the filter almost all values of sp are smaller than 0.1, see Fig. 11. The filtering will affect the bottom topography, and thereby the topographic steer-

(18)

The Hardanger fjord hydrography and circulation is strongly influenced by fresh water runoff. There is an annual variation of the stratification through the fjord system. The vertical density profile, however, often resembles an exponential profile which is used for the present tests. Because of the strong stratification the density 3 difference ρ is set to 5.0 kg/m . It may be noted that even if the initial density profile is a function of the vertical coordinate z only, the density field will evolve in space and time due to numerical errors. All boundaries are closed and there is no flow through the boundaries. Initially all components of the velocity are set to zero. There is no external forcing, and with the initial stratification (18), the true solution is no flow. The model is mode split (Berntsen (2000)), and the threedimensional time step t is 1.0 s and there are 30 two-dimensional time steps per three-dimensional step. The time stepping is done with a predictor-corrector method with the leapfrog method as the predictor and the fully implicit method is used as corrector in the present set of experiments. The horizontal and vertical diffusivities acting on the density field are set to zero. The horizontal viscosity is set to 10 m2 /s, and the vertical viscosity is set to 1 × 10−3 m2 /s. The density field is allowed to be advected with the velocity field generated by the errors in the internal pressure gradients and the calculations are accordingly prognostic, see Mellor et al. (1994). The bottom stress vector τb (x, y, t ) is specified by

 τb (x, y, t ) = ρ0CD u2b + v2b ub (x, y, t )

(19)

where ρ 0 is the reference density that is set to 1025 kg /m3 , ub = (ub , vb ), where ub and vb are the velocity components in the bottom most grid cell in the x-direction and the y-direction, respectively. The

52

J. Berntsen et al. / Ocean Modelling 92 (2015) 42–55

drag coefficient CD is given by



CD = max 0.0025,

κ2 (ln(zb /z0 ))2

(20)

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 z0 is set to 0.001 m. The model is started for five days with the standard POM-method and the SOPT-method for estimating the internal pressure gradients. For each method, four versions of the bottom topography are applied. The unfiltered topography, the one time filtered bottom matrix, the five times filtered bottom matrix, and the 20 times filtered bottom matrix. It may be mentioned that the use of z-level methods for the estimation of density gradients will involve vertical interpolation since we have a time dependent density field. In particular, there are problems related to the estimation of density gradients near the bottom in areas with steep topography, see for instance Ciappa (2006), when using z-level methods. We have therefore not included tests of such methods for the Hardanger fjord case. In the experiments, estimates of the volume averaged kinetic energy

Ekin =

1 2V



V

Table 1 Time to reach Umax . In the table the times to reach the threshold value Umax = 0.01 m/s1 are given for the experiments with the POM-method and the SOPTmethod using topographies that are Shapiro filtered N times. The value ∗ indicates that the threshold value is not reached after five days simulation.

(u2 + v2 )dV ,

(21)

and maximum values of the horizontal velocity components

Umax = max(u, v ) ,

(22)

are recorded. In these equations, u is the velocity component in the xdirection, v is the velocity component in the y-direction, V is the total volume of the wet cells in the computational domain. The maximum value in Eq. (22) is taken over all wet velocity points. As demonstrated for the seamount case in Beckmann and Haidvogel (1993) and Mellor et al. (1998), there may be growing modes in experiments like the present. With no flow initially, there is no initial friction to balance the growth. The friction will increase with the erroneous flow field, and depending on the magnitude of the horizontal and vertical viscosities, the growth may be balanced by friction. If not, the velocities will grow towards infinity. The times to reach a threshold value of Umax are given in Table 1. The numerical experiments based on the unfiltered topography and the POM-method were stopped when Umax exceeded 10 m/s1 after 10 s. Time series of Ekin and Umax for the experiment with the

N

POM-method

SOPT-method

0 1 5 20

2s 4s 32 s 116 s

6.6 hours ∗ ∗ ∗

unfiltered topography and the SOPT-method are given in Fig. 13. For this case with ρ = ρ (z ) initially, the SOPT-method is perfectly balanced in the sense that all density and pressure gradients should be zero given exact arithmetic. The σ SOPT values are computed from σ x Eq. (2) using the division σ SOPT = ( δxδρx )/( 1 δδxxH δσδσρ ). With exact σ

H

δx H δσ ρ should accordingly be zero. Howarithmetic δxδρx − σ δ x δσ H ever, the computations are performed with double precision FORTRAN and the relative errors are of the order 10−16 in the computations of the density gradients with the SOPT-method. The round-off errors in the density gradients accelerate slowly the velocity field, and the maximum velocities and kinetic energies increase accordingly. The time series of Ekin and Umax for the experiments with smoothed topographies are given in Fig. 14. For the SOPT-method, the velocities after five days are very small. Filtering the topography once removes the 2 x and 2 t noise and stops the associated growth of the errors. Further filtering has only minor effects on Ekin and Umax when using the SOPT-method in the present set of experiments. It may be noticed from the time series of Umax for the SOPT-method and the one time filtered topography that Umax increases slightly around day one. This is due to a slower growth of larger scale modes. The energy of these modes is, however, low, see the time series of Ekin , and further filtering of the topography stops the growth of these modes. In the experiments with the POM-method, the growth rate is reduced with repeated filtering of the topography, see Table 1. After spin-up, further acceleration of the flow is balanced by the Coriolis force and friction. However, the errors in Ekin and Umax in the experiments with the POM-method are still orders of magnitude larger than the corresponding results produced with the SOPT-method. SOPT

x

0

10

0

10 −5

10

−10

[m/s]

−5

−15

10

max

10

−20

10

U

E

kin

2

[m /s2]

10

−10

10

−25

10

−30

10

−15

10

−35

10

0

1

2

3

4

5

0

1

2

3

Time [days]

Time [days]

(a)

(b)

Fig. 13. Time series of (a) Ekin and (b) Umax for the experiment with the unfiltered topography and the SOPT-method.

4

5

J. Berntsen et al. / Ocean Modelling 92 (2015) 42–55

53

0

10

0

10 −5

10

−10

SOPT−1

SOPT−1

−5

POM−1

SOPT−5 SOPT−20 POM−1

−10

POM−5

10

POM−5

−25

10

10

max

−20

10

SOPT−20

[m/s]

SOPT−5

−15

10

U

E

kin

2

[m /s2]

10

POM−20

POM−20

−30

10

−15

10

−35

10

0

1

2

3

4

5

0

1

2

3

Time [days]

Time [days]

(a)

(b)

4

5

Fig. 14. Time series of (a) Ekin and (b) Umax for the SOPT-method (solid lines) and the POM-method (dashed lines). The number of times the topography has been filtered in each experiment is given in the legends. 0

17

−1e −−21 e7− 1

depth [m]

depth [m]

−200

−300

−1

U (ci=1E−17 m s )

−400 −500

−300

−500 −600

−700

−700

0

200

400

600

800

1000

1200

1400

1600

1800

−800

2000

−1

U (ci=0.01 m s )

−400

−600

−800

.01 0432 −.0 .0 −−0−0050.0 −0.

−100

−200

−0. .07 34561 .0 .0 2 0.0 1.0 −−0 0−0 −0−.0

−−0−0− .0 0.0 5402 .0 0. 3.01 −0

e−

−3

.01

7

−2e− 17−17 −1e

−0

0 −100

0

200

400

600

distance [m]

800

(a)

1800

1400

1600

1800

2000

−1e−17

−1

−500

−300

−500 −600

−700

−700

400

600

800

1000

1200

1400

1600

1800

2000

distance [m]

(c)

−1

W (ci=0.005 m s )

−400

−600

200

.01 −0.00 5

−1e−17

−1e−17

−0

−200

depth [m]

depth [m]

−100

W (ci=1E−17 m s )

0

1600

(b)

−400

−800

1400

.0 1505 02 .0 −−00−.0

−200 −300

1200

0

−1e−17

0 −100

1000

distance [m]

−800

0

200

400

600

800

1000

1200

2000

distance [m]

(d)

Fig. 15. The horizontal velocity component in the x-direction (into the fjord) and the vertical velocity component across the section of the Hardanger fjord near Utne, see the line across the fjord in Fig. 10, after five days. The horizontal velocity components in the x-direction are given in the upper panel and the vertical velocity components in the lower panel. The results produced with the SOPT-method are given to the left and the results from the POM-method to the right. The results are for the experiments with a one time filtered topography.

To illustrate the magnitude and nature of the errors, the along fjord and the vertical velocity components from the experiments with a one time filtered topography over a cross-section of the fjord near Utne after five days are given in Fig. 15. The errors in the results for the POM-method are substantially larger than corresponding errors produced with the SOPT-method, and they are linked to the steep nearshore slopes. For the SOPT-method with filtered topographies, the velocities are too small to affect the density field significantly. In these prognostic runs with the POM-method, the density field is advectively adjusted and the adjustment is illustrated in Fig. 16. The structure of the errors in the results produced with the POM-

method over the Utne section is very similar to the errors for the test basin given in Mellor et al. (1994), see their Figs. 5 and 7. These are in Mellor et al. (1998) denoted as SEFK-errors (Sigma Errors of the First Kind). In the experiments from Mellor et al. (1994) they die out prognostically. In the present experiments with the POM-method, and in the experiment with the SOPT-method and an unfiltered topography, they do not. The errors grow to an equilibrium level that may be reduced by repeated filtering of the topography. Analysis of the SEFKerrors as in Mellor et al. (1998) may be used to investigate possible growing modes. However, we may have in mind that the test case in Mellor et al. (1994) and the subsequent analysis of SEFK-errors were

54

J. Berntsen et al. / Ocean Modelling 92 (2015) 42–55 0

−0.1

−0.1 −0.1

−100 −200

depth [m]

−300

ρ′ (ci=0.1 kg m−3)

−400 −500 −600 −700 −800

0

200

400

600

800

1000

1200

1400

1600

1800

2000

distance [m] Fig. 16. The difference between the density field produced with the POM-method and a one time filtered topography after five days and the initial density field across the section of the Hardanger fjord near Utne.

representative for oceanic scale two-dimensional cases whereas the present experiments are for a three-dimensional fjord with very complex topography.

Acknowledgement The authors thank two reviewers and the editor for constructive remarks.

5. Discussion The results from the two idealized test cases and the Hardanger fjord case support that for cases with steep topography it is very beneficial to apply a numerical method that is balanced such that the errors are zero for ρ = ρ (z ) to reduce the internal pressure gradient errors and the growth of these errors. Furthermore, the results support that the SOPT-method may produce accurate results also when the slope parameter is larger than the suggested value of 0.2. It may be especially important to have this in mind, when using terrain following ocean models in high resolution fjord studies. The terms in the error expansion, see Eq. (4), are balanced perfectly for ρ = ρ (z ). The idealized experiments with tilting iso-lines of density show that as the tilt increases the benefits from using the SOPT-method are gradually reduced. However, for small values of the density interface tilt, which is the general case, the errors in the results from the SOPT-method are substantially smaller than corresponding errors from the POM-method. For the Hardanger fjord case, the errors may be orders of magnitude smaller when applying the SOPT-method rather than the POMmethod, and less filtering of the topography is required to reduce the errors. The use of the SOPT-method accordingly facilitates experiments with more realistic topography. The vertical density profiles are in reality also functions of the horizontal coordinates and time. That is: ρ = ρ (x, y, z, t ), see Pickard (1961), Hansen and Rattray Jr. (1966), Stigebrandt (1981), Farmer and Freeland (1983). For such cases, the σ SOPT values may be calibrated using areal means

ρ B (x, y, z, t ) =

the background stratification changes through the season. Apart from this there are no extra computational costs per time step by replacing the POM-method with the SOPT-method. The 4th order method of McCalpin (1994) is included in tests for the seamount case and the Nordic Seas given in Berntsen (2011). At closed boundaries, the order of the method is reduced to second order, and for long and narrow fjords it may be less beneficial to apply higher order methods. One may, however, have in mind that the balancing of the terms of Eq. (2) may also be done if all discrete differences and averages are taken to 4th order before computing σ SOPT . This will give a 4th order balanced method and the possible benefits remain to be investigated. In this paper, the SOPT-method is tested using the BOM. The method is, however, generic to all σ -coordinate ocean models such as the Princeton Ocean Model (www.ccpo.odu.edu/POMWEB/) and the Regional Ocean Modeling System (ROMS, www.myroms.org). The idea of balancing the errors is also applicable to unstructured grid models such as the Finite Volume Community Ocean Model (FVCOM, fvcom.smast.umassd.edu/fvcom/).

1 ρ (x , y , z, t )dx dy Area R

(23)

where R for each value of z is a circle around lateral position (x, y) with radius R. The choice of radius may be related to the length scales over which the density fields changes. A few km may be a good choice for fjords. Since, the background stratification changes through the season, the σ SOPT may need to be recomputed accordingly, see the discussion on computational costs below. This areal mean method is tested for the Nordic Seas in Berntsen (2011). When using the SOPT-method, there is an initial cost involved to compute the optimal weights σ SOPT in all wet velocity points. It is necessary to compute all terms of the discrete density gradient, see Eq. (2), and a division is required to obtain σ SOPT . In long time integrations, it may be necessary to re-calibrate the weights σ SOPT as

References Adcroft, A., Hallberg, R., Harrison, M., 2008. A finite volume discretization of the pressure gradient force using analytic integration. Ocean Modell. 22, 106–113. Ali, A., Thiem, Ø., Berntsen, J., 2011. Numerical modelling of organic waste dispersion from fjord located fish farms. Ocean Dyn. 61, 977–989. doi:10.1007/s10236-0110393-8. Ali, A., Thiem, Ø., Berntsen, J., 2013. Numerical simulation of flow and aquaculture organic waste dispersion in a curved channel. Ocean Dyn. 63, 1073–1082. doi:10.1007/s10236-013-0638-9. Asplin, L., Salvanes, A., Kristoffersen, J., 1999. Non-local wind-driven fjord-coast advection and its potential effect on plankton and fish recruitment. Fisheries Oceanography 8, 255–263. Avlesen, H., Berntsen, J., Espelid, T., 2001. A convergence study of two prognostic, sigma coordinate ocean models on a density driven flow in a quadratic basin. Int. J. Numerical Methods Fluids 36, 639–657. Barnier, B., Marchesiello, P., de Miranda, A., Molines, J.-M., Coulibaly, M., 1998. A sigmacoordinate primitive equation model for studying the circulation in the South Atlantic. Part I:model configuration with error estimates. Deep-Sea Res. I 45, 543– 572. Beckmann, A., Haidvogel, D., 1993. Numerical simulation of flow around a tall isolated seamount. Part I: problem formulation and model accuracy. J. Phys. Oceanography 23, 1736–1753. Berntsen, J., 2000. Users guide for a modesplit σ -coordinate numerical ocean model. Technical Report 135. Dept. of Applied Mathematics, University of Bergen, Johs. Bruns gt.12, N-5008 Bergen, Norway.p. 48 Berntsen, J., 2002. Internal pressure errors in sigma-coordinate ocean models. J. Atmos. Oceanic Technol. 19 (9), 1403–1414. Berntsen, J., 2011. A perfectly balanced method for estimating the internal pressure gradients in σ -coordinate ocean models. Ocean Modell. 38, 85–95. doi:10.1016/j.ocemod.2011.02.006. Berntsen, J., Oey, L.-Y., 2010. Estimation of the internal pressure gradients in σ coordinate ocean models: comparison of second, fourth, and sixth order schemes. Ocean Dyn. 60, 317–330. doi:10.1007/s10236-009-0245-y. Berntsen, J., Thiem, Ø., 2007. Estimating the internal pressure gradient errors in a σ -coordinate ocean model for the Nordic Seas. Ocean Dyn. 57, 417–429. doi:10.1007/s10236-007-0118-1. Berntsen, J., Xing, J., Davies, A., 2008. Numerical studies of internal waves at a sill: sensitivity to horizontal size and subgrid scale closure. Cont. Shelf Res. 28, 1376–1393. doi:10.1016/j.csr.2008.03.029. Berntsen, J., Xing, J., Davies, A., 2009. Numerical studies of flow over a sill: sensitivity of the non-hydrostatic effects to the grid size. Ocean Dyn. 59, 1043–1059. doi:10.1007/s10236-009-0227-0. Blumberg, A., Mellor, G., 1987. A description of a three-dimensional coastal ocean circulation model. In: Heaps, N. (Ed.), Three-Dimensional Coastal Ocean Models. American Geophysical Union, pp. 1–16. Chu, P., Fan, C., 1997. Sixth-order difference scheme for sigma coordinate ocean models. J. Phys. Oceanography 27, 2064–2071. Ciappa, A., 2006. An operational comparative test of z-levels PGF schemes to reduce pressure gradient errors of the ocean model POM. Ocean Modell. 12, 80–100. Ciappa, A., 2008. A method for reducing pressure gradient errors improving the sigma coordinate stretching function: an idealized flow patterned after the Libyan nearshore region with the POM. Ocean Modell. 23, 59–72.

J. Berntsen et al. / Ocean Modelling 92 (2015) 42–55 Cummins, P., 2000. Stratified flow over topography: time-dependent comparisons between model solutions and observations. Dyn. Atmos. Oceans 33, 43–72. Eliassen, I., Heggelund, Y., Haakstad, M., 2001. A numerical study of the circulation in Saltfjorden, Saltstraumen and Skjerstadfjorden. Cont. Shelf Res. 21, 1669– 1689. Ezer, T., Arango, H., Shchepetkin, A., 2002. Developments in terrain-following ocean models: intercomparisons of numerical aspects. Ocean Modell. 4, 249– 267. Farmer, D., Freeland, H., 1983. The physical oceanography of fjords. Prog. Oceanography 12, 147–220. Fortunato, A., Baptista, A., 1996. Evaluation of horizontal gradients in sigma-coordinate shallow water models. Atmos. Ocean 34, 489–514. Griffies, S., 2004. Fundamentals of ocean climate models. Princeton University Press. Haidvogel, D., Beckmann, A., 1999. Numerical ocean circulation modeling. vol. 2 of Series on Environmental Science and Management. Imperial College Press. Haney, R., 1991. On the pressure gradient force over steep topography in sigma coordinate ocean models. J. Phys. Oceanography 21, 610–619. Hansen, D., Rattray Jr., M., 1966. New dimensions in estuary classification. Limnology and Oceanography 11, 319–326. Kantha, L., Clayson, C., 2000. Numerical models of oceans and oceanic processes. International Geophysics Series vol. 66. Academic Press. Lamb, K., 2004. On boundary-layer separation and internal wave generation at the Knight Inlet sill. In: Proceedings of the Royal Society of London A, 460, pp. 2305– 2337. Lynge, B., Berntsen, J., Gjevik, B., 2010. Numerical studies of dispersion due to tidal flow through Moskstraumen, northern Norway. Ocean Dyn. 60, 907–920. doi:10.1007/s10236-010-0309-z. Martinho, A., Batteen, M., 2006. On reducing the slope parameter in terrain-following numerical ocean models. Ocean Modell. 13, 166–175. McCalpin, J., 1994. A comparison of second-order and fourth-order pressure gradient algorithms in a σ -coordinate ocean model. Int. J. Numerical Meth. Fluids 18, 361– 383. Mellor, G., Ezer, T., Oey, L.-Y., 1994. The pressure gradient conundrum of sigma coordinate ocean models. J. Atmos. Oceanic Technol. 11, 1126–1134.

55

Mellor, G., Häkkinen, S., Ezer, T., Patchen, R., 2002. A generalization of a sigma coordinate ocean model and an intercomparison of model vertical grids. In: Ocean Forecasting: Conceptual Basis and Applications. Springer, pp. 55–72. Mellor, G., Oey, L.-Y., Ezer, T., 1998. Sigma coordinate pressure gradient errors and the seamount problem. J.Atmos. Oceanic Technol. 15, 1122–1131. Mesinger, F., Arakawa, A., 1976. Numerical methods used in atmospheric models, vol I. WMO/ICSU Joint Organizing Committee, Garp Publication Series No. 17. Myksvoll, M., Sandvik, A., Asplin, L., Sundby, S., 2014. Effects of river regulations on fjord dynamics and retention of coastal cod eggs. ICES J. Mar. Sci. 71, 943–956. doi:10.1093/icesjms/fst113. Pickard, G., 1961. Oceanographic features of inlets in the british columbia mainland coast. J. Fisheries Res. Board Can. 18, 907–998. Shapiro, R., 1975. Linear filtering. Math. Comput. 29, 1094–1097. Shchepetkin, A., McWilliams, J., 2003. A method for computing horizontal pressuregradient force in an oceanic model with a non-aligned vertical coordinate. J. Geophys. Res. 108 (C3, 3090). doi:10.1029/2001C001047. ´ M., Janekovic, ´ I., Kuzmic, ´ M., 2009. A new approach to bathymetry smoothing Sikiric, in sigma-coordinate ocean models. Ocean Modell. 29, 128–136. Slørdal, L., 1997. The pressure gradient force in sigma-co-ordinate ocean models. Int. J. Numerical Meth. Fluids 24, 987–1017. Stelling, G., Van Kester, J., 1994. On the approximation of horizontal gradients in sigma coordinates for bathymetry with steep bottom slopes. Int. J. Numerical Meth. Fluids 18, 915–935. Stigebrandt, A., 1981. A mechanism governing the estuarine circulation in deep, strongly stratified fjords. Estuarine, Coastal and Shelf Sci. 13, 197–211. Sutherland, D., MacCready, P., Banas, N., Smedstad, L., 2011. A Model Study of the Salish Sea Estuarine Circulation. J. Phys. Oceanography 41, 1125–1143. doi:10.1175/2011JPO4540.1. Thiem, Ø., 2013. Numerical study of dispersion from a fjord located fish farm when the net pen drag is included. Technical Report 31. UNI Computing, UNI research, Bergen, Norway. Thiem, Ø., Ravagnan, E., Fosså, J., Berntsen, J., 2006. Food supply mechanisms for cold water corals along the continental shelf edge. J. Mar. Syst. 60, 207–219. doi:10.1016/j.jmarsys.2005.12.004.