Setting tidal forcing for regional modelling of internal waves

Two methods for setting tides in regional numerical models focussed on replication of internal waves are discussed in this paper. A traditional procedure prescribes the tidal forcing at the boundaries of the model domain. As it follows from our analysis, this method works well in some cases but does not allow to control identical generation conditions in all points of the model grid. Tidal ellipses are spatially dependent on the relative distance to the model boundaries. An alternative method for activation of tidal motions in general circulation models is to set a tidal body force on the right-hand side of the momentum balance equations. The procedure of how to do this correctly is in the focus of the present study. It is based on the analytical solution of the shallow water equations and includes compensating terms providing blockage of the generation of spurious inertial oscillations. The method was tested in a wide range of input parameters and showed a good match of the model-predicted tidal currents in all points of the model domain. The method was verified against an analytical solution by comparing model-predicted and theoretical internal wave amplitudes in case of realistic pycnocline-type oceanic stratification.


Introduction
The generation of internal waves by tides that streamline bottom topography in the stratified ocean is of great interest to the scientific community.According to Munk and Wunsch (1998), the ocean without tides would turn into a stagnant pool of salty water.The numerous measurements of internal waves in the ocean demand proper three dimensional modelling that can be used as a helping tool complementing sparse observations.The question about setting tides in the models is very important for the correct interpretation of in-situ collected data.
At present there two approaches for setting tides in the models: (i) from the boundaries and (ii) as a forcing in momentum balance equations.In the first approach, three types of boundary conditions are applied, i.e. fixed elevation, fixed velocity, or combination of both with an account of radiation conditions, (A., 1976).Their performance was investigated by Carter and Merrifield (2007) using the Princeton Ocean Model (POM).To avoid reflection of generating internal waves from the model boundaries, they used a group of areas with high viscosity, often referred as sponge layers.Munroe and Lamb (2005), Ma and Mandes (2011) and Rainville et al. (2010) also initiated tidal forcing from boundaries in modelling with POM; similar method was applied by Klymak et al. (2016) in experiments with the Massachusetts Institute of Technology general circulation model (MITgcm, Marshall et al., 1997).
The second method in setting tides is to include a body force into the momentum balance equations.S. (2003) using the MITgcm as a tool for modelling of internal tides added tidal body force to the zonal momentum balance equation.The model output was validated against an analytical solution obtained in a two-dimensional case for monotonously stratified fluid.The comparison was focused on the estimation of integral wave characteristics such as total energy fluxes and tidal energy conversion rate without analysis of the parameters of every particular baroclinic mode.Later, Lorenzo et al. (2006) considering similar problem of internal waves generated by a zonal tidal flow, included also a Coriolis parameter dependent body force term in the second momentum balance equations which allowed to suppress generation of spurious inertial oscillations.The analytical expression of the rate at which energy is converted from barotropic tide to internal waves was provided by Pétrélis and Llewellyn Smith (2006).
In this paper, we outline a new method for setting the tidal forcing in regional numerical models.It is based on introducing two periodic terms in each of the momentum balance equations to control the parameters of tidal oscillations and suppress the generation of spurious inertial oscillations in all possible oceanic conditions.The method is tested by comparing the amplitudes of tidally generated internal waves with the analytical solution.
The paper is organized as follows.A general form of tidal ellipses is presented in Section 2. Section 3 focuses on two methods of setting the tidal forcing in a general three-dimensional case.The water depth is assumed variable, although no baroclinic effect is considered at this stage (barotropic tidal forcing).The comparative analysis of both settings, i.e. from boundaries or using the tidal body force, is the key element here.In Section 4, we present an analytical solution for the internal modes generated by tides in the continuously stratified ocean.This section includes also the validation of the body force method by considering different oceanic scenarios.Some wider applications of the proposed methodology are discussed in Section 5.

Tidal ellipses
Theory predicts elliptical trajectories of water particles produced by tides in the rotating ocean.The ellipse parameters, i.e. the length of major and minor axis, and the inclination angle are the functions of the tidal flow velocity components.To find this relationship, consider the tidal flow in the Cartesian coordinate system (0) with -axis and -axis directed to the East and to the North, respectively.For each of tidal harmonics the eastward  and northward  tidal velocities are the periodical functions: (1) Here  0 and  0 are amplitudes,  is the tidal frequency,  is the time,   and   are  and  the tidal phases.Superposition of these two velocity components in space and time produce tidal ellipses as that shown in Fig. 1.The major axis  and minor axis  are expressed in terms of the tidal parameters as follows (Stashchuk et al., 2017): The inclination angle  (with respect to the eastward direction) is (3) Here the phase lag  =   −  .Taking into account that the tidal ellipse parameters (2)-(3) depend only on the phase lag , formula (1) can be simplified as follows: The tidal currents in numerical models can be introduced using (4) by setting the amplitudes  0 ,  0 of the eastward (zonal) and northward (meridional) tidal currents and the phase lag  between them that control the axes of the tidal ellipse ,  and its inclination angle .
Having the relationship between tidal currents (4) and the parameters of tidal ellipses (2)-(3), the next question is how to introduce tidal currents into the numerical models technically.There are two commonly used methods of barotropic tides initiation.One of them assumes setting the periodic flow conditions (4) at the outer boundaries of the model domain.The adjustable parameters  0 ,  0 , and  can be tuned in such a way to achieve the best fit of the model-produced tidal flow with observations in all points inside the model domain.However, this is not always an easy task having in mind various interference conditions of two tidal components in different parts of the considered area and possible reflection from the boundaries.
An alternative method of tidal motion activation is to set periodical forcing terms directly in the momentum balance equations.Theoretically, this method should provide similar generation conditions in all grid points.The great advantage of such setting is it does not produce spurious boundary artefacts.Note, although, that the application of this method can be limited in shallow water areas with complex topography and islands.Below we consider both methods in detail, focusing on their performance, accuracy in the reproduction of tidal flow and possible limitations.Methodology wise, the comparative analysis can be made using any oceanic general circulation model.Our choice is the Massachusetts Institute of Technology general circulation model, MITgcm, (Marshall et al., 1997) which has a number of external users' packages and allows free operation with both model settings.

Setting tides in the model
In this section, we compare two methods of setting tides into the numerical model without considering any baroclinic effects, specifically, excluding the generation of internal waves from the analysis.In doing so, we assume the fluid is homogeneous in this series of experiments.
Consider a rotating basin in the Cartesian coordinate system (0) with the -axis and -axis directed to the East and to the North, respectively, and the -axis is vertically upward ( = 0 is the undisturbed free surface).The localized bottom topography (, ) in these experiments has a form of a Gaussian-type underwater bank: where  0 =200 m is the reference water depth,  max =70 m is the bank's height, and  = 20 km is the topography horizontal scale.Methodology wise, these numerical experiments are aimed to study the establishment of tidal flows around and over the bank during ten tidal periods and to quantify the difference between two methods of tidal forcing.
The central part of the model grid was 1414 × 1414 nodes in  and  directions with 100 m grid step.In vertical direction the grid has 41 layers with 5 m resolution.In case of the stratified ocean when the generation of internal waves is considered, additional 93 grid points with a telescopic increase of the horizontal resolution from 100 m in the central part to 3.5 km at the periphery were added to the lateral boundaries.Such a grid extension excludes backward reflection of internal waves during ten tidal cycles, the time required internal waves to reach the boundaries.
The M 2 tidal harmonic is considered in this study for setting the frequency , although the method is applicable to any other tidal constituent.The velocities  0 ,  0 in (4) are assumed to be spatially variable (the functions of  and ), and no restriction on the latitude is applied.The only condition is the requirement  < , where  is the Coriolis parameter,  = 2 sin ,  is the angular rotation rate of the Earth, and  is the latitude.
An analytical solution for topographic generation of internal waves in the continuously stratified ocean is used in this study for the method validation.This solution is linear, therefore in the numerical experiments we keep all the parameters in the range of linearity assumption.A quantitative measure of the contribution of the non-linear effects into the process of tide-topography interaction reads:  = tidal excursion amplitude topographic length scale × topographic height fluid depth .
It was shown by Gerkema and Zimmerman (1995) that the problem of tide-topography interaction in a stratified fluid can be regarded as linear when this parameter is small, i.e.  ≪ 1.In our further numerical analysis, we will keep all the parameters to be in this range.

Tidal forcing from boundaries
The simplest way to initiate tides in the models is to prescribe periodical tidal motions on boundaries.In the present series of experiments, we set the tidal velocities (4) at the southern and western boundaries, Fig. 2 (the MITgcm solves hydrodynamics equations using the finite volume technique and requires to specify prescribed horizontal velocities at the boundaries without setting sea surface displacements).Note, however, in case two other boundaries are kept closed, the tidal energy unable to escape from the model domain which results in the generation of standing waves.To cope with the wave reflection from the model boundaries, we use Orlanski boundary conditions (I., 1976) at the northern and eastern boundaries.They are based on the equation of a free propagation of any quantity  through the boundary: where  is the outward normal, and  is the wave phase speed.In case of the shallow water systems A series of numerical experiments has been conducted with different forcing parameters  0 ,  0 , and .We start our analysis with the simplest case when the phase lag between eastward and northward velocities is equal to zero.With =0 the tidal ellipse turns into a straight line presented by only major axis  (see formula (3)).The orientation of the line is controlled by the ratio of the velocity amplitudes  0 and  0 .For instance, by setting  0 =  0 = 0.005 m s −1 the expected direction of particle propagation should be 45 • with respect to the -axis.
The experiment was designed to examine the influence of the Earth rotation (latitude) on the resulting tidal ellipses produced by the model.They are presented in Fig. 3 for three latitudes,  = 0 • , 30 • , and 60 • (hereafter the blue colour is used for original ellipses defined by formula (4), and the red colour illustrates the model output).For better presentation all ellipses were scaled by a factor 1000.The principal conclusion that can be drawn from Fig. 3 is that with zero phase lag,  = 0, the model produces tidal ellipses with the minor axis anyway instead of generation straight lines.Their shape and orientation strongly depend on the rotation.The model-produced ellipses are closer to straight lines in case without rotation (left panel).However, with the increase of rotation, the orientation of the ellipses dramatically changes compared to the original direction.Consider now the model results obtained for latitude  = 30 • when the tidal lag  ≠ 0 ( ≠ 0, formula (3)).Fig. 4 presents original (blue) and the model obtained (red) tidal ellipses for four values of the phase lag, i.e.  = 0; ∕6; ∕3; ∕2.The amplitudes of tidal velocities  0 and  0 in these experiments were taken 0.005 m s −1 .It is clear from Fig. 4 that the increase of the phase lag between the velocities at the boundaries results in a better match between red and blue ellipses.They are almost identical when  = ∕2, although the coincidence is still not perfect.
The considered examples reveal the fact that in a rotating ocean, the model-generated tidal velocities within the model domain are different from those prescribed at the boundaries.The following model run explains the reason for this difference.In this experiment the model is forced only from the western boundary with the velocity amplitude  0 = 0.005 m s −1 leaving the southern boundary unforced, i.e.  0 = 0. Quite a moderate rotation was set in this experiment (latitude  =30 • ).Fig. 5 shows time series of eastward (in blue colour) and northward (in red colour) velocities recorded at points 1,2,3 and 4 depicted in Fig. 2. As one can see from Fig. 5, being forced from just one western boundary, the model produces also north-south oscillations (and thus  velocity) due to the rotation.Important is that these newly produced northward velocities are out of phase to the eastward ones at all considered points (phase shift is close to ).In addition, the model-produced -velocity is also different from those prescribed at the western boundary.Its amplitude within the model domain exceed the initial value 0.005m s 1 in all four points.
Thus the tidal forcing from only western boundary generates also northward velocities  which, in turn, have a positive feedback to the -velocity.The generation of extra -velocities by the model explains the decrease of ellipses inclination with the increase of the rotation, Fig. 3.By setting the tidal forcing at the model boundaries, one must be prepared to the generation of some extra spurious currents with spatially variable phase leg that can lead to the spatial variability of the tidal signal within the model domain.

Tidal forcing as a body force in momentum balance equations
An alternative approach for setting tides in the ocean numerical models is to introduce a tidal body force on the right-hand side (RHS) of the momentum balance equations.Below we present the idea of how the expression for body force can be found.
Consider a system of shallow water equations for homogeneous fluid.An assumption of the rigid lid at the free surface is applied.The simplest way to initiate tidal velocities in the momentum balance equations using formula (1) is as follows: Without rotation,  = 0,  0 and  0 are the amplitudes of the  and  velocities that satisfy the shallow water equations.The Coriolis force introduces a lateral acceleration that makes these two equations coupled.To clarify the role of the rotation some conversion of ( 7) is required.Applying the cross-differentiation procedure system (7) can be rewritten in the form of two equations, one for  velocity only: and another for  velocity: It is seen that Eqs. ( 8) and ( 9) contain two forcing terms in their RHS.Both terms are periodic functions with the frequency  but have different amplitudes and the phase shifts   and   .Applying simple trigonometrical procedures to the RHS of Eqs. ( 8) and ( 9), they can be converted into the following expressions  2 ũ0 sin( − φ ) and  2 ṽ0 sin( − φ ), respectively, where ũ0 , φ , ũ0 , and φ are as follows: Then the solutions of differential Eqs. ( 8) and ( 9) are: It is clear that velocities (11) are quite different from the expected form (1).They coincide at  = 0, but the rotation introduces substantial distortion into solution (11).As a result, by setting the tidal forcing according to (7), it is difficult to control the parameters of the resulting tidal ellipses.Moreover, the velocity value enormously increases with the increase of rotation when  → .
To suppress the influence of the rotation, we propose to add additional terms to the tidal body force in the RHS of the momentum balance equations, as shown below: With these additional terms after the routine mathematical procedure, the system (12) takes the form: The solution of the system of differential Eqs. ( 13) is: This solution is quite simple (compare with ( 10) and ( 11)) in which the resulting amplitudes  0 and  0 of zonal and meridional currents, and their phases   and   can be taken as input model parameters for tidal ellipses (1), from observations, for instance.Note that solution (14) was found in the basin of constant depth.In case of variable depth, we can assume constant water discharge over varying topography which is applicable for the regional scale modelling.Then the tidal forcing   and   in the momentum balance equations for the case of variable depth (, ) can be generalized as follows: Here  0 (, ) =   0 ∕(, ), and  0 (, ) =   0 ∕(, );   0 = const and   0 = const are the water discharges in  and  directions.Some further simplifications can be introduced taking into account that only the phase lag  =   −   (it was shown above) matters for setting the ellipses, not the absolute values of the phases   and   .The examples below show how it works in real model experiments with the MITgcm.Fig. 6 shows the scheme with the arrangement of the model forcing and boundary conditions.The blue colour of the model domain means that  0 (, ) and  0 (, ) tidal velocities are activated in the whole area using ( 15).Technically, an extra subroutine with the defining of tidal forcing can be added to the MITgcm code using mypackage option.With such a setting, all boundaries should be open with zero derivatives for all variables, Fig. 6.
Consider the difference between Figs. 2 and 6.In setting the telescopic extension of the grid near boundaries, in the case with the stratified ocean, one should provide a smooth telescopic grid step increase both for (i) surface and (ii) internal waves.Really, in a stratified ocean of variable depth, both long barotropic and shorter internal tidal waves are generated.Physical properties of these waves, i.e. the vertical structure, wavelength, phase speed, are quite different.Barotropic tidal waves are much longer and faster than generated by tides internal waves.A two-step telescopic grids should provided to avoid reflection of both types of waves from boundaries.The first additional grid a quite gradual increase of steps that works for internal waves.The second grid makes a more rapid change of the steps to provide a substantial increase of a distance that required for fast surface waves.This combination of the grid stretching provides quite a long time of propagation, both surface and internal waves from the generation site to the boundary.Within this time interval, there will be no backward reflection of both types of waves.This methodology works quite well for localized bottom topographies like seamounts or oceanic ridges, and also in case of restricted areas of slope-shelf topography.
The best way to illustrate the performance of an alternative method of the tidal forcing in the numerical model is to conduct similar experiments that were done with the forcing from boundaries and to compare the model results.Fig. 7 presents the velocity time series recorded at the same four control points and parameters ( 0 = 0.005 m s 1 ,  0 = 0,  =30 • ) as it was in Fig. 5.
The principal difference between the time series shown in Figs. 5 and 7 is that with the tidal body force (15) the model perfectly reproduces velocities that coincide with their initial presentations ( 14).With tidal forcing in the momentum balance equations, the model does not generate spurious meridional velocities  that are visible in Fig. 5.Moreover, at all considered points, the resulting zonal velocities  have the amplitudes  0 = 0.005m s −1 .
We repeated all previous experiments using tidal body forcing (15).The model-produced tidal ellipses perfectly coincided with those shown in Figs.3-4 by blue colours.

Generation of internal waves by tides in stratified fluid (analytical and numerical solutions)
Two methods for settings the tidal forcing in ocean models presented above revealed inconsistency of their model outputs.The better performance of the 'tidal body force' procedure in comparison to the 'boundary forcing' method for homogeneous water requires its further testing in case of stratified fluid when the tidal flow interacting with bottom topography generates internal waves.
The best way to validate a numerical model is to compare its outputs against any analytical solution.In the present study, we use the Fourier transforms and complex analysis to derive a solution for internal waves generated by tide over an underwater ridge.A threeparametric buoyancy frequency profile is used in the analysis, which allows a reasonably realistic smoothed oceanic pycnocline.

Analytic solution
Consider a rotating basin of variable depth filled with stratified fluid.The bottom topography () is a ridge extended in the meridional direction (  = 0): Here 2 is the ridge width and  max is its height.Fig. 8 shows the bottom profile along with the buoyancy frequency profile () = (− 0 ∕  ) 1∕2 (discussed below), where  0 () is the undisturbed background density,   is the mean value of density,  is the acceleration due to gravity.A linearized two-dimensional (, ) system of Euler equations in terms of the stream function   =   ,  = −  . ( and vorticity  =   +   reads: +    2 ()  ∕ = 0.
Here (, , ) is the wave-induced density deviation.For the baroclinic response of the ocean to the periodic tidal forcing with frequency , a 'rigid-lid' condition at the free surface is used: At the bottom the slip conditions are applied: The periodic solutions of the problem with frequency  are considered: With these assumptions Eqs. ( 18)-( 20) are reduced to: where Introducing a new variable problem ( 22)-( 23) takes the form where  For the polynom () =  1 ( +  2 ) 2 +  3 ( 1 ,  2 ,  3 are arbitrary constants), the buoyancy frequency () is a three-parametric family of curves The first equation in (25) reduces to the linear Klein-Gordon equation which has an analytical solution expressed in terms of the Bessel function  0 .The coefficients  1 ,  2 and  3 can be expressed in terms of the maximum buoyancy frequency   , the depth of pycnocline   , and its width   , Fig. 8, as follows: In the ocean of constant depth  0 the analytical solution of problem ( 22)-( 23) for stratification (26) has a form of standard orthogonal baroclinic modes (Vlasenko et al., 2005): where  = 1, 2, 3, … is the number of baroclinic mode, and ℎ 0 =  1 (− 0 ).
For finding the amplitudes of internal waves generated by a periodic tidal flow over ridge ( 16), a perturbation method is used.A small parameter  0 is introduced as a relative height of the ridge in (,  1 ) space as follows: Here ℎ max = ∫  max − 0 − 0 ∕().Note, that small height ℎ max of the ridge in (,  1 ) space does not mean its smallness in (, ) coordinates.For instance, for the ridge with height  max = 70 m in a basin of total depth  0 = 200 m the relative height  max ∕ 0 = 0.35 in (, ) coordinates, whereas  0 = 0.0430 in (,  1 ) coordinates.Thus, the transition from (, ) to (,  1 ) space increases the range of applicability of this method to real oceanographic problems.
In the perturbation method the solution (,  1 ) is considered as an asymptotic expansion: Substitution ( 29) into ( 25) leads to a sequence of boundary value problems which are solved by the application of the Fourier transforms and the integral Cauchy theorem (the details are discussed in Vlasenko et al. (2005)).Here we present the final form of the solution: where   = ( + 1∕2),  = 1, 2, 3, … is the mode number and   is the amplitude of −mode: Scrutiny of the solution shows that the wave motion consists of a superposition of the barotropic tidal flow (first terms in ( 30)) and the sum of generated baroclinic modes (the second term in ( 30)).Taking into account that  = −  (formula ( 17)), by differentiating (30) with respect to x, one can find the amplitude of the vertical velocities where of − is the mode number.

Numerical solution
Analytical solution (30) gives the vertical velocity amplitude of every particular internal mode, (32).Note that the model-predicted vertical velocity −field is a superposition of all generated internal modes, i.e.  = ∑ ∞ =1   .For the comparison analysis, it is necessary to separate baroclinic modes in the model output and to compare the model-predicted and analytical amplitudes of every particular mode separately.The method of their separation is based on the orthogonality property of baroclinic modes.The mathematical procedure of the amplitude calculation is briefly outlined below.
In terms of vertical velocity  the wave equation describing internal waves in three dimensions (, , ) reads: In a two-dimensional case considered here the -direction is ignored and the periodic internal waves radiated from the ridge can be presented as follows: Substitution ( 34) into (33) results in ordinary differential equation for  : For internal waves the rigid-lid condition  = 0 is usually applied at the surface  = 0 and at the bottom  = −.A general solution of ( 33) is a superposition of baroclinic modes: where  ±  is the amplitudes of −mode propagating in both directions from the ridge, and   is the wave number of −th mode.
The orthogonality condition for internal modes reads: Multiplying ( 36) by ( 2 ()− 2 )∕( 2 − 2 )  () and applying integration from the surface to the bottom, with orthogonality condition (37) we have: From this formula, two methods for the definition of the wave amplitude  ±  are evident.Having a spatial distribution of the velocity  at a fixed moment of time  one can estimate the wave amplitude  ±  considering a spatial -transect of the -field from the ridge.Ideally, without dissipation, the spatial profile of the waves radiated from topography should be periodical at any distance from the source of generation having the same amplitude  ±  at all distances.Alternatively, one can consider  time series at a fixed point taking the maximum value of oscillation as an amplitude.

Comparison of analytical and numerical solutions
Validation of the model settings with the tidal body force in the RHS of momentum balance equations was conducted for the following stratification parameters:   = 30 m,   = 50 m and   = 0.02 s −1 , Fig. 8.The bottom topography is the ridge specified by ( 16) with the height  max =70 m in a basin with a total water depth  0 =200 m.The half-width of the ridge  in the experiments varies from 2 km to 65 km, depending on the latitude,  = 0 • ; 30 • ; 60 • .
For the considered model setting the parameter  0 = 0.043 is really small ( 0 ≪ 1), so the analytical solution ( 30) is a good first guess for the model validation.In the numerical model the horizontal resolution  =20 m was small enough for replication of the generated waves and calculation of their amplitudes.The method of the telescopic increase of the grid steps near the model domain boundaries described above was used.In the vertical direction, the grid resolution was  = 1 m.The amplitude of the water discharge  0 in  and  directions was   0 =1 m 2 s −1 ( 0 () =   0 ∕()), and   0 = 0 (thus  0 = 0), and there was no phase lag  in (12),  = 0.
Theory predicts that the generation mechanism of internal waves by a tidal flow over the ridge ( 16) has a resonance nature.Really, the analytical solution for the amplitudes of the generated internal modes (31)-(32) includes the following periodical term: sin(  )∕[(  ) 2 − 2 ].It predicts zero −mode wave amplitude when    = , for  = 2, 3, … .However, the amplitude of th generated mode is maximal for Fig. 9 shows the analytical dependence of the wave amplitude on the ridge width for the first three baroclinic modes (32) (solid lines) and the model-predicted values presented by the dots of the corresponding colour.They were calculated using model time series of the vertical velocity recorded at some distance from the ridge using (37).The dependencies   (),  = 1, 2, 3 are shown for three latitudes: It is seen from the graphs that the amplitude of the first baroclinic mode dominates over all other modes for wide ridge.However, for narrow ridges, the amplitudes of the higher modes are larger.This is the case when the topography is getting near-critical in terms of the bottom inclination (critical conditions ∕ = ()).Note also that the generation efficiency decreases with the increase of rotation (compare maximum amplitudes in all three graphs).
Analysis of Fig. 9 shows also that the numerical model captures all properties of the analytical solution very well.The coincidence of the model-predicted wave amplitudes with the analytical solution is relatively good in the whole range of the bottom width and considered latitudes (the largest discrepancy between the wave amplitudes does not exceed 10%).

Setting the tidal forcing
In the considered above experiments the body force terms   and   in the RHS of the momentum balance equations were used in the form as they stand in (15), i.e. including compensating terms (the balanced forcing case).In the present section we address the question what would be the model output in the 'unbalanced forcing case', i.e. when compensating terms  0 (, ) sin() in   , and  0 (, ) sin(−) in   are ignored.
The difference between the balanced and unbalanced forcing cases is illustrated here for the same stratification parameters as that used in the previous subsection.The width of the ridge for two 'rotational' experiments with  = 30 • and  = 60 • was taken in such a way that only first baroclinic mode is generated, which makes our analysis simpler.These two cases for  =21 km (latitude  = 30 • ) and  =40 km (latitude  = 60 • ) are marked in Fig. 9 by black circles.Apart from the fact that the first mode predominates in the model and analytical solutions, important is that with these widths of the ridge the coincidence between numerically obtained and the analytical amplitudes is just perfect (see Fig. 9).There was no meridional forcing in these experiments, i.e.  0 = 0 and  = 0.For low latitudes the difference between model-produced -velocities in balanced and unbalanced cases can be not significant.For instance, this inconsistency accounts for only 25% at  = 30 • latitude (black and magenta lines in Fig. 10 a).The contribution of spurious inertial oscillations to the model velocities strengthens with the increase of the latitude.As an example, Fig. 11 shows the model time series recorded near the ridge with  =40 km width located at  = 60 • latitude.It is clear from Fig. 11 that unbalanced terms in the  model forcing (15) at high latitude,  = 60 • , lead to 2.5 times increase of the amplitude of vertical velocity (panel a) and enormously high level of the spurious inertial oscillations in the horizontal velocity time series (panels b and c).Quantitatively this result is confirmed by the spectral level calculated from the model output.The energy spectrum level of the inertial frequency can even exceed the energy level at tidal frequency, Fig. 11 d.Figs. 10 and 11 show that in case of unbalanced body force, the numerical model produces both tidal and inertial oscillations.Their superposition in complex form can be presented as follows: Here the amplitudes  1 ,  2 ,  1 , and  2 can be both real or complex values.Introducing the new frequencies  =  −  and   = ( +  )∕2, solution (39) can be rewritten as follows: Formula (40) shows that both velocities,  and  are present in the solution.The high-frequency carrier wave with frequency   is modulated by a low-frequency envelope wave with frequency ∕2.Thus, in case of unbalanced body force, the model does not reproduce precisely tidal ellipses.The velocity vectors in the model domain exhibit variations of amplitude and phase in space and time.The superposition of two waves with slightly different frequencies, results in high-frequency wave with modulation at the lower frequency.It is clearly seen from Fig. 11, when the Coriolis parameter ( = 1.2627 × 10 −4 s −1 ) is close to the M2-tide frequency ( = 1.4053 × 10 −4 s −1 ).This superposition creates a beat envelope around the original sine wave.
In the previous examples, the northward velocity was not activated in the body force ( 0 = 0).To generalize the obtained results, let us investigate how the forcing compensation works in three dimensions.In doing so, in the next series of experiments the parameters  0 and  in (15) are taken as non-zero.The bottom topography ( 16) is still twodimensional extended in the meridional direction (shown as shaded area in Fig. 12 a), although the major axis of tidal ellipses is not

Summary and conclusions
Two methods for setting tidal currents in numerical models aimed at replication of internal tides are considered in this paper.The TPXO model predictions (Egbert and Erofeeva, 2002) can be used in both methods for setting the tidal forcing to describe the internal tide response generated over regional size domains (several hundred kilometres horizontal scale).Tidal currents in the numerical model can be initiated in two ways.A traditional method is to prescribe periodic oscillations at the domain boundaries.This procedure works well in some cases, although it requires controlling free radiation of the generated waves through the model boundaries, which is not always granted.An alternative way is to introduce a periodic body force in the RHS of the momentum balance equations.Such forcing initiates the tidal oscillations providing same generation conditions in all grid points.The comparison of these two methods is in the focus of the present study.
The main conclusion from this paper is that tidal forcing from boundaries is not the best option for achieving the consistency of the wave signal in the whole area.Figs. 3 shows some possible side effects.Specifically, the tidal ellipses are not identical in the central part and at the periphery of the model domain.The rotation introduces another uncontrolled effect.It generates an additional spurious signal that is spatially variable and can completely change the tidal dynamics in the area, Figs.4-5.The principal outcome of our study is that using a tidal body force in the RHS of the momentum balance equations, (15), is free from any side effects and provides same generation conditions in all points of the model domain.We believe that the suggested method works well for regional modelling of internal tidal waves.In restricted domains with horizontal scales of tens or hundreds of kilometres, the generation conditions can be considered as spatially consistent.However, the applicability of this method to the larger oceanic scales must be reconsidered.
Note that in this paper, tidal forcing (15) was checked not only in the barotropic case.It was also verified comparing the model predicted internal wave amplitudes  ±  , (38), with the analytical solution (32) obtained for the internal waves generated in a stratified ocean.A schematic diagram of the generation of the internal tides over a localized bottom topography is shown in Fig. 8.A realistic continuous fluid stratification (26) was used in this analysis, which allowed analytical solution for internal modes.The comparison analysis conducted for a wide range of the bottom topography width and latitudes (from  = 0 • to 60 • ) has shown a consistency of numerical and analytical solutions, Fig. 9.
It is worth mentioning, that tidal forcing (15) assumes a linear problem statement.Generally speaking, the consistency of numerical and analytical solutions in the considered case does not necessarily assume that this method can be extended to all oceanographic conditions.Specifically, the question remains to what extent it is good enough for setting tidal forcing when the Froude number is comparable or larger than 1.A more detailed analysis for finding the answer is required.However, some heuristic speculations can be given considering the time variation of the tidal forcing  , which is assumed to be the water discharge.In the Euler variables the time variations of  can be expressed as follows:

Fig. 1 .
Fig. 1.Tidal ellipse (red) presented as a superposition a clockwise (blue) and counterclockwise (green) tidal velocity components.Angle  shows the ellipse's inclination. 0 and  0 are the tidal velocity amplitudes,  and  are the major and minor axes.

Fig. 2 .
Fig. 2. Schematic diagram of the model domain with the tidal forcing initiated from boundaries.The dotted lines show the area from which the grid resolution telescopically increases to the boundaries.

Fig. 3
Fig. 3 also shows that the model-produced tidal ellipses are sensitive to their spatial position.They are different in different parts of the model domain.For instance, the ellipses are seen as straight lines in the central part of the domain along the diagonal at  = 0 • (left panel).However, in the top left and bottom right corners, the presence of minor axes is quite apparent.The difference between ellipses in the whole area increases with the increase of latitude (see central panel for  = 30 • and right panel for  = 60 • ).Most of the distortion takes place just over the seamount summit where model tidal ellipses have a well developed minor axis.

Fig. 4 .
Fig. 4. The same as in Fig. 3 but for the latitude  = 30 • and different phase lags  depicted in the subplot titles.

Fig. 5 .
Fig. 5. Tidal velocities  (blue) and  (red) recorded at points 1-4, Fig. 2, in case the tidal forcing is activated from boundaries.The latitude  = 30 • .Dotted black lines show the amplitude of the −velocity prescribed at the western boundary.

Fig. 6 .
Fig. 6.Schematic diagram of the model domain with tidal forcing activated in the momentum balance equations.The light blue colour denotes that ,  velocities, and the phase lag  are set equal in all points of the model grid.

Fig. 8 .
Fig. 8.The scheme of the bottom topography with the buoyancy frequency profile.
Fig. 10 is based on the model-predicted time series recorded at a distance   =3 km from the ridge during 20 tidal periods (240 h of the model time) for  = 30 • .The balanced and unbalanced forcing cases are given in black and magenta colours, respectively.Applying orthogonality condition (37) to the time series of vertical velocity , we calculated  + 1 cos( 1   − ) according to formula (38), Fig. 10 a.It is seen that the amplitudes of vertical velocity in the unbalanced case (magenta line) exceed those obtained with the balanced forcing (black line) and the analytically derived amplitude (32).Figs. 10 b and c, present the time series of vertically integrated horizontal velocities: In fact, the horizontal velocities produced by the model include both barotropic and baroclinic signals.The vertical integration removes internal waves from the model output.Thus Figs. 10 b and c illustrate the difference between balanced and unbalanced cases for barotropic tidal motions.The model output without rotational compensation in the body force reveal intermittency in the  ()-velocity signal, Fig. 10 b, with the generation of the  ()-velocity component, Fig. 10 c, which initially was not set in the model (black horizontal line in panel c).It is essential that being unbalanced the tidal forcing leads to a generation of inertial oscillations.The energy spectrum in Fig. 10 d clearly shows their presence.

Fig. 10 .
Fig. 10.Model-predicted time series of vertical baroclinic velocity  (panel a) and vertically integrated barotropic velocities  and  (panels b and c, respectively) sampled 3 km from ridge (18).Magenta colour depicts the numerical solution when compensating terms  0 (, ) () for   , and  0 (, ) sin( − ) are ignored in the body force (15), and the black line shows a fully balanced solution.The width of the ridge  =21 km, the latitude  = 30 • and the forcing tidal parameters were  0 = 0.005 ms −1 ,  0 = 0,  = 0.The ridge height and the stratification parameters are the same as discussed in Section 4.3.(d) Energy spectrum at the sampling point.The meaning of the black and magenta colours is as above.

Fig. 12 .
Fig. 12.(a) Original tidal ellipses set in the model by the body force (15) with parameters  0 =  0 = 0.005 m s −1 and the phase lag  = ∕3.The ridge width scale  = 20.7 km, the latitude  = 30 • and all other model parameters are as that used for Fig. 9.(b) The model predicted vertical velocity of the first internal mode at the cross-ridge transect after ten tidal cycles of the model run.(c) Vertically integrated velocities  and  sampled 3 km from the ridge.(d) The model produced tidal ellipses.