- Email: [email protected]
Two- and three-dimensional numerical models of internal tide generation at a continental slope
Ocean Modelling 12 (2006) 32–45 www.elsevier.com/locate/ocemod
Two- and three-dimensional numerical models of internal tide generation at a continental slope K. Katsumata
*
School of Physical, Environmental and Mathematical Sciences, University of New South Wales, Australian Defence Force Academy, Canberra, ACT 2602, Australia Received 25 October 2004; received in revised form 23 March 2005; accepted 23 March 2005 Available online 22 April 2005
Abstract Some numerical models of internal tide generation at a continental slope are two-dimensional where the along-slope variation is neglected. The energy flux carried by internal tides computed using such twodimensional models is often underestimated, compared with three-dimensional simulations of the same region, by a factor of 10 or more. The reason for this difference is investigated using both numerical and analytical models. It is shown that in numerical models, it is not the lack of the along-shelf forcing but the use of sponge or radiating conditions at the cross-shelf boundaries that leads to the severe underestimate of the offshore flux. To obtain realistic estimates of energy flux a three-dimensional model with an along-shelf scale of at least 5 internal tide wave lengths at the depth of maximum forcing is necessary. 2005 Elsevier Ltd. All rights reserved. Keywords: Internal tide generation; Numerical model; Energy flux
1. Introduction Generation of internal tides (internal waves at tidal frequencies) at continental shelves has been a dynamically interesting problem in its own right. Theoretical models have evolved from a simple *
Tel.: +61 2 6268 8290; fax: +61 2 6268 8002. E-mail address: [email protected]
1463-5003/$ - see front matter 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.ocemod.2005.03.001
K. Katsumata / Ocean Modelling 12 (2006) 32–45
33
step topography (Rattray, 1960), increasing in complexity from gentle slopes (Baines, 1973) to steep slopes (Baines, 1974; Sandstrom, 1976; Craig, 1987). These models are all two-dimensional, where the physical quantities vary in the across-shelf and vertical directions only. Although these analytical or semi-analytical models are powerful in understanding the mechanism of the internal tide generation, application to a realistic ocean is difficult because of the restrictive assumptions on the form of the bottom topography or the density stratification. Therefore numerical models have been used to examine the effects of more realistic stratification and cross-shore bathymetry. A two-dimensional model of New (1988) successfully reproduced observed amplitudes and phases of internal tides in the Bay of Biscay (Pingree and New, 1989, 1991). Craig (1988) applied a characteristics-following method (Craig, 1987) to a cross-section on the Australian North West Shelf. The vertical structure of the internal tides predicted in this model agrees well with observations. Holloway (1996) applied a primitive equation numerical model to several cross-sections on the Australian North West Shelf and found reasonable agreement in terms of wave structure, propagation direction and regions of generation and dissipation with observations. In the well-known Abyssal Recipes, Munk (1966) and Munk and Wunsch (1998) discuss the role of internal tides in maintaining the world ocean circulation, in which a key quantity is the energy converted from the barotropic surface tide to the baroclinic internal tides (Llewellyn Smith and Young, 2002) and carried away by the internal tides. Holloway (2001) using a fully threedimensional model of the Australian North West Shelf estimated the total offshore energy flux per unit length of the shelf by the M2 semi-diurnal tide to be 600–1000 W/m depending on the stratification. Similar estimates using the two-dimensional model (Holloway, 1996) are 20–300 W/m depending on the stratification and topography. A close comparison between observations and the two-dimensional model result shows a tendency of the model to underestimate the amplitude of internal tides. The tendency was also found in the model of Craig (1988) and it was speculated that the most severe assumption of the model is that of the two-dimensionality. This conspicuous enhancement of the energy flux in a three-dimensional model is also noted in a threedimensional numerical model of northern British Columbia by Cummins and Oey (1997). This three-dimensional effect, as manifested by the underestimate of the baroclinic energy flux in two-dimensional models, is the motivation of the present work. (Related subject of internal tide generation by along-slope topographic variations is a topic of recent numerical studies (Legg, 2004a,b).) In models of internal tide generation, the difference between two-dimensional and three-dimensional models can be attributed to two different processes; the difference in forcing, and the difference in response to the forcing. The tidal forcing is a vertically uniform barotropic flow which is driven by surface displacement (Baines, 1973). In reality, this surface displacement is set up by the gravitational attraction from the sun and the moon, but in numerical models the displacement is usually prescribed at the model boundary. Even when the prescribed displacement is the same, the barotropic flow is not necessarily the same because the barotropic flow depends on the bottom topography. This is most clearly demonstrated by Holloway and Merrifield (1999). In their numerical study, a seamount and a ridge are used with the same boundary forcing and it is shown that the resulting baroclinic energy flux changes by nearly an order of magnitude because the barotropic flow results in stronger vertical displacement of isopycnals for the ridge topography. In Section 2, this relationship between the forcing and the resulting flux is quantitatively examined with a realistic numerical model. It will be shown that the difference in forcing does not explain the underestimate of energy flux by two-dimensional models. In Sections 3 and 4, it will be
34
K. Katsumata / Ocean Modelling 12 (2006) 32–45
-5000
0
00
-10
-5000
00
-5
-15˚
0 00 200 -10 -
0
-10 -3 00 -4000 000
00
-2000
-4000 -50 00
-3
-50 -3 00 000
-400 0
-10˚
-40000 0 -30 -2000
-20˚
W
-50
-4 00 0
00
0 -3
00
y
x
-25˚ 110˚
115˚
120˚
125˚
130˚
Fig. 1. Model domain. Three model domains with different widths (W) are shown by three different lines. Also shown is the definition of x and y, which are approximately along-slope and across-slope directions.
shown that the effect of across-shelf sponge layers or the radiation boundary condition is the cause of the underestimate. Although the underestimate in numerical models can be explained as above, the reason for the underestimate in analytical two-dimensional models remains unclear. For example, the twodimensional analytical model by Baines (1982) predicted the flux to be 174 W/m on the Australian North West Shelf. Using more realistic topography and stratification (but two-dimensional), Craig (1988) estimated the flux to be less than 100 W/m. A possible explanation is the sensitivity of the flux to such conditions as the topographic slope a, stratification, and the forcing barotropic flow. Craig (1987) found that on a subcritical linear shelf slope and with linear stratification, the flux scales as a5. Also the model by Baines (1973) shows an order-of-magnitude difference in offshore energy flux depending on the steepness of the continental shelf. The dependence of the flux on the stratification might not be as simple, but given the nonlinear behaviour of internal waves (Maas and Lam, 1995), it is possible that a choice of a ÔrepresentativeÕ cross-section and stratification from a highly-complicated topography (e.g. see Fig. 1) and stratification, may explain the order-of-magnitude underestimate of these analytical models.
2. Numerical experiments with realistic bathymetry In this section, sensitivity of a three-dimensional model to the size of the domain is examined by varying the model width. The relation between the forcing and the resulting flux is investigated in detail. The Princeton Ocean Model (Blumberg and Mellor, 1987) was used to simulate the internal tide in the region shown in Fig. 1. The model setup is almost identical to that of Holloway (2001). The
K. Katsumata / Ocean Modelling 12 (2006) 32–45
35
full domain size excluding the sponge layer is 1540 km · 624 km, the grid size is 4 km · 4 km, and the vertical sigma grid has 51 levels. Forcing at the three open boundaries is provided by the M2 barotropic tidal elevation predicted by the global tidal model of Le Provost et al. (1998). The same initial horizontally uniform stratification as Holloway (1996) is used, that is, the average of February, March, and April salinity and temperature data taken from the Levitus atlas (Levitus, 1982). A relaxation layer exists at the outer 10 grids of the across-slope boundaries and the offshore boundaries. The model solves the barotropic (vertically averaged) and baroclinic flow components seperately. In the relaxation layer, the barotropic velocity is damped by increased horizontal viscosity of 500 m2/s, and the baroclinic velocity is relaxed gradually to zero by Martinsen and Engedahl (1987)Õs method, and the temperature and salinity are relaxed similarly to values at the outermost grids calculated by the upwind differencing scheme. The horizontal viscosity is fixed at zero in the interior. The model width W (Fig. 1) is varied from 120 km to 1540 km (excluding the relaxation layers) by changing the number of grids. The forcing is always by the prescribed tidal elevation at all three open boundaries. The models are run for 12 days and a tidal analysis is performed on the data from the last day. The energy flux carried by the internal tides is calculated and integrated along the outermost grids (but not in the relaxation layer) as Z ZZ 1 T vn p dl dz dt; T 0 where T is the semi-diurnal tidal period, vn is the baroclinic velocity normal to the boundary, p is the baroclinic component of the pressure, l is the along-boundary coordinate, and t is the time. The detail of the calculation is explained in the appendix of Holloway and Merrifield (1999). The integrated energy flux is shown in Fig. 2(a). Note that vn is the component crossing the
Flux
1.6 1.4
Flux / Width
1200
x-flux y-flux
x-flux y-flux
1000 800
1
Flux/Width [W/m]
Integrated flux [GW]
1.2
0.8 0.6 0.4
600 400 200
0.2 0
0 -0.2
-200 0
(a)
500
1000 1500 W [km]
2000
0
(b)
500
1000 1500 2000 W [km]
Fig. 2. (a) Energy flux across the boundaries integrated along the outermost grids as a function of the model width. y-flux is the offshore flux. (b) The integrated energy flux divided by the length of the shelf as a function of the model width.
36
K. Katsumata / Ocean Modelling 12 (2006) 32–45
boundary so that x-flux in Fig. 2 is the flux evaluated on the across-slope boundaries and y-flux is the offshore flux. The along-shelf flux (x-flux) does not bear a simple relationship with the model width, but the across-shelf flux (y-flux) shows an almost linear increase for W > 500 km. For W < 500 km, this linear relation does not hold, which is emphasized in Fig. 2(b), where the flux is normalised by the model width W. When the model is wider, the flux per unit shelf length is approximately 1000 W/m, but when the model region is narrower than approximately 500 km, the outgoing flux is much reduced. This situation is similar to the underestimate by two-dimensional models mentioned in the introduction, although the topography and the barotropic tidal forcing are not strictly two-dimensional for small W. For example, if the 120-km wide model is used to simulate the internal tides at the middle cross-section of the shelf, the model would predict the outgoing offshore flux of 130 W/m, which is much smaller than the offshore flux predicted by wider (W > 500 km) models (approximately 1000 W/m). It is possible that the centre of the model domain just so happens to coincide with a region of weak generation due to local details of the barotropic forcing and bathymetry. This possibility will be examined in Section 3, but for now it is assumed that the narrow-width limit represents the two-dimensional case. It is conjectured that a strong along-shelf barotropic tidal flow, combined with along-shelf topographic variation, produces additional forcing of the internal tide (Legg, 2004a,b). Baines (1973) defines the internal tide forcing as og o u ov z þ ; ð1Þ F ¼ N2 ot ox oy where N is the buoyancy frequency, g is the surface displacement, and the overbar denotes vertical averaging. This term appears in the density conservation equation after separating the velocity into the barotropic tidal field of an equivalent unstratified ocean and the internal tidal field (Baines, 1973, also see Eq. (10)). In the two-dimensional model, the ou=ox term is absent. An interpretation of the forcing can be given in the following way. The potential energy change E of a fluid parcel of unit volume by a barotropic tidal flow in a stratified ocean is approximately dq0 2 ð2Þ f; E¼g dz where f is the vertical displacement of the fluid parcel caused by the barotropic flow, and q0(z) is the stratification of the undisturbed ocean. The rate of work to cause this energy increase is dE dq0 df ¼g ð3Þ 2f ¼ 2fq0 F ; dt dt dz where F is the forcing defined in (1). Using this interpretation, it is possible to provide a measure of the internal tide amplitude in the model. By integrating (3) in space and time and assuming that all the perturbation energy is radiated as the internal tide, (i.e. dE dx dy dz/dt = Flux dl dz) one can define RR Flux dl dz RRR ; ð4Þ f0 ¼ 2q0 Forcing dx dy dz where all terms are time averaged. Alternatively, f0 can be regarded as the ‘‘efficiency’’ of the model as an internal tide generator, for it yields the output flux when multiplied by the input forcing.
K. Katsumata / Ocean Modelling 12 (2006) 32–45 ζ0
Integrated forcing
0.9
37
1.2
Amplitude [m]
Forcing [x106 m4s-3]
0.9 0.6
0.6
0.3 0.3
0
0 0
(a)
500
1000 1500 W [km]
2000
0
(b)
500
1000 1500 W [km]
2000
Fig. 3. (a) Forcing (1) integrated in the model domain as a function of the model width. (b) ‘‘Efficiency’’ or the total outgoing flux divided by the total forcing (see (4)) as a function of the model width.
Fig. 3 shows the forcing integrated over the model region and the corresponding f0 as functions of the model width. It is seen in Fig. 3(a) that the forcing increases almost linearly with the model width W, and as a result the average tidal amplitude, or the efficiency, f0, increases for W < 800 km but saturates for W > 800 km. From Fig. 3(b), it is implied that the narrow-width models are not as ÔefficientÕ as the wide models in generating the internal tide. This result suggests that the reason for the underestimate in the narrow-width (and two-dimensional) models is twofold; one is the reduced forcing and the other is the reduced efficiency. A closer look at Figs. 2(b) and 3(b) suggests that there are two different dynamical regimes in the transition between the two-dimensional and three-dimensional models. One is between W = 0 km and W 500 km, where the normalised flux increases with the model width, and the other is W > 500 km, where the normalised flux is approximately constant. In the large scale regime, the difference in the energy flux is mostly explained by the difference in the forcing because f0 is almost constant and the integrated flux (Flux) is almost proportional to the integrated forcing (Forcing). In the small scale regime, however, not only the normalised flux, but also f0 increases with the model width. This small scale regime is further investigated in the next two sections.
3. Numerical experiments with idealised bathymetry In the experiments described in the previous section, when the model width is narrower than approximately 500 km, the resulting flux is not proportional to the forcing. In order to ascertain that this phenomenon is not caused by the irregular topography and the tidal flow of the particular Australian North West Shelf simulation, an idealised situation is considered in this section.
38
K. Katsumata / Ocean Modelling 12 (2006) 32–45
An idealised shelf topography y ð5Þ H ðyÞ ¼ D þ A tanh b with D = 128 m, A = 16 m, b = 2 km is used. The model grid size is fixed at 0.8 km · 0.8 km, and the across-slope size of the model domain is 51.2 km (64 grids) including 8 km (10 grids) wide relaxation layers on the across-slope boundaries and the on/offshore boundaries. The relaxation scheme is identical to the one used in the realistic model in Section 2. Tidal forcing is given by the surface elevation with an amplitude of 1 cm at the M2 tidal frequency x and with opposite phases at the two long-shelf boundaries. Fig. 4 shows the forcing, which is almost perfectly along-slope and confined near the slope. Constant stratification of N = 1.07 · 102 s1 is set up with the salinity while the temperature is fixed constant. The Coriolis parameter is fixed at the value of 30 north. The model was run for four days from quiescent and horizontally uniform initial conditions and results from the tidal analysis from the last day are examined. The same experiment as in Section 2 was carried out and the result is shown in Fig. 5. The horizontal length scale is nondimensionalised by the internal tide characteristic wave length L0 = 2D/c (approximately 23 km here), and the forcing and the flux are nondimensionalised by ð6Þ F 0 ¼ N 2 VD=L0 and G0 ¼ q0 A2 V 2 N 2 H =ð2xL0 Þ; ð7Þ 1 respectively, where V is the typical across shelf barotropic flow speed (here 3.5 cm s ), and c is the characteristic slope (here 1.1 · 102) of the internal tide defined by x2 f 2 c2 ¼ 2 . ð8Þ N x2
y
1
1 5 5 1
0
1.0
-1 -1
0 x
1
Fig. 4. Forcing field in the simple numerical model when the model width is 51.2 km. The ellipses show the nondimensionalised barotropic forcing flow. The scale of ellipses is shown on the bottom right. The contours show the nondimensionalised strength of the forcing defined in (1) with the contour interval of 2 unit forcing. The size of the ellipses and the forcing are nondimensionalised by V and DF0, respectively (see text for definition). The shaded area is the relaxation layers.
K. Katsumata / Ocean Modelling 12 (2006) 32–45 1
39
0.8
0.8 0.6
× G0
× DF0
0.6 0.4
0.4 0.2 0.2
0
0 0
(a)
5 10 15 Model width (× L0)
20
0
(b)
5 10 15 Model width (× L0)
20
Fig. 5. (a) Across-slope offshore energy flux as a function of the along-slope model width. The flux is integrated over the on- and off-shore boundaries and normalised per unit length of the shelf. (b) Forcing as a function of the alongslope model width. The forcing is integrated in depth and averaged in area except the relaxation layers.
The forcing scale (6) is taken from (1), and the flux scale (7) is based on the dimensional form of the analytic solution (17) discussed in the following section. It is demonstrated that even though the forcing per unit area is almost constant, the resulting flux drastically decreases when the model width is less than approximately 5 internal tide wave lengths. Fig. 6 is an example of the model output, showing the vertically integrated and tidally averaged energy flux when the model width is 51.2 km. The flux field shows strong along-slope variability,
y
1
1.1 1 0.9
0
1.0
-1 -1
0
1
x Fig. 6. Vertically integrated and tidally averaged energy flux when the model width is 51.2 km. The shaded zone indicates the relaxation layers. The contour is the depth with contour interval of 0.1. All the quantities are nondimensional.
40
K. Katsumata / Ocean Modelling 12 (2006) 32–45
not seen in the forcing distribution shown in Fig. 4. The fact that the offshore flux is the strongest near the middle (x 0) suggests that this along-slope variability comes from the across-shelf relaxation layers. Having tracked down the problem to a simple setting, we can now use an even simpler semianalytical model to understand the reason for the underestimate in the narrow-width models.
4. A simple model of three-dimensional internal tide generation 4.1. Formulation Governing equations of an inviscid and linear internal wave field with the Boussinesq and hydrostatic approximations are p ot u fv ¼ ox ; 0 q p ot v þ fu ¼ oy ; 0 q p qg ð9Þ 0 ¼ oz ; 0 q 0 q dq dq ot q þ w 0 ¼ W 0 ; dz dz ox u þ oy v þ oz w ¼ 0; where (u, v, w) are the velocities associated with the internal waves, p is the pressure perturbation 0 is a constant reference density, q0(z) is the unperturbed density, and q is caused by the waves, q the density perturbation. W is the vertical velocity associated with the barotropic tide defined by og o u ov z þ W ¼ ð10Þ ot ox oy (cf. Eq. (1)). For simplicity, we assume the buoyancy frequency N is constant. The boundary conditions in the vertical direction are w ¼ 0 at z ¼ 0; w ¼ uox H voy H
at z ¼ H.
ð11Þ
The other boundary condition is the radiation condition for energy at open horizontal boundaries. These equations and the boundary conditions are identical to those used by Baines (1973). The topography H(x, y) is prescribed as (5) and the barotropic tide ðuðx; y; tÞ; vðx; y; tÞÞ is given by Q ð12Þ u ¼ 0; v ¼ cos xt; H where Q is a constant. With the rigid-lid approximation, (10) becomes W ¼ AQ
z 1 y sech2 cos xt; b H2 b
K. Katsumata / Ocean Modelling 12 (2006) 32–45
for the bottom topography (5). In the narrow-shelf limit, b ! 0 and the forcing becomes z W ¼ 2AQ dðyÞ cos xt. ðD þ A sgn ðyÞÞ2
41
ð13Þ
Offshore (y > 0), the bottom is well approximated by D1 = D + A = constant in the narrow-shelf limit and the solution can be obtained using a vertical mode expansion. The solution for w, subject to the boundary condition (11) is given by 1 X npz ixt wn ðx; yÞ sin e . ð14Þ wðx; y; zÞ ¼ D1 n¼1 Henceforth, the physical quantities are nondimensionalised by D1 (vertical length), L0 = 2D1/c (horizontal length), and 1/x (time) unless otherwise noted. It is also assumed that all quantities have time-dependence eit. Then the formal solution for the vertical nth mode wn is Z 1 Z 1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð1Þ in2 p2 H 0 2np ðx nÞ2 þ ðy gÞ2 W n ðn; gÞ dn dg; ð15Þ wn ðx; yÞ ¼ W n 1
1
ð1Þ H0
is the Hankel function of the first kind and Wn is the vertical nth mode of the barowhere ð1Þ tropic forcing W given by (10). The radiation condition of the energy is satisfied by taking H 0 ð2Þ rather than H 0 . Essentially the same solution has been given by Llewellyn Smith and Young (2002). The physical interpretation of the solution (15) is that the internal tide response at superposition of (x, y) is the response to the local forcing, Wn, plusqthe the response from the ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð1Þ
forcing at (n, g) with the GreenÕs function H 0
2np ðx nÞ2 þ ðy gÞ2 .
Since only the offshore region y > 0 is considered, the integration in g should be limited to g > , where is an infinitesimal positive number. The experiment in Section 3 is replicated by limiting the along-shelf extent of the source within [S, S]. Thus, Z S qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð1Þ n aH 0 2np ðx nÞ2 þ y 2 dn; wn þ W n ¼ 4iðÞ np ð16Þ S
where a is the nondimensional forcing strength a ¼ DA1 D1QLx . From (16) and (9), it is possible, by numerically integrating (16), to calculate the energy flux, 2 2 2 iG oG f G oG ! Z 0 u 1 p N A Q ox oy x R ; dz ¼ R f oG oG 2 q x D D Lx v 1 1 iG oy þ x G ox p 0
ð17Þ
where R; ð Þq means the real part and the complex conjugate, respectively, and G ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R S ð1Þ H 0 2np ðx nÞ2 þ y 2 dn. Fig. 7 shows an example with the same source width as S
Fig. 6. Because of the idealised d-function forcing, the magnitude of the energy flux (17) does not explicitly depend on the vertical mode n, but it is seen that the higher modes have smaller spatial scales. Gravest 10 modes are chosen to plot Fig. 7 to compare the analytical solution with the numerical simulation with expectation that higher modes are prone to the stronger damping. This choice is based on a suggestion by St. Laurent and Garrett (2002) that internal waves with equivalent modes less than 10 are likely to radiate away. Despite the many simplifying assumptions, a
42
K. Katsumata / Ocean Modelling 12 (2006) 32–45
y
1.0
0.5
-1
0
x
1
1.0
Fig. 7. Offshore flux calculated by (17) for the solution (16). The source width is S = 1.0, which corresponds to Fig. 6.
qualitative similarity can be found between Figs. 6 and 7 such as the magnitude of the flux and weak radial refraction. The stronger damping in the numerical solution might be caused by incomplete absorption in the relaxation layers. Two obvious limits exist for (16). First, the point source limit can be obtained by replacing the pffiffiffiffiffiffiffiffiffiffiffiffiffiffi constant a with d(n) and it can be seen that all the quantities are now functions of r ¼ x2 þ y 2 only. It is also possible to show that the outgoing energy flux is constant when integrated across the surface r = constant, which means the energy flux per unit length (energy flux density) decreases as r1. Second, the two-dimensional limit is obtained by S ! 1, where it is possible to analytically integrate (16) such that Z S qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð1Þ n aH 0 2np ðx nÞ2 þ y 2 dn; ¼ 4iðÞn nei2npy ð18Þ wn þ W n ¼ 4iðÞ np S
(see e.g. Moriguchi et al., 1987). This limit can also be obtained by neglecting oxo terms in (9). In fact, these two limits give an explanation of the result of the numerical experiments in Section 3. Assuming that the relaxation layers work perfectly and no spurious reflection occurs, the narrow width limit in the numerical experiments corresponds to the point source limit above. In this limit, the energy flux per unit length of the offshore boundary decreases as r1. On the other hand, in the two-dimensional limit, the energy flux per unit length of the offshore boundary is constant for all y > 0. The dependence of the flux on the source width S is evaluated and shown in Fig. 8. The flux increases as S increases and saturates around 2S 5, in reasonable agreement with the result of the numerical experiments shown in Fig. 5. Thus it is shown that the underestimation in narrow-width models is a result of the artificial across-shelf relaxation (or radiation) boundaries.
5. Discussion An obvious remedy to avoid the effect of the across-shelf artificial boundaries is to use periodic boundary conditions in the along-shelf direction. In fact, when this boundary condition is used, the flux estimates by numerical simulations and theoretical estimates agree well (Khatiwala, 2003).
K. Katsumata / Ocean Modelling 12 (2006) 32–45
43
Offshore flux (× F0)
0.8
0.6
0.4
0.2
0 0
5
10 15 Source length 2S (× L0)
20
Fig. 8. Offshore flux integrated along y = 1.0 and normalised per unit length as a function of the source width S.
In realistic applications, however, the use of the periodic conditions is complicated by two factors. First, in most applications both the tidal forcing and the topography show significant along-slope variation (for example Holloway, 2001; Cummins and Oey, 1997), and oversimplification of such variation may result in unreliable estimates. Second, in order to accommodate the periodic boundary condition, the forcing must be in phase along the shelf, producing unrealistic barotropic tidal flow and energy flux. The experiment in Section 2 was repeated with the narrow model width of 52 km, whose forcing and topography are taken near the middle of the model region (x = W/2 in Fig. 1), but with the periodic boundary condition in the along-shelf direction. The resultant offshore flux was unrealistically small (less than 1% of the wide scale models in Fig. 2). Therefore for a realistic estimate, a fully-three-dimensional numerical simulation seems to be the appropriate method to use. In Figs. 5 and 8, it was found that the saturation of the normalised offshore energy flux occurs when the model width is approximately 5 internal wave lengths. The result shown in Fig. 2 can be interpreted similarly if an appropriate depth of the shelf is chosen. It is possible to calculate the forcing (1) from the stratification and the tidal data for the region of interest without actually
Fig. 9. Vertically integrated forcing (1) on Australian North West Shelf. The depth contours are shown every 100 m up to 1000 m depth, then every 1000 m.
44
K. Katsumata / Ocean Modelling 12 (2006) 32–45
running a numerical model. Fig. 9 shows the forcing thus calculated for the experiment in Section 2 with the full domain size. It is seen that the maximum forcing is localised along the steep shelf slope around 1000 m depth and weaker peaks near the shelf break at 200 m depth. Below the thermocline, the buoyancy frequency is approximately N = 7 · 103 s1, and the corresponding characteristic slope for the M2 tide is 2 · 102. Therefore at 1000 m depth, the internal tide wave length as defined by the distance between two surface reflections of an internal tide ray is approximately 100 km. In Fig. 2, the saturation occurs approximately at W = 500 km, which is about 5 wave lengths. Additional tests were also performed with the simple numerical model in Section 3 with along-slope sinusoidal variation of the topography such that ! y sinð2px=LÞ ; ð19Þ H ðx; yÞ ¼ D þ A tanh þ b cosh2 ðy=bÞ in place of (5). The along-slope model widths were changed as in Section 3. The results for different along-slope variations (L = 1.0L0, 1.2L0, 1.8L0, 2.3L0) were all similar to Fig. 5, that is, when the model width is 5L0, the offshore flux per unit length is more than 90% of the corresponding wide (W = 18L0) models. It is thus confirmed that to obtain realistic tidal simulation and energy fluxes, the model domain must be greater than 5 internal wave lengths.
Acknowledgments The author thanks Drs. G. Symonds, P. Craig, and M. Merrifield for helpful discussions. This work was initiated from a chat with late Dr. Peter Holloway. Comments from two anonymous reviewers helped to improve the manuscript and are gratefully acknowledged.
References Baines, P.G., 1973. The generation of internal tides by flat-bump topography. Deep-Sea Res. 20, 179–205. Baines, P.G., 1974. The generation of internal tides over steep continental slopes. Philos. Trans. R. Soc. London A 277, 27–58. Baines, P.G., 1982. On internal tide generation models. Deep-Sea Res. 29, 307–338. Blumberg, A.F., Mellor, G.L., 1987. A description of a three-dimensional coastal ocean model. In: Heaps, N. (Ed.), Three-dimensional Coastal Ocean Models, Coastal and Estuarine Science. American Geophysical Union, AGU, Washington, DC, pp. 1–16. Craig, P.D., 1987. Solutions for internal tidal generation over coastal topography. J. Mar. Res. 45, 83–105. Craig, P.D., 1988. A numerical model study of internal tides on the Australian Northwest Shelf. J. Mar. Res. 46, 59–76. Cummins, P.F., Oey, L.-Y., 1997. Simulation of barotropic and baroclinic tides off northern British Columbia. J. Phys. Oceanogr. 27, 762–781. Holloway, P.E., 1996. A numerical model of internal tides with application to the Australian North West Shelf. J. Phys. Oceanogr. 26, 21–37. Holloway, P.E., 2001. A regional model of the semidiurnal internal tide on the Australian North West Shelf. J. Geophys. Res. 106, 19625–19638. Holloway, P.E., Merrifield, M.A., 1999. Internal tide generation by seamounts, ridges, and islands. J. Geophys. Res. 104, 25937–25951.
K. Katsumata / Ocean Modelling 12 (2006) 32–45
45
Khatiwala, S., 2003. Generation of internal tides in an ocean of finite depth: Analytical and numerical calculations. Deep-Sea Res. 50, 3–21. Legg, S., 2004a. Internal tides generated on a corrugated continental slope. Part I: Cross-slope barotropic forcing. J. Phys. Oceanogr. 34, 156–173. Legg, S., 2004b. Internal tides generated on a corrugated continental slope. Part II: Along-slope barotropic forcing. J. Phys. Oceanogr. 34, 1824–1838. Le Provost, C., Lyard, F., Molines, J.M., Genco, M.L., Rabilloud, F., 1998. A hydrodynamic ocean tide model improved by assimilating a satellite altimeter-derived data set. J. Geophys. Res. 103, 5513–5529. Levitus, S., 1982. Climatological atlas of the world ocean. NOAA Prof. Rep. 13. Environ. Res. Lab., Geophys. Fluid Dyn. Lab., US Department of Commer., Princeton, NJ, 173pp. Llewellyn Smith, S.G., Young, W.R., 2002. Conversion of the barotropic tide. J. Phys. Oceanogr. 32, 1554–1566. Maas, L.R.M., Lam, F.-P.A., 1995. Geometric focusing of internal waves. J. Fluid Mech. 300, 1–41. Martinsen, E.A., Engedahl, H., 1987. Implementation and testing of a lateral boundary scheme as an open boundary condition in a barotropic ocean model. Coast. Eng. 11, 603–627. Moriguchi, S., Udagawa, K., Hitotsumatsu, S., 1987. In: Sugaku Koshiki II (mathematical Formulae II). Iwanami Shoten, Tokyo, p. 340 (in Japanese). Munk, W., 1966. Abyssal recipes. Deep-Sea Res. 13, 707–730. Munk, W., Wunsch, C., 1998. Abyssal recipes II: Energetics of tidal and wind mixing. Deep-Sea Res. 45, 1977–2010. New, A.L., 1988. Internal tidal mixing in the Bay of Biscay. Deep-Sea Res. 35, 691–709. Pingree, R.D., New, A.L., 1989. Downward propagation of internal tidal energy into the Bay of Biscay. Deep-Sea Res. 36, 735–758. Pingree, R.D., New, A.L., 1991. Abyssal penetration and bottom reflection of internal tidal energy in the Bay of Biscay. J. Phys. Oceanogr. 21, 28–39. Rattray Jr., M, 1960. On the coastal generation of internal tides. Tellus 12, 54–62. Sandstrom, H., 1976. On topographic generation and coupling of internal waves. Geophys. Fluid Dyn. 7, 231–270. St. Laurent, L., Garrett, C., 2002. The role of internal tides in mixing the deep ocean. J. Phys. Oceanogr. 32, 2882–2899.
Two- and three-dimensional numerical models of internal tide generation at a continental slope K. Katsumata
*
School of Physical, Environmental and Mathematical Sciences, University of New South Wales, Australian Defence Force Academy, Canberra, ACT 2602, Australia Received 25 October 2004; received in revised form 23 March 2005; accepted 23 March 2005 Available online 22 April 2005
Abstract Some numerical models of internal tide generation at a continental slope are two-dimensional where the along-slope variation is neglected. The energy flux carried by internal tides computed using such twodimensional models is often underestimated, compared with three-dimensional simulations of the same region, by a factor of 10 or more. The reason for this difference is investigated using both numerical and analytical models. It is shown that in numerical models, it is not the lack of the along-shelf forcing but the use of sponge or radiating conditions at the cross-shelf boundaries that leads to the severe underestimate of the offshore flux. To obtain realistic estimates of energy flux a three-dimensional model with an along-shelf scale of at least 5 internal tide wave lengths at the depth of maximum forcing is necessary. 2005 Elsevier Ltd. All rights reserved. Keywords: Internal tide generation; Numerical model; Energy flux
1. Introduction Generation of internal tides (internal waves at tidal frequencies) at continental shelves has been a dynamically interesting problem in its own right. Theoretical models have evolved from a simple *
Tel.: +61 2 6268 8290; fax: +61 2 6268 8002. E-mail address: [email protected]
1463-5003/$ - see front matter 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.ocemod.2005.03.001
K. Katsumata / Ocean Modelling 12 (2006) 32–45
33
step topography (Rattray, 1960), increasing in complexity from gentle slopes (Baines, 1973) to steep slopes (Baines, 1974; Sandstrom, 1976; Craig, 1987). These models are all two-dimensional, where the physical quantities vary in the across-shelf and vertical directions only. Although these analytical or semi-analytical models are powerful in understanding the mechanism of the internal tide generation, application to a realistic ocean is difficult because of the restrictive assumptions on the form of the bottom topography or the density stratification. Therefore numerical models have been used to examine the effects of more realistic stratification and cross-shore bathymetry. A two-dimensional model of New (1988) successfully reproduced observed amplitudes and phases of internal tides in the Bay of Biscay (Pingree and New, 1989, 1991). Craig (1988) applied a characteristics-following method (Craig, 1987) to a cross-section on the Australian North West Shelf. The vertical structure of the internal tides predicted in this model agrees well with observations. Holloway (1996) applied a primitive equation numerical model to several cross-sections on the Australian North West Shelf and found reasonable agreement in terms of wave structure, propagation direction and regions of generation and dissipation with observations. In the well-known Abyssal Recipes, Munk (1966) and Munk and Wunsch (1998) discuss the role of internal tides in maintaining the world ocean circulation, in which a key quantity is the energy converted from the barotropic surface tide to the baroclinic internal tides (Llewellyn Smith and Young, 2002) and carried away by the internal tides. Holloway (2001) using a fully threedimensional model of the Australian North West Shelf estimated the total offshore energy flux per unit length of the shelf by the M2 semi-diurnal tide to be 600–1000 W/m depending on the stratification. Similar estimates using the two-dimensional model (Holloway, 1996) are 20–300 W/m depending on the stratification and topography. A close comparison between observations and the two-dimensional model result shows a tendency of the model to underestimate the amplitude of internal tides. The tendency was also found in the model of Craig (1988) and it was speculated that the most severe assumption of the model is that of the two-dimensionality. This conspicuous enhancement of the energy flux in a three-dimensional model is also noted in a threedimensional numerical model of northern British Columbia by Cummins and Oey (1997). This three-dimensional effect, as manifested by the underestimate of the baroclinic energy flux in two-dimensional models, is the motivation of the present work. (Related subject of internal tide generation by along-slope topographic variations is a topic of recent numerical studies (Legg, 2004a,b).) In models of internal tide generation, the difference between two-dimensional and three-dimensional models can be attributed to two different processes; the difference in forcing, and the difference in response to the forcing. The tidal forcing is a vertically uniform barotropic flow which is driven by surface displacement (Baines, 1973). In reality, this surface displacement is set up by the gravitational attraction from the sun and the moon, but in numerical models the displacement is usually prescribed at the model boundary. Even when the prescribed displacement is the same, the barotropic flow is not necessarily the same because the barotropic flow depends on the bottom topography. This is most clearly demonstrated by Holloway and Merrifield (1999). In their numerical study, a seamount and a ridge are used with the same boundary forcing and it is shown that the resulting baroclinic energy flux changes by nearly an order of magnitude because the barotropic flow results in stronger vertical displacement of isopycnals for the ridge topography. In Section 2, this relationship between the forcing and the resulting flux is quantitatively examined with a realistic numerical model. It will be shown that the difference in forcing does not explain the underestimate of energy flux by two-dimensional models. In Sections 3 and 4, it will be
34
K. Katsumata / Ocean Modelling 12 (2006) 32–45
-5000
0
00
-10
-5000
00
-5
-15˚
0 00 200 -10 -
0
-10 -3 00 -4000 000
00
-2000
-4000 -50 00
-3
-50 -3 00 000
-400 0
-10˚
-40000 0 -30 -2000
-20˚
W
-50
-4 00 0
00
0 -3
00
y
x
-25˚ 110˚
115˚
120˚
125˚
130˚
Fig. 1. Model domain. Three model domains with different widths (W) are shown by three different lines. Also shown is the definition of x and y, which are approximately along-slope and across-slope directions.
shown that the effect of across-shelf sponge layers or the radiation boundary condition is the cause of the underestimate. Although the underestimate in numerical models can be explained as above, the reason for the underestimate in analytical two-dimensional models remains unclear. For example, the twodimensional analytical model by Baines (1982) predicted the flux to be 174 W/m on the Australian North West Shelf. Using more realistic topography and stratification (but two-dimensional), Craig (1988) estimated the flux to be less than 100 W/m. A possible explanation is the sensitivity of the flux to such conditions as the topographic slope a, stratification, and the forcing barotropic flow. Craig (1987) found that on a subcritical linear shelf slope and with linear stratification, the flux scales as a5. Also the model by Baines (1973) shows an order-of-magnitude difference in offshore energy flux depending on the steepness of the continental shelf. The dependence of the flux on the stratification might not be as simple, but given the nonlinear behaviour of internal waves (Maas and Lam, 1995), it is possible that a choice of a ÔrepresentativeÕ cross-section and stratification from a highly-complicated topography (e.g. see Fig. 1) and stratification, may explain the order-of-magnitude underestimate of these analytical models.
2. Numerical experiments with realistic bathymetry In this section, sensitivity of a three-dimensional model to the size of the domain is examined by varying the model width. The relation between the forcing and the resulting flux is investigated in detail. The Princeton Ocean Model (Blumberg and Mellor, 1987) was used to simulate the internal tide in the region shown in Fig. 1. The model setup is almost identical to that of Holloway (2001). The
K. Katsumata / Ocean Modelling 12 (2006) 32–45
35
full domain size excluding the sponge layer is 1540 km · 624 km, the grid size is 4 km · 4 km, and the vertical sigma grid has 51 levels. Forcing at the three open boundaries is provided by the M2 barotropic tidal elevation predicted by the global tidal model of Le Provost et al. (1998). The same initial horizontally uniform stratification as Holloway (1996) is used, that is, the average of February, March, and April salinity and temperature data taken from the Levitus atlas (Levitus, 1982). A relaxation layer exists at the outer 10 grids of the across-slope boundaries and the offshore boundaries. The model solves the barotropic (vertically averaged) and baroclinic flow components seperately. In the relaxation layer, the barotropic velocity is damped by increased horizontal viscosity of 500 m2/s, and the baroclinic velocity is relaxed gradually to zero by Martinsen and Engedahl (1987)Õs method, and the temperature and salinity are relaxed similarly to values at the outermost grids calculated by the upwind differencing scheme. The horizontal viscosity is fixed at zero in the interior. The model width W (Fig. 1) is varied from 120 km to 1540 km (excluding the relaxation layers) by changing the number of grids. The forcing is always by the prescribed tidal elevation at all three open boundaries. The models are run for 12 days and a tidal analysis is performed on the data from the last day. The energy flux carried by the internal tides is calculated and integrated along the outermost grids (but not in the relaxation layer) as Z ZZ 1 T vn p dl dz dt; T 0 where T is the semi-diurnal tidal period, vn is the baroclinic velocity normal to the boundary, p is the baroclinic component of the pressure, l is the along-boundary coordinate, and t is the time. The detail of the calculation is explained in the appendix of Holloway and Merrifield (1999). The integrated energy flux is shown in Fig. 2(a). Note that vn is the component crossing the
Flux
1.6 1.4
Flux / Width
1200
x-flux y-flux
x-flux y-flux
1000 800
1
Flux/Width [W/m]
Integrated flux [GW]
1.2
0.8 0.6 0.4
600 400 200
0.2 0
0 -0.2
-200 0
(a)
500
1000 1500 W [km]
2000
0
(b)
500
1000 1500 2000 W [km]
Fig. 2. (a) Energy flux across the boundaries integrated along the outermost grids as a function of the model width. y-flux is the offshore flux. (b) The integrated energy flux divided by the length of the shelf as a function of the model width.
36
K. Katsumata / Ocean Modelling 12 (2006) 32–45
boundary so that x-flux in Fig. 2 is the flux evaluated on the across-slope boundaries and y-flux is the offshore flux. The along-shelf flux (x-flux) does not bear a simple relationship with the model width, but the across-shelf flux (y-flux) shows an almost linear increase for W > 500 km. For W < 500 km, this linear relation does not hold, which is emphasized in Fig. 2(b), where the flux is normalised by the model width W. When the model is wider, the flux per unit shelf length is approximately 1000 W/m, but when the model region is narrower than approximately 500 km, the outgoing flux is much reduced. This situation is similar to the underestimate by two-dimensional models mentioned in the introduction, although the topography and the barotropic tidal forcing are not strictly two-dimensional for small W. For example, if the 120-km wide model is used to simulate the internal tides at the middle cross-section of the shelf, the model would predict the outgoing offshore flux of 130 W/m, which is much smaller than the offshore flux predicted by wider (W > 500 km) models (approximately 1000 W/m). It is possible that the centre of the model domain just so happens to coincide with a region of weak generation due to local details of the barotropic forcing and bathymetry. This possibility will be examined in Section 3, but for now it is assumed that the narrow-width limit represents the two-dimensional case. It is conjectured that a strong along-shelf barotropic tidal flow, combined with along-shelf topographic variation, produces additional forcing of the internal tide (Legg, 2004a,b). Baines (1973) defines the internal tide forcing as og o u ov z þ ; ð1Þ F ¼ N2 ot ox oy where N is the buoyancy frequency, g is the surface displacement, and the overbar denotes vertical averaging. This term appears in the density conservation equation after separating the velocity into the barotropic tidal field of an equivalent unstratified ocean and the internal tidal field (Baines, 1973, also see Eq. (10)). In the two-dimensional model, the ou=ox term is absent. An interpretation of the forcing can be given in the following way. The potential energy change E of a fluid parcel of unit volume by a barotropic tidal flow in a stratified ocean is approximately dq0 2 ð2Þ f; E¼g dz where f is the vertical displacement of the fluid parcel caused by the barotropic flow, and q0(z) is the stratification of the undisturbed ocean. The rate of work to cause this energy increase is dE dq0 df ¼g ð3Þ 2f ¼ 2fq0 F ; dt dt dz where F is the forcing defined in (1). Using this interpretation, it is possible to provide a measure of the internal tide amplitude in the model. By integrating (3) in space and time and assuming that all the perturbation energy is radiated as the internal tide, (i.e. dE dx dy dz/dt = Flux dl dz) one can define RR Flux dl dz RRR ; ð4Þ f0 ¼ 2q0 Forcing dx dy dz where all terms are time averaged. Alternatively, f0 can be regarded as the ‘‘efficiency’’ of the model as an internal tide generator, for it yields the output flux when multiplied by the input forcing.
K. Katsumata / Ocean Modelling 12 (2006) 32–45 ζ0
Integrated forcing
0.9
37
1.2
Amplitude [m]
Forcing [x106 m4s-3]
0.9 0.6
0.6
0.3 0.3
0
0 0
(a)
500
1000 1500 W [km]
2000
0
(b)
500
1000 1500 W [km]
2000
Fig. 3. (a) Forcing (1) integrated in the model domain as a function of the model width. (b) ‘‘Efficiency’’ or the total outgoing flux divided by the total forcing (see (4)) as a function of the model width.
Fig. 3 shows the forcing integrated over the model region and the corresponding f0 as functions of the model width. It is seen in Fig. 3(a) that the forcing increases almost linearly with the model width W, and as a result the average tidal amplitude, or the efficiency, f0, increases for W < 800 km but saturates for W > 800 km. From Fig. 3(b), it is implied that the narrow-width models are not as ÔefficientÕ as the wide models in generating the internal tide. This result suggests that the reason for the underestimate in the narrow-width (and two-dimensional) models is twofold; one is the reduced forcing and the other is the reduced efficiency. A closer look at Figs. 2(b) and 3(b) suggests that there are two different dynamical regimes in the transition between the two-dimensional and three-dimensional models. One is between W = 0 km and W 500 km, where the normalised flux increases with the model width, and the other is W > 500 km, where the normalised flux is approximately constant. In the large scale regime, the difference in the energy flux is mostly explained by the difference in the forcing because f0 is almost constant and the integrated flux (Flux) is almost proportional to the integrated forcing (Forcing). In the small scale regime, however, not only the normalised flux, but also f0 increases with the model width. This small scale regime is further investigated in the next two sections.
3. Numerical experiments with idealised bathymetry In the experiments described in the previous section, when the model width is narrower than approximately 500 km, the resulting flux is not proportional to the forcing. In order to ascertain that this phenomenon is not caused by the irregular topography and the tidal flow of the particular Australian North West Shelf simulation, an idealised situation is considered in this section.
38
K. Katsumata / Ocean Modelling 12 (2006) 32–45
An idealised shelf topography y ð5Þ H ðyÞ ¼ D þ A tanh b with D = 128 m, A = 16 m, b = 2 km is used. The model grid size is fixed at 0.8 km · 0.8 km, and the across-slope size of the model domain is 51.2 km (64 grids) including 8 km (10 grids) wide relaxation layers on the across-slope boundaries and the on/offshore boundaries. The relaxation scheme is identical to the one used in the realistic model in Section 2. Tidal forcing is given by the surface elevation with an amplitude of 1 cm at the M2 tidal frequency x and with opposite phases at the two long-shelf boundaries. Fig. 4 shows the forcing, which is almost perfectly along-slope and confined near the slope. Constant stratification of N = 1.07 · 102 s1 is set up with the salinity while the temperature is fixed constant. The Coriolis parameter is fixed at the value of 30 north. The model was run for four days from quiescent and horizontally uniform initial conditions and results from the tidal analysis from the last day are examined. The same experiment as in Section 2 was carried out and the result is shown in Fig. 5. The horizontal length scale is nondimensionalised by the internal tide characteristic wave length L0 = 2D/c (approximately 23 km here), and the forcing and the flux are nondimensionalised by ð6Þ F 0 ¼ N 2 VD=L0 and G0 ¼ q0 A2 V 2 N 2 H =ð2xL0 Þ; ð7Þ 1 respectively, where V is the typical across shelf barotropic flow speed (here 3.5 cm s ), and c is the characteristic slope (here 1.1 · 102) of the internal tide defined by x2 f 2 c2 ¼ 2 . ð8Þ N x2
y
1
1 5 5 1
0
1.0
-1 -1
0 x
1
Fig. 4. Forcing field in the simple numerical model when the model width is 51.2 km. The ellipses show the nondimensionalised barotropic forcing flow. The scale of ellipses is shown on the bottom right. The contours show the nondimensionalised strength of the forcing defined in (1) with the contour interval of 2 unit forcing. The size of the ellipses and the forcing are nondimensionalised by V and DF0, respectively (see text for definition). The shaded area is the relaxation layers.
K. Katsumata / Ocean Modelling 12 (2006) 32–45 1
39
0.8
0.8 0.6
× G0
× DF0
0.6 0.4
0.4 0.2 0.2
0
0 0
(a)
5 10 15 Model width (× L0)
20
0
(b)
5 10 15 Model width (× L0)
20
Fig. 5. (a) Across-slope offshore energy flux as a function of the along-slope model width. The flux is integrated over the on- and off-shore boundaries and normalised per unit length of the shelf. (b) Forcing as a function of the alongslope model width. The forcing is integrated in depth and averaged in area except the relaxation layers.
The forcing scale (6) is taken from (1), and the flux scale (7) is based on the dimensional form of the analytic solution (17) discussed in the following section. It is demonstrated that even though the forcing per unit area is almost constant, the resulting flux drastically decreases when the model width is less than approximately 5 internal tide wave lengths. Fig. 6 is an example of the model output, showing the vertically integrated and tidally averaged energy flux when the model width is 51.2 km. The flux field shows strong along-slope variability,
y
1
1.1 1 0.9
0
1.0
-1 -1
0
1
x Fig. 6. Vertically integrated and tidally averaged energy flux when the model width is 51.2 km. The shaded zone indicates the relaxation layers. The contour is the depth with contour interval of 0.1. All the quantities are nondimensional.
40
K. Katsumata / Ocean Modelling 12 (2006) 32–45
not seen in the forcing distribution shown in Fig. 4. The fact that the offshore flux is the strongest near the middle (x 0) suggests that this along-slope variability comes from the across-shelf relaxation layers. Having tracked down the problem to a simple setting, we can now use an even simpler semianalytical model to understand the reason for the underestimate in the narrow-width models.
4. A simple model of three-dimensional internal tide generation 4.1. Formulation Governing equations of an inviscid and linear internal wave field with the Boussinesq and hydrostatic approximations are p ot u fv ¼ ox ; 0 q p ot v þ fu ¼ oy ; 0 q p qg ð9Þ 0 ¼ oz ; 0 q 0 q dq dq ot q þ w 0 ¼ W 0 ; dz dz ox u þ oy v þ oz w ¼ 0; where (u, v, w) are the velocities associated with the internal waves, p is the pressure perturbation 0 is a constant reference density, q0(z) is the unperturbed density, and q is caused by the waves, q the density perturbation. W is the vertical velocity associated with the barotropic tide defined by og o u ov z þ W ¼ ð10Þ ot ox oy (cf. Eq. (1)). For simplicity, we assume the buoyancy frequency N is constant. The boundary conditions in the vertical direction are w ¼ 0 at z ¼ 0; w ¼ uox H voy H
at z ¼ H.
ð11Þ
The other boundary condition is the radiation condition for energy at open horizontal boundaries. These equations and the boundary conditions are identical to those used by Baines (1973). The topography H(x, y) is prescribed as (5) and the barotropic tide ðuðx; y; tÞ; vðx; y; tÞÞ is given by Q ð12Þ u ¼ 0; v ¼ cos xt; H where Q is a constant. With the rigid-lid approximation, (10) becomes W ¼ AQ
z 1 y sech2 cos xt; b H2 b
K. Katsumata / Ocean Modelling 12 (2006) 32–45
for the bottom topography (5). In the narrow-shelf limit, b ! 0 and the forcing becomes z W ¼ 2AQ dðyÞ cos xt. ðD þ A sgn ðyÞÞ2
41
ð13Þ
Offshore (y > 0), the bottom is well approximated by D1 = D + A = constant in the narrow-shelf limit and the solution can be obtained using a vertical mode expansion. The solution for w, subject to the boundary condition (11) is given by 1 X npz ixt wn ðx; yÞ sin e . ð14Þ wðx; y; zÞ ¼ D1 n¼1 Henceforth, the physical quantities are nondimensionalised by D1 (vertical length), L0 = 2D1/c (horizontal length), and 1/x (time) unless otherwise noted. It is also assumed that all quantities have time-dependence eit. Then the formal solution for the vertical nth mode wn is Z 1 Z 1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð1Þ in2 p2 H 0 2np ðx nÞ2 þ ðy gÞ2 W n ðn; gÞ dn dg; ð15Þ wn ðx; yÞ ¼ W n 1
1
ð1Þ H0
is the Hankel function of the first kind and Wn is the vertical nth mode of the barowhere ð1Þ tropic forcing W given by (10). The radiation condition of the energy is satisfied by taking H 0 ð2Þ rather than H 0 . Essentially the same solution has been given by Llewellyn Smith and Young (2002). The physical interpretation of the solution (15) is that the internal tide response at superposition of (x, y) is the response to the local forcing, Wn, plusqthe the response from the ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð1Þ
forcing at (n, g) with the GreenÕs function H 0
2np ðx nÞ2 þ ðy gÞ2 .
Since only the offshore region y > 0 is considered, the integration in g should be limited to g > , where is an infinitesimal positive number. The experiment in Section 3 is replicated by limiting the along-shelf extent of the source within [S, S]. Thus, Z S qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð1Þ n aH 0 2np ðx nÞ2 þ y 2 dn; wn þ W n ¼ 4iðÞ np ð16Þ S
where a is the nondimensional forcing strength a ¼ DA1 D1QLx . From (16) and (9), it is possible, by numerically integrating (16), to calculate the energy flux, 2 2 2 iG oG f G oG ! Z 0 u 1 p N A Q ox oy x R ; dz ¼ R f oG oG 2 q x D D Lx v 1 1 iG oy þ x G ox p 0
ð17Þ
where R; ð Þq means the real part and the complex conjugate, respectively, and G ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R S ð1Þ H 0 2np ðx nÞ2 þ y 2 dn. Fig. 7 shows an example with the same source width as S
Fig. 6. Because of the idealised d-function forcing, the magnitude of the energy flux (17) does not explicitly depend on the vertical mode n, but it is seen that the higher modes have smaller spatial scales. Gravest 10 modes are chosen to plot Fig. 7 to compare the analytical solution with the numerical simulation with expectation that higher modes are prone to the stronger damping. This choice is based on a suggestion by St. Laurent and Garrett (2002) that internal waves with equivalent modes less than 10 are likely to radiate away. Despite the many simplifying assumptions, a
42
K. Katsumata / Ocean Modelling 12 (2006) 32–45
y
1.0
0.5
-1
0
x
1
1.0
Fig. 7. Offshore flux calculated by (17) for the solution (16). The source width is S = 1.0, which corresponds to Fig. 6.
qualitative similarity can be found between Figs. 6 and 7 such as the magnitude of the flux and weak radial refraction. The stronger damping in the numerical solution might be caused by incomplete absorption in the relaxation layers. Two obvious limits exist for (16). First, the point source limit can be obtained by replacing the pffiffiffiffiffiffiffiffiffiffiffiffiffiffi constant a with d(n) and it can be seen that all the quantities are now functions of r ¼ x2 þ y 2 only. It is also possible to show that the outgoing energy flux is constant when integrated across the surface r = constant, which means the energy flux per unit length (energy flux density) decreases as r1. Second, the two-dimensional limit is obtained by S ! 1, where it is possible to analytically integrate (16) such that Z S qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð1Þ n aH 0 2np ðx nÞ2 þ y 2 dn; ¼ 4iðÞn nei2npy ð18Þ wn þ W n ¼ 4iðÞ np S
(see e.g. Moriguchi et al., 1987). This limit can also be obtained by neglecting oxo terms in (9). In fact, these two limits give an explanation of the result of the numerical experiments in Section 3. Assuming that the relaxation layers work perfectly and no spurious reflection occurs, the narrow width limit in the numerical experiments corresponds to the point source limit above. In this limit, the energy flux per unit length of the offshore boundary decreases as r1. On the other hand, in the two-dimensional limit, the energy flux per unit length of the offshore boundary is constant for all y > 0. The dependence of the flux on the source width S is evaluated and shown in Fig. 8. The flux increases as S increases and saturates around 2S 5, in reasonable agreement with the result of the numerical experiments shown in Fig. 5. Thus it is shown that the underestimation in narrow-width models is a result of the artificial across-shelf relaxation (or radiation) boundaries.
5. Discussion An obvious remedy to avoid the effect of the across-shelf artificial boundaries is to use periodic boundary conditions in the along-shelf direction. In fact, when this boundary condition is used, the flux estimates by numerical simulations and theoretical estimates agree well (Khatiwala, 2003).
K. Katsumata / Ocean Modelling 12 (2006) 32–45
43
Offshore flux (× F0)
0.8
0.6
0.4
0.2
0 0
5
10 15 Source length 2S (× L0)
20
Fig. 8. Offshore flux integrated along y = 1.0 and normalised per unit length as a function of the source width S.
In realistic applications, however, the use of the periodic conditions is complicated by two factors. First, in most applications both the tidal forcing and the topography show significant along-slope variation (for example Holloway, 2001; Cummins and Oey, 1997), and oversimplification of such variation may result in unreliable estimates. Second, in order to accommodate the periodic boundary condition, the forcing must be in phase along the shelf, producing unrealistic barotropic tidal flow and energy flux. The experiment in Section 2 was repeated with the narrow model width of 52 km, whose forcing and topography are taken near the middle of the model region (x = W/2 in Fig. 1), but with the periodic boundary condition in the along-shelf direction. The resultant offshore flux was unrealistically small (less than 1% of the wide scale models in Fig. 2). Therefore for a realistic estimate, a fully-three-dimensional numerical simulation seems to be the appropriate method to use. In Figs. 5 and 8, it was found that the saturation of the normalised offshore energy flux occurs when the model width is approximately 5 internal wave lengths. The result shown in Fig. 2 can be interpreted similarly if an appropriate depth of the shelf is chosen. It is possible to calculate the forcing (1) from the stratification and the tidal data for the region of interest without actually
Fig. 9. Vertically integrated forcing (1) on Australian North West Shelf. The depth contours are shown every 100 m up to 1000 m depth, then every 1000 m.
44
K. Katsumata / Ocean Modelling 12 (2006) 32–45
running a numerical model. Fig. 9 shows the forcing thus calculated for the experiment in Section 2 with the full domain size. It is seen that the maximum forcing is localised along the steep shelf slope around 1000 m depth and weaker peaks near the shelf break at 200 m depth. Below the thermocline, the buoyancy frequency is approximately N = 7 · 103 s1, and the corresponding characteristic slope for the M2 tide is 2 · 102. Therefore at 1000 m depth, the internal tide wave length as defined by the distance between two surface reflections of an internal tide ray is approximately 100 km. In Fig. 2, the saturation occurs approximately at W = 500 km, which is about 5 wave lengths. Additional tests were also performed with the simple numerical model in Section 3 with along-slope sinusoidal variation of the topography such that ! y sinð2px=LÞ ; ð19Þ H ðx; yÞ ¼ D þ A tanh þ b cosh2 ðy=bÞ in place of (5). The along-slope model widths were changed as in Section 3. The results for different along-slope variations (L = 1.0L0, 1.2L0, 1.8L0, 2.3L0) were all similar to Fig. 5, that is, when the model width is 5L0, the offshore flux per unit length is more than 90% of the corresponding wide (W = 18L0) models. It is thus confirmed that to obtain realistic tidal simulation and energy fluxes, the model domain must be greater than 5 internal wave lengths.
Acknowledgments The author thanks Drs. G. Symonds, P. Craig, and M. Merrifield for helpful discussions. This work was initiated from a chat with late Dr. Peter Holloway. Comments from two anonymous reviewers helped to improve the manuscript and are gratefully acknowledged.
References Baines, P.G., 1973. The generation of internal tides by flat-bump topography. Deep-Sea Res. 20, 179–205. Baines, P.G., 1974. The generation of internal tides over steep continental slopes. Philos. Trans. R. Soc. London A 277, 27–58. Baines, P.G., 1982. On internal tide generation models. Deep-Sea Res. 29, 307–338. Blumberg, A.F., Mellor, G.L., 1987. A description of a three-dimensional coastal ocean model. In: Heaps, N. (Ed.), Three-dimensional Coastal Ocean Models, Coastal and Estuarine Science. American Geophysical Union, AGU, Washington, DC, pp. 1–16. Craig, P.D., 1987. Solutions for internal tidal generation over coastal topography. J. Mar. Res. 45, 83–105. Craig, P.D., 1988. A numerical model study of internal tides on the Australian Northwest Shelf. J. Mar. Res. 46, 59–76. Cummins, P.F., Oey, L.-Y., 1997. Simulation of barotropic and baroclinic tides off northern British Columbia. J. Phys. Oceanogr. 27, 762–781. Holloway, P.E., 1996. A numerical model of internal tides with application to the Australian North West Shelf. J. Phys. Oceanogr. 26, 21–37. Holloway, P.E., 2001. A regional model of the semidiurnal internal tide on the Australian North West Shelf. J. Geophys. Res. 106, 19625–19638. Holloway, P.E., Merrifield, M.A., 1999. Internal tide generation by seamounts, ridges, and islands. J. Geophys. Res. 104, 25937–25951.
K. Katsumata / Ocean Modelling 12 (2006) 32–45
45
Khatiwala, S., 2003. Generation of internal tides in an ocean of finite depth: Analytical and numerical calculations. Deep-Sea Res. 50, 3–21. Legg, S., 2004a. Internal tides generated on a corrugated continental slope. Part I: Cross-slope barotropic forcing. J. Phys. Oceanogr. 34, 156–173. Legg, S., 2004b. Internal tides generated on a corrugated continental slope. Part II: Along-slope barotropic forcing. J. Phys. Oceanogr. 34, 1824–1838. Le Provost, C., Lyard, F., Molines, J.M., Genco, M.L., Rabilloud, F., 1998. A hydrodynamic ocean tide model improved by assimilating a satellite altimeter-derived data set. J. Geophys. Res. 103, 5513–5529. Levitus, S., 1982. Climatological atlas of the world ocean. NOAA Prof. Rep. 13. Environ. Res. Lab., Geophys. Fluid Dyn. Lab., US Department of Commer., Princeton, NJ, 173pp. Llewellyn Smith, S.G., Young, W.R., 2002. Conversion of the barotropic tide. J. Phys. Oceanogr. 32, 1554–1566. Maas, L.R.M., Lam, F.-P.A., 1995. Geometric focusing of internal waves. J. Fluid Mech. 300, 1–41. Martinsen, E.A., Engedahl, H., 1987. Implementation and testing of a lateral boundary scheme as an open boundary condition in a barotropic ocean model. Coast. Eng. 11, 603–627. Moriguchi, S., Udagawa, K., Hitotsumatsu, S., 1987. In: Sugaku Koshiki II (mathematical Formulae II). Iwanami Shoten, Tokyo, p. 340 (in Japanese). Munk, W., 1966. Abyssal recipes. Deep-Sea Res. 13, 707–730. Munk, W., Wunsch, C., 1998. Abyssal recipes II: Energetics of tidal and wind mixing. Deep-Sea Res. 45, 1977–2010. New, A.L., 1988. Internal tidal mixing in the Bay of Biscay. Deep-Sea Res. 35, 691–709. Pingree, R.D., New, A.L., 1989. Downward propagation of internal tidal energy into the Bay of Biscay. Deep-Sea Res. 36, 735–758. Pingree, R.D., New, A.L., 1991. Abyssal penetration and bottom reflection of internal tidal energy in the Bay of Biscay. J. Phys. Oceanogr. 21, 28–39. Rattray Jr., M, 1960. On the coastal generation of internal tides. Tellus 12, 54–62. Sandstrom, H., 1976. On topographic generation and coupling of internal waves. Geophys. Fluid Dyn. 7, 231–270. St. Laurent, L., Garrett, C., 2002. The role of internal tides in mixing the deep ocean. J. Phys. Oceanogr. 32, 2882–2899.