# Open boundary conditions for internal gravity wave modelling using polarization relations

Patrick Marsaleix, Caroline Ulses, Ivane Pairaud, Marine Julie Herrmann, Jochem Willem Floor, Claude Estournel, Francis Auclair

► To cite this version:

Patrick Marsaleix, Caroline Ulses, Ivane Pairaud, Marine Julie Herrmann, Jochem Willem Floor, et al.. Open boundary conditions for internal gravity wave modelling using polarization relations. Ocean Modelling, 2009, 29 (1), pp.27-42. 10.1016/j.ocemod.2009.02.010. hal-02110177

HAL Id: hal-02110177

https://hal.science/hal-02110177

Submitted on 28 Sep 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Open boundary conditions for internal gravity wave modelling using polarization relations

Patrick Marsaleix *, Caroline Ulses, Ivane Pairaud, Marine Julie Herrmann, Jochem Willem Floor, Claude Estournel, Francis Auclair

Laboratoire d'Aérologie – CNRS et Université de Toulouse, 14, Avenue Edouard Belin, 31400 Toulouse, France

This paper proposes an original approach of the open boundary condition problem, within the framework of internal hydrostatic wave theory. These boundary conditions are based on the relations of polarization of internal waves. The method is presented progressively, beginning with a simple case (non-rotating regime, propagation direction normal to the open boundary), ending with a more general situation (rotating regime, multimodal & multi-dimensional propagations and variable background field). In the non-rotating case and as far as we assume that the direction of propagation is locally normal to the open boundary, the so-called PRM (polarization relation method) scheme can be seen as a three-dimensional version of the barotropic Flather boundary conditions. The discrete form of the scheme is detailed. Numerical stability issues proper to leap-frog time stepping are in particular discussed. It is shown that errors on phase speed prescribed in the boundary conditions can notably deteriorate radiation properties. The normal mode approach is introduced to identify coherent structures of propagation and their corresponding phase speed. A simple and robust multi-dimensional propagation scheme can easily be derived from polarization relations. The rotating case is more difficult but it is possible, to some extent, to get around the dependency of phase speed on wave frequency and to keep the non-rotating formulation of the PRM conditions almost unchanged. The PRM scheme being applied to field anomalies, the question of the background reference state is addressed. The latter can be used to introduce incoming waves across the open boundaries or, alternatively, to represent the low-frequency variability of the model itself. The consistency of the pressure and tracer boundary conditions is finally discussed.

## 1. Introduction

There are basically two ways to prescribe the boundary conditions in regional and coastal oceanic circulation models: radiation boundary conditions or absorbing boundary conditions. The FRS method (Martinsen and Engedahl, 1987) is an absorbing condition frequently used by oceanic modellers and still subject to new investigations (Lavelle and Thacker, 2008). This method has interesting radiation properties, however, a particular attention should be paid to the reference state toward which the modelled inner solution is relaxed. Indeed, though it is numerically stable, the FRS condition will not produce a realistic behaviour in the case of an outgoing flux if the reference state does not match the modelled inner solution. In order to counteract this, several authors suggested to use a nudging layer with a large time constant (Marc-hesiello et al., 2001). In a way, this type of boundary condition is similar to the FRS condition. The difference between the FRS meth-

od of Martinsen and Engedahl (1987) which imposes a strong constraint on the boundary solution, and that of Jensen (1998) and Marchesiello et al. (2001), who imposes a weak constraint, is that the region near the boundary in the latter cases acts like a zone of nudging towards an imposed boundary condition. Since this kind of hybrid radiation-absorbing boundary condition seems to be widely used, using sponge or nudging layers does not exempt modellers from improving the performances of radiation conditions.

In free surface ocean models, internal and external modes are generally computed separately (Blumberg and Mellor, 1987; Shchepetkin and McWilliams, 2005). Authors, therefore, distinguish between the boundary conditions schemes for the external mode (barotropic conditions) and for the internal mode (3D or baroclinic conditions). The Flather radiation boundary condition (Flather, 1976), or characteristics condition (for a detailed description, see Blayo and Debreu, 2005) is frequently used for the external mode. Beside its interesting stability and conservation properties (Marsaleix et al., 2006), it is indeed a way to include efficiently a forcing term, and is consequently particularly appropriate for tide simulations (Oey and Chen, 1992).

<footnote>* Corresponding author. Tel.: +33 5 61 33 27 63; fax: +33 5 61 33 27 90. E-mail address: marp@aeroobs-mip.fr (P. Marsaleix).</footnote>
When applied to the internal mode, the efficiency of such radiation conditions is less obvious, and, within the framework of baroclinic studies, absorbing boundary conditions were proposed as an alternative (Carter and Merrifield, 2007). One of the main reasons for that is related to the fact that internal waves phase speed is difficult to estimate precisely, whereas surface waves phase speed is unambiguous ($c = \sqrt{gh}$) under the hydrostatic assumption. The internal waves phase speed depends not only on the stratification but also on the propagation mode. Moreover, the phase speed associated with a given baroclinic mode is related to the signal frequency. Finally, the possible mode coupling induced by the bathymetric slopes limits anyway the thrust of the modal approach in terms of the propagation speed. Orlanski (1976) proposed a method based on the analysis of the solution near the open boundaries in order to compute a local and time-evolving value for the wave phase speed. However, further studies demonstrated the weakness of the performances (Palma and Matano, 1998) as well as the unstable nature (Treguier et al., 2001) of the Orlanski type methods. Consequently, radiation conditions based on a fixed phase speed, a-priori representative of the main propagation mode (the first baroclinic mode), may appear as a reasonable option (Kourafalou et al., 1996).

Whether the phase speed is defined a-priori or computed using the Orlanski method, the radiation conditions for the internal mode are usually based on a Sommerfeld type wave equation. Although Nycander et al. (2008) recently proposed a set of characteristic conditions for the barotropic and baroclinic modes, studies dealing with a three-dimensional version of the barotropic Flather condition are, to our knowledge, rather scarce. This may be due to several technical problems, the resolution of which constitutes the scope of our study.

We actually propose a set of open boundary conditions based on the relations of polarization of internal waves. This approach has some similarities with the upper boundary condition for atmospheric limited-height models proposed by Bougeault (1983) and Klemp and Durran (1983). The method is presented progressively. We start with a simple case, that is, a single wave propagating in a non-rotating ocean in the direction normal to the open boundary (Section 2). The stability of different possible numerical schemes is discussed in Section 3.

More general cases are then addressed. It is known that the dispersive nature of internal waves is barely compatible with simplistic hypothesis on phase speed, generally found in usual 3D boundary schemes (Bennett and Chua, 1994). This problem can be partly solved by higher order schemes (Higdon, 1994) but we note that Nycander and Döös (2003) recommended against using the second-order Higdon's condition. As discussed in Section 4, the modal approach (Jensen, 1998, 1993) is a reliable alternative.

Currently used radiation conditions assume that the waves propagate in a direction normal to the open boundaries and are therefore not relevant for realistic multidirectional studies. Approximate boundary conditions for multidirectional propagations have been proposed by Engquist and Mayda (1977). Raymond and Kuo (1984) developed a method based on a Sommerfeld condition for multidirectional cases, however, several authors (Marchesiello et al., 2001; Barnier et al., 1998) questioned its stability in some cases. A stable method aimed at adapting Flather condition to three-dimensional cases is presented in Section 5. Finally, the definition of the reference state can also be an obstacle when adapting Flather condition to three-dimensional cases. The impact of this question is double: it concerns the potentially active nature of the boundary condition in the case of incoming flux, and can be a source of errors in the case of outgoing flux if the low-frequency evolution of the ambient stratification is not correctly taken into account. Those two points are detailed in Section 6. Numerical tests are performed in order to estimate the performances of the

3D Flather condition developed in this paper. The limitations of this method, in particular those induced by the rotating case, are discussed in Section 7. The consistency of pressure and tracer boundary conditions is finally discussed in Section 8.

## 2. Basic schemes

### 2.1. The non-rotating case

In this study, we use the regional and coastal oceanic circulation model SYMPHONIE. This model, fully described in Marsaleix et al. (2008), uses the Boussinesq and hydrostatic approximations. In the first case, and for the sake of clarity, we use a simple configuration where propagation is in the $Oxz$ vertical plane. Practically, this means that we neglect hereafter the Coriolis term (note that the latter will be considered in Section 7) as well as the derivative with respect to $y$ in the model equations.

First, we consider the equation proposed by Gill (1982, Section 6.4):

$$ \frac{\partial u'}{\partial t} = - \frac{1}{\rho_0} \frac{\partial p'}{\partial x} \qquad (1) $$

where $(u', p')$ are the horizontal current and pressure perturbations associated with the wave. The total current and pressure are $(u, p) = (u' + u_{ref}, p' + p_{ref})$ where $(u_{ref}, p_{ref})$ corresponds to an ambient reference state evolving at a much smaller frequency than the perturbations. Then we consider a simple gravity wave propagating in the positive $x$ direction:

$$ (u', p') = (u_0, p_0) \times \cos(\omega t - kx) \qquad (2) $$

This wave must leave the domain through the open boundary located downstream the propagation, where we propose to apply a boundary condition that satisfies (1) and (2):

$$ p' = \rho_0 cu' \qquad (3) $$

where $c = \omega/k$ is the phase speed. A similar reasoning applied to the waves propagating in the negative $x$ direction provides the boundary condition at the other open boundary: $p' = -\rho_0 c u'$. Replacing $p'$ by the barotropic pressure $g\rho_0\eta$, $u'$ by $\bar{u}$ and $c$ by the surface wave phase speed in the long wave approximation, $\sqrt{gh}$, where $h$, $\eta$ and $\bar{u}$ correspond, respectively, to the bathymetry, the free surface elevation and to the mean transport, this condition is analogous to the classical Flather barotropic condition, $\bar{u} \pm (g/h)^{1/2}\eta = 0$ (Johns et al., 1983) with no relaxation term. Eq. (3) can also be seen as a simple deduction of the relation of polarization proposed by Gill (1982, p. 263, Eq. 8.5.3). In the following, this scheme will be referred as a polarization relation method (PRM).

We note that $u_{ref}$ is not necessarily negligible compared to typical phase speeds $c$ and thus Eq. (1) could have included an advection term involving the reference field, namely $\frac{\partial u'}{\partial t} = -u_{ref} \frac{\partial u'}{\partial x} - \frac{1}{\rho_0} \frac{\partial p'}{\partial x}$. Morel et al., 2008). This would have led to replacing the boundary condition (3) by $p' = \rho_0(c + u_{ref})u'$. However, we did not retain this formulation, as it is not compatible with the normal modes approach proposed in Section 4, where it is presumed that only $p'$ and $u'$ depend on $z$. Thus, using the PRM method requires that $u_{ref} \ll c$ to be valid.

### 2.2. Numerical implementation

The ocean model is based on a C-grid (Arakawa and Lamb, 1977). Fig. 1 shows the discrete distribution of the variables near the open boundary involved in the boundary condition defined by Eq. (3). In anticipation of the following sections dealing with fully three-dimensional cases, Fig. 1 also indicates the position of
<figure><img src="image_2.png" /><figcaption>Fig. 1. Horizontal location of the state variables on the C-grid. Black circles: pressure on inner grid points. White circles, pressure on open boundary grid points. Black triangles pointing to the right, Ox current component, u, on inner grid points. White triangles pointing to the right, Ox current component, u, on open boundary grid points. Black triangles pointing upward, Oy current component, v, on inner grid points. White triangles pointing upward, Oy current component, v, on open boundary grid points.</figcaption></figure>

the Oy-component of the current: v. Regarding the Ox direction, the momentum equation for the u component is computed between i = 3/2 and i = M - 1/2, and i = 1/2 and i = M + 1/2 correspond to the open boundaries. For now, the pressure at i = 1 and i = M is considered as located on the open boundary. This raises the question at which point the boundary condition (3) should be applied: should one compute $p'(i = M)$ or $u'(i = M + 1/2)$? As $p'$ is not a prognostic variable of the model, a straightforward reasoning would rather lead to compute an open boundary current. However, we chose the solution recommended by Marsaleix et al. (2006) for the barotropic condition, i.e. we apply the boundary condition at the pressure point. This recommendation stems from considerations associated with the model energy balance. This balance is obtained by multiplying (1) by $u'$ and by integrating the resulting equation over the whole domain. In practice, the discrete form of this balance is $\sum_{i=1,M-1} u'_{i+1/2} \frac{\partial u'_{i+1/2}}{\partial t} \Delta x = -\sum_{i=1,M-1} u'_{i+1/2} \frac{p'_{i+1} - p'_{i}}{\rho_0 \Delta x} \Delta x$, where the right-hand term is equal to the sum of a term of potential energy production, $P_{EP}$ (details in Marsaleix et al., 2008), and of a residual term corresponding to the energy flux through the open boundary: $\sum_{i=1,M-1} u'_{i+1/2} \frac{\partial u'_{i+1/2}}{\partial t} \Delta x = P_{EP} - \frac{p'_{M} u'_{M-1/2} - p'_{1} u'_{3/2}}{\rho_0}$. Applying the boundary condition on $p'_{M}$ and $p'_{1}$ enables the control of the discrete form of the energy flux through the open boundary and, therefore, to control the stability of the model. In particular, it clearly appears that the boundary conditions $p'_{1} = -\rho_{0} cu'_{3/2}$ and $p'_{M} = \rho_{0} cu'_{M-1/2}$ ensure that the net energy flux through the open boundaries, $p'_{1} u'_{3/2} - p'_{M} u'_{M-1/2}$, is systematically negative, which prevents the energy from increasing artificially inside the domain. Besides, we note that several authors (Shulman, 1997; Shulman et al., 1998) have demonstrated the interest of optimized Flather type conditions constraining the flux of energy through the open boundary.

Following a method proposed by Trefethen (1982), Nycander and Döös (2003) have demonstrated that a time lag between pressure and current (pressure being evaluated at t and current at t - Δt) enforces the numerical stability of the Flather condition when a leap-frog scheme is used to update the inner grid model equations (see also Marsaleix et al., 2008). We thus applied the following discretization:

$$p_1^t = -\rho_0 cu_{3/2}^{t-\Delta t} \quad (4a)$$

$$p_M^t = \rho_0 cu_{M-1/2}^{t-\Delta t} \quad (4b)$$

Accuracy and robustness of (4) are discussed in the following section.

# 3. PRM based modelling of internal gravity waves

## 3.1. Description of the numerical experiment

The PRM boundary conditions given by (4) are now implemented in order to study the propagation of internal gravity waves through a continuously stratified, non-rotating, ocean. We use a simple two-dimensional Oxz numerical domain with a constant depth of 5000 m. The initial density field, horizontally homogeneous, is such that the Brunt-Väisälä frequency, $N = (-g\rho_0^{-1}d\rho/dz)^{0.5}$, is equal to $1.4 \times 10^{-3} \text{ s}^{-1}$. The analytical theory (Gill, 1982) gives the vertical shape of baroclinic modes of propagation, as well as their respective phase speed. We use this information to replace radiation conditions at upstream open boundaries by Dirichlet conditions clamped onto incoming wave fields. Practically, the condition (4a) is replaced by a current boundary condition

$$u_{1/2}^t = u_0 \cos \left( \frac{\pi z}{h} \right) \sin(\omega t - kx_{1/2}), \quad (5)$$

with $u_0 = 0.01 \text{ ms}^{-1}$, the pressure $p'_{1}$ being thus computed through the inner domain equations. The frequency has been fixed to $\omega = 1.45 \times 10^{-4} \text{ s}^{-1}$ in order for the experiment to be representative of internal waves excited by semidiurnal tides. The phase speed, corresponding to the internal wave Mode 1, is $c \approx 2.2 \text{ ms}^{-1}$. Since boundary condition (4) is applied to variable anomalies (and not to global variables) the background reference state, $(u_{ref}, p_{ref})$ must be defined. For this first experiment, it is maintained constant throughout the simulation, and it is given by the initial fields. The main model parameters are given in Table 1. Details concerning numerical schemes of the model can be found in Marsaleix et al., 2008.

Along the downstream open boundary, the pressure anomaly is given by Eq. (4b). We can note, however, that applying a boundary condition on $p'_{M}$ does not exclude to define a boundary condition for $u'_{M+1/2}$, which is actually required for computation of the advective and diffusive fluxes in the momentum equation at the grid point node $i = M - 1/2$ (hereafter referred to as the flux condition). The flux condition can be an important issue in boundary layer modelling (Chini and Leibovich, 2003). Here, as we do not expect advection and diffusion to be the dominant processes, boundary value for $u'_{M+1/2}$ should hardly affect the momentum equation computed at $i = M - 1/2$. The latter should consequently be much more sensitive to the boundary conditions on the pressure anomaly $p'_{M}$ through the pressure gradient force. To remain consistent with the following sections of this paper, we chose to compute $u'_{M+1/2}$ with Sommerfeld conditions. A sensitivity test is moreover presented in a following section in order to clarify the influence of the flux condition.

Along these lines, tracer equations at the grid nodes $i = 1$ and $i = M$ require the knowledge of the outer domain temperature

<table><caption>**Table 1**
Model parameters for 2D Oxy simulations, small and large domains.</caption><tr><td>Δx</td><td>3000 m</td></tr><tr><td>Horizontal extent</td><td>L = 3000 km(large domain) L/2 = 1500 km(Small domain)</td></tr><tr><td>Δz</td><td>167 m</td></tr><tr><td>Bathymetry</td><td>Flat, 5000 m</td></tr><tr><td>Internal time step</td><td>216 s</td></tr><tr><td>Turbulence closure</td><td>Gaspar et al. (1990)</td></tr><tr><td>Initial state</td><td>u = 0, N = 1.4 × 10<sup>-3</sup> s<sup>-1</sup></td></tr><tr><td>Background field</td><td>Initial state</td></tr><tr><td>(u<sub>ref</sub>, p<sub>ref</sub>)</td><td></td></tr><tr><td>Coriolis parameter</td><td>f = 0</td></tr><tr><td>Duration of simulations</td><td>15 days</td></tr></table>
and salinity in the case of entering flux (Marsaleix et al., 2006). In the frame of the following experiments, the latter are simply given by the initial fields. Until Section 8 the pressure boundary conditions (4) remains independent of the local values of the temperature and salinity, but we note that a method is proposed in Section 8 to make them be consistent with each other.

Two domains are investigated (Table 1). The reference domain is 3000 km long and the simulation is 15 days long. The initial current is set to zero. As previously mentioned, internal waves (baroclinic Mode 1) are excited at the upstream boundary, corresponding to the grid node i = 1/2. Then the waves can propagate through the domain with the expected phase speed c = 2.2 ms⁻¹. Considering the time required by waves to cross the entire domain (roughly 16 days) this simulation is not affected by wave reflection processes on the downstream open boundary. To clearly investigate boundary condition radiation properties, a second, smaller configuration is investigated. It is basically identical to the reference simulation except that the numerical domain is 1500 km long, (half of the reference domain), so that the internal waves reach the downstream open boundary after about 8 days. The radiation properties of the boundary conditions are computed based on the global kinetic energy of the residual current. The latter is defined as the difference between the current computed over the smaller domain, u'_{small}, and the current of the reference simulation, u'_{large}, after 15 days, practically

$$ \int_{-h}^{0} \int_{0}^{L/2} \frac{(u'_{\text{large}} - u'_{\text{small}})^2}{2} dzdx, \quad (6) $$

where -h is the bottom depth and (0, L/2) are the horizontal locations of the upstream and downstream boundaries of the small domain. In the ideal case of perfect boundary conditions, this residual energy is zero. As previously mentioned, the choice of the phase speed can have a major impact on the baroclinic radiation conditions. The following numerical test consists in evidencing the sensitivity of radiation properties to errors on this phase speed. An ensemble of simulations has thus been performed with phase speed value used in Eq. (4b) ranging from 0 to 13.9 m/s (this upper value corresponding to an upper limit for Sommerfeld-type conditions considered in the next section). Fig. 2 shows the global kinetic energy of the final residual current as a function of the value of the phase speed specified in the boundary condition (4b): the low residual

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Residual kinetic energy as a function of the phase speed (in ms⁻¹). Curves are normalized by the residual kinetic energy obtained with a zero phase speed. Solid line, PRBM scheme. Dotted line, Sommerfeld equation. The straight line indicates the level of the residual energy obtained when the phase speed is computed according to the Orlanski method.</figcaption></figure>

energy values correspond to relatively good radiation properties and conversely high residual energy values correspond to relatively bad radiation properties. Curves are normalized by E₀, the residual kinetic energy (6) obtained with a zero phase speed. This case corresponds to the fully reflective “clamped” condition p'_{M} = 0. Besides, we have checked that E₀ is actually equal to the global kinetic energy of the second half of the reference domain, namely E₀ = ∫₀^L ∫_L/2 u'_{large}^2 dzdx. Fig. 2 shows, as expected, that the residual energy is minimum (about 1.1 × 10⁻³ × E₀) when the phase speed used in the boundary conditions is close to the theoretical value of 2.2 ms⁻¹. We also see that a 10%-error on phase speed increases the residual energy up to 4 × 10⁻³ × E₀.

## 3.2. Comparison with Sommerfeld radiation conditions

We now compare the PRM condition (4) to a classical Sommerfeld scheme of the form ∂u'/∂t ± c∂u'/∂x = 0. As recommended by Blayo and Debreu (2005), we applied these conditions to the current anomaly rather than to the global current. By doing this, we also maintain some consistency with the PRM conditions (4). We use the discrete form recommended by Orlanski (1976). At the downstream open boundary, the PRM condition is thus replaced by the following boundary condition on the current:

$$ \frac{u_{M+1/2}^{t+\Delta t} - u_{M+1/2}^{t-\Delta t}}{2\Delta t} = -c \frac{0.5(u_{M+1/2}^{t+\Delta t} + u_{M+1/2}^{t-\Delta t}) - u_{M-1/2}^{t}}{\Delta x} \quad (7) $$

The phase speed c can be fixed to an a priori value, or computed according to the Orlanski's (1976) method, that consists in reversing (7) at the neighbouring node and previous iteration where the current velocity field is known. Namely:

$$ c = - \frac{u_{M-1/2}^{t} - u_{M-1/2}^{t-2\Delta t}}{u_{M-1/2}^{t} + u_{M-1/2}^{t-2\Delta t} - 2u_{M-3/2}^{t-\Delta t}} \times \frac{\Delta x}{\Delta t} \quad (8) $$

Eq. (8) normally comes with a numerical stability condition, 0 < cΔt/Δx < 1. As far as the simulation parameters of Table 1 are concerned, this is equivalent to 0 < c < 13.9 ms⁻¹.

In a first implementation, the numerical experiments are based on (7) with a constant phase speed rather than on (7) and (8). As for the PRM scheme, an ensemble of simulations, each simulation differing from the others by the phase speed value used in (7), is now performed, in order to evidence the sensitivity of the radiation properties to errors on the phase speed. The latter ranges from 0 to 13.9 ms⁻¹, the lower and upper bounds of the aforementioned stability criterion. We also follow the same approach based on two domains with different spatial extent, and the kinetic energy of the residual current is systematically computed. Fig. 2 shows that the results are similar to those obtained with conditions (4). The minimum residual energy, corresponding to c = 2.2 ms⁻¹, is about 2.1 × 10⁻³ × E₀.

These tests include one simulation with phase speed computed according to the Orlanski's method (8). We note that the radiation properties of these boundary conditions (Fig. 2, straight dashed line) are rather disappointing since residual energy is about 0.07 × E₀. The simple Sommerfeld scheme remains better than the Orlanski scheme as long as the fixed phase speed used in (7) belongs to [1.5, 4] ms⁻¹, which represents quite a large relative error on the phase speed (30% and 80%, respectively, for the lower and upper bounds). Fig. 3 shows an example of celerity computed by Eq. (8) (solid line) and the same calculus but applied at an inner grid point, far enough from downstream boundary not to be perturbed by reflected waves (dashed line). The phase speed computed next to the boundary is far from equal to the expected theoretical value of 2.2 ms⁻¹. Although the instability problem reported by Treguier et al. (2001, p. 22,127) does not appear on Fig. 3, the phase speed computed next to the open boundary exhibits
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Phase speed (in ms⁻¹) computed according to the Orlanski method, as a function of time (in days) between Days 13 and 14. Solid line, the scheme is (as expected) computed on the inner grid point located next to the open boundary. Dotted line, the scheme is computed far from the open boundary so that the calculus are not affected by the wave reflections. The straight line indicates the expected analytical value: 2.2 ms⁻¹.</figcaption></figure>

abrupt changes of regime (negative values bounded by zero alternate with series of non zero values), revealing the complexity of dealing simultaneously with incoming and reflected waves. But the coexistence of incoming and reflected waves is not the only reason explaining the relative failure of the boundary conditions (8). Indeed, phase speed computed at the inner grid node (Fig. 3 dashed line) also shows large errors. The periodic occurrence of a division by zero every time (8) is computed on a wave crest, actually leads to a "blow-up" which is only partially controlled by the stability criterion.

### 3.3. Sensitivity to the flux condition

We now return to the PRM condition and test the influence of the flux condition. As previously mentioned, the advection and diffusion terms require a boundary condition (the flux condition) for $u'_{M+1/2}$, in addition to the PRM condition for $p'_M$. Here, this flux condition is simply given by the Sommerfeld condition (7) with an a-priori phase speed, equal to that used for the PRM condition. The numerical experiment of Section 3.1 is repeated but the Sommerfeld condition used for the flux condition is now replaced by a zero velocity condition, i.e. $u'_{M+1/2} = 0$. We note that this condition corresponds to a vanishing $c$ in (7). The PRM condition (4) is applied to $p'_M$ using a phase speed equal to the optimal value, i.e. $c - 2.2$ ms⁻¹. The residual energy is unchanged (about $1.1 \times 10^{-3} \times E_0$), showing that the flux condition has little influence on the model solution. This expected result is largely explained by the fact that the amplitude of the forcing current at the upstream boundary (5) is relatively small ($u_0 = 0.01$ ms⁻¹). The same experiment is thus repeated but using (5) with a stronger forcing current ($u_0 = 0.1$ ms⁻¹). The flux condition is given by the Sommerfeld equation. An ensemble of simulations is performed with a phase speed value ranging from 0 to $13.9$ ms⁻¹ in the flux condition, while the phase speed appearing in the PRM condition (4) is kept constant and equal to the optimal value $c = 2.2$ ms⁻¹. Fig. 4 shows the global kinetic energy of the residual current, as defined in (6), as a function of the phase speed used in (7). The curve is normalized by $E_0$, the residual energy obtained with a fully reflective condition (i.e. $u'_{M+1/2} = p'_M = 0$). As expected the minimum residual energy (about $1.2 \times 10^{-3} \times E_0$), is obtained when the phase speed used in the flux condition is close to the theoretical value $c = 2.2$ ms⁻¹. Now the sensitivity to the flux condition is more visible since the residual energy is about four times larger than the minimum residual energy when the flux condition uses $c = 0$. However, the radiation properties of the boundary conditions remain correct even in this worst case. When the phase speed is equal to its maximum value ($c = \Delta x / \Delta t = 13.9$ ms⁻¹), the Sommer-

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Residual kinetic energy as a function of the phase speed (in ms⁻¹) used in the flux condition (7). Note that the phase speed used in the pressure condition (4) is set to the optimal value $c = 2.2$ ms⁻¹. The curve is normalized by the residual kinetic energy obtained with a reflective condition, i.e. $c = 0$ in both flux and pressure conditions.</figcaption></figure>

field condition is equivalent to a zero-gradient condition. We note that in this case, the residual energy is still close to the minimum residual energy and thus we conclude that a zero-gradient condition may be a reasonable choice for the flux condition.

### 3.4. Truncation error

We now examine the discretization of the PRM condition. We note that the current and the pressure involved in (4) are evaluated neither at the same grid node nor at the same time. As underlined by Nycander and Döös (2003), this shortcoming has a little adverse effect on the model solution as long as the grid resolution and the time step are, respectively, much smaller than the length scale and the period of the outgoing waves. Generally, in this type of model, the time step is much smaller than the period of internal gravity waves, but the grid spacing is not negligible compared to the length scale of internal waves, notably if high baroclinic modes of propagation are expected to be present in the solution. A discretization estimating the pressure and the current at the same grid node and at the same time is now proposed and then compared to (4). This new discretization is based on a linear extrapolation
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Extrapolation of the inner grid velocity field used by scheme (9) at the downstream boundary. Step 1 (leaning dashed line): $\dot{u}_{M+1/2}^t$ is extrapolated from $u_{M-1/2}^{t-\Delta t}$ and $u_{M-3/2}^{t-2\Delta t}$, namely $\dot{u}_{M+1/2}^t = 2u_{M-1/2}^{t-M} - u_{M-3/2}^{t-2\Delta t}$. Step 2 (horizontal dashed line): $\dot{u}_M^t$ is interpolated from $u_{M-3/2}^t$ and $\dot{u}_{M+1/2}^t$, namely $\dot{u}_M^t = 0.75\dot{u}_{M+1/2}^t + 0.25u_{M-3/2}^{t-\Delta t} = 1.5u_{M-1/2}^{t-\Delta t} - 0.75u_{M-3/2}^{t-2\Delta t} + 0.25u_{M-3/2}^t$.</figcaption></figure>

of the inner grid velocity field, illustrated by Fig. 5. The new scheme is then given by:

$$p_1^t = -\rho_0 c (1.5 \times u_{3/2}^{t-\Delta t} - 0.75 \times u_{5/2}^{t-2\Delta t} + 0.25 \times u_{5/2}^t) \quad (9a)$$

$$p_M^t = \rho_0 c (1.5 \times u_{M-1/2}^{t-\Delta t} - 0.75 \times u_{M-3/2}^{t-2\Delta t} + 0.25 \times u_{M-3/2}^t) \quad (9b)$$

The previous numerical test is repeated. The residual energy is computed for an ensemble of simulations, each simulation differing from the others by the phase speed value used in Eq. (9). The results (given in Fig. 6) are compared to the same experiment performed with the PRM conditions (4). Fig. 6 shows significant improvement of radiation properties when the phase speed used in (9) is close to its theoretical value $2.2 \text{ ms}^{-1}$. Indeed, the minimum residual energy is $5 \times 10^{-6} \times E_0$, that is about two hundred time smaller than the minimum residual energy obtained with (4). However, we noticed that this new scheme is less stable than (4) for high values of the

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Residual kinetic energy as a function of the phase speed (in ms⁻¹). Curves are normalized by the residual kinetic energy obtained with a zero phase speed. The open boundary conditions are given by the extrapolation method (9b) (dotted line) or by the basic scheme (4b) (solid line).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Maximum permitted phase speed as a function of the time step. Solid line, the PRM condition is given by (4). Dashed line, the PRM condition is given by (9).</figcaption></figure>

prescribed phase speed. The stability condition for (4) and (9) (i.e. the highest permitted prescribed phase speed) is shown in Fig. 7. These curves were obtained from an ensemble of runs using different model time steps $\Delta t$. For a given $\Delta t$, increasing values of the prescribed phase speed in (4) or (9) have been tested until the inner solution becomes unstable (critical values given by Fig. 7). For a wide range of time steps, Fig. 7 shows that (4) remains stable as long as the prescribed phase speed does not exceed $\Delta x/\Delta t$. When $\Delta t$ becomes close to its maximum theoretical value, the highest prescribed phase speed is a bit smaller. We note that the highest possible time step is given by $\Delta t_{max} = 0.5\Delta x/c_1$ where $c_1 = 2.2 \text{ ms}^{-1}$ corresponds to the phase speed of the first baroclinic mode (see Nycander and Döös, 2003 for further details on the C.F.L. condition). In the case of (9), the condition of stability is globally obtained for $c < \Delta x/(1.5\Delta t)$, in other words, the highest prescribed phase speeds permitted by (9) are roughly 1.5 times smaller than those allowed by (4). Indeed, as the first inner current grid point is particularly sensitive to boundary instabilities, the RHS of (9b) is likely dominated by $u_{M-1/2}^{t-\Delta t}$ under unstable conditions, leading (9b) to be equivalent to $p_M^t = 1.5\rho_0 c u_{M-1/2}^{t-\Delta t}$ in this special case, i.e. 1.5 times the RHS of (4). Thus we conclude that the scheme (9) is a better option than the scheme (4) (since the former is more accurate than the latter) as long as we use a time step that satisfies the following stability condition: $\Delta t < \Delta x/(1.5c_1)$.

Our conclusions are however tempered by the fact that results obtained with (9) show relatively small differences with those obtained with (4) when phase speed differs significantly from its theoretical value. The benefit of the improved discrete form (9) becomes negligible if we make a large error on phase speed. Thus, in the following section, we focus on an accurate determination of the phase speed which appears as a key parameter.

# 4. Normal modes approach

## 4.1. Description

This approach, deduced from the theory of long waves propagating through a continuously stratified ocean (see Gill, 1982 (p. 159); Blayo and Debreu, 2005), is rather well adapted to the basic equations of our hydrostatic ocean model.

First of all, a vertical normal mode decomposition of the model fields leads to $w' = \dot{w}(z)\dot{w}(x, t)$, $u' = \dot{u}(z)\dot{u}(x, t)$ and $p' = \dot{p}(z)\dot{p}(x, t)$ where $w'$, the vertical velocity anomaly, is related to the horizontal current through the continuity equation:
$$ \frac{\partial u'}{\partial x} + \frac{\partial w'}{\partial z} = 0 \qquad (10) $$

Taking (1) and (10) into account, the hydrostatic equilibrium, $\frac{\partial \rho'}{\partial z} = -g\rho'$, and the hypothesis that $\rho'$, the density anomaly, is mainly driven by vertical advection of the background density field $\rho$, i.e. $\frac{\partial \rho'}{\partial t} + w' \frac{d\rho}{dz} = 0$, we obtain (see Gill, 1982, for details):

$$ \frac{d^2 \hat{w}}{dz^2} + \frac{N^2}{c_e^2} \hat{w} = 0 \qquad (11) $$

with surface and bottom boundary conditions $\hat{w}_{z=0} = \hat{w}_{z=-h} = 0$. N is the Brunt-Väisälä frequency. Eigenvalues and eigenfunctions of this problem give the possible values of phase speed $c_e$ and their corresponding vertical mode. Although an analytical solution is available in this constant N case, more general density profiles require a numerical solution. This can be performed on the model grid itself. If n is the number of vertical levels, including the two levels at the surface and at the bottom, where the boundary condition $\hat{w} = 0$ is applied, we can theoretically obtain the baroclinic phase speeds $c_q$ of the first n-2 modes ($q = 1, n-2$), arranged in decreasing order, and their corresponding vertical velocity profile $\hat{w}_q(z)$. The number of modes accessible to the model is normally much smaller since a correct representation of the periodic signals requires a minimum number of grid nodes per wavelength. As far as the grid characteristics of Table 1 are concerned, Marsaleix et al. (2008) showed that the first seven baroclinic modes of propagation can reasonably be modelled. Fig. 8 shows the first three baroclinic vertical velocity profiles and the associated phase speeds, obtained from a numerical solution of (11). We note that these results are confirmed by the analytical solution: $\hat{w}_q = \sin(q\pi z/h)$, $c_q = 2.2/q$ ms$^{-1}$. Taking (1) and (10) into account, vertical profiles for other variables are defined as $\hat{u}_q(z) = \hat{p}_q(z) = d\hat{w}_q/dz$.

The next step consists in a decomposition of the form, $u'(x, z, t) = \sum_{q=1, n'} [\hat{u}_q(z)\hat{u}_q(x, t)] + \varepsilon$, where $n^* \le n-2$ is the number of the last normal mode retained, $\hat{u}_q$ are the amplitudes of the projection,

$$ \hat{u}_q(x, t) = \left( \int_{-h}^{0} \hat{u}_q^2(z) dz \right)^{-1} \times \int_{-h}^{0} \hat{u}_q(z) u'(x, z, t) dz \qquad (12) $$

and $\varepsilon$ is a residual current. The normal mode method enables to adapt the phase speed in the boundary conditions (3) to each of the vertical modes composing the full solution. In practice, the PRM boundary condition (3) becomes:

<figure><img src="image_9.png" /><figcaption>Fig. 8. The first 3 baroclinic normal modes for vertical velocity and their corresponding phase speeds.</figcaption></figure>

$$ p' = \pm \sum_{q=1, n^*} \rho_0 c_q \hat{u}_q \hat{u}_q \qquad (13) $$

Concerning numerical implementation, we use the same extrapolation form as in Eq. (9). As an example, the horizontal velocity $\hat{u}_q$ at the open downstream boundary is given by

$$ \left( \int_{-h}^{0} \hat{u}_q^2(z) dz \right)^{-1} \times \int_{-h}^{0} \hat{u}_q(z) \left( 1.5 \times u_{M-1/2}^{t-\Delta t} - 0.75 \times u_{M-3/2}^{t-2\Delta t} + 0.25 \times u_{M-3/2}^{t} \right) dz \qquad (14) $$

## 4.2. Limitations of the linear theory

The experiment of Section 3.1 is repeated but the downstream boundary condition is now given by the modal condition (13). Vertical profile and phase speed of baroclinic modes are estimated according to (11). Considering that the forcing wave at the upstream boundary corresponds to Mode 1, we start with a simple configuration of (11) for which the current is projected on the first mode only ($n^* = 1$). As the PRM condition is theoretically limited to cases of small wave amplitude cases, we performed an ensemble of simulations in order to evidence the sensitivity of the radiation properties of (13) to $u_0$ the amplitude of the forcing current at the upstream boundary (5). As in Section 3, the kinetic energy of the reflected current is computed according to (6) and then normalized by $E_0$, the kinetic energy of the reflected current obtained with the fully reflecting boundary condition. Fig. 9 shows that the radiation properties remain satisfactory (reflected energy lower than $10^{-3} \times E_0$) as long as $u_0$ is lower than 0.4 ms$^{-1}$ (about 15% of the wave phase speed). Beyond this value, the boundary condition does not perform correctly (results not shown). When $u_0 = 0.01$ ms$^{-1}$, the maximum vertical displacement of isopycnals at mid-height of the water column is about 15 m. When $u_0 = 0.4$ ms$^{-1}$, the amplitude of the internal wave is of the order of 600 m, quite a strong value if one considers internal wave amplitudes observed in the Bay of Biscay, a region of the North Eastern Atlantic known for the strength of its internal tides (Pichon and Correard, 2006).

We observed that changing the Brunt-Väisälä frequency had an impact on the maximum amplitude $u_0$ that can be applied without

<figure><img src="image_10.png" /><figcaption>Fig. 9. Residual kinetic energy as a function of the upstream forcing current $u_0$ (in ms$^{-1}$). Curves are normalized by the residual kinetic energy obtained with a reflective condition. The open boundary conditions are given by (13) using $n^* = 1$ (solid line) or using $n^* = 6$ (dashed line).</figcaption></figure>
unreasonable deterioration of the radiation properties. In our numerical experiments (results not shown), higher phase speeds enable larger amplitudes, leading us to conclude that the PRM condition performs well as long as the amplitude $u_0$ remains one order of magnitude smaller than the wave phase speed.

Simulations also show that high values of $u_0$ tend to excite higher baroclinic modes through non-linear processes and unsurprisingly tests reveal that boundary condition (13) performs better if current is projected on several modes rather than on Mode 1 only (Fig. 9).

### 4.3. A multimodal propagation test

In this section, the waves excited at the upstream boundary include the first three baroclinic modes, namely $u_{1/2}^t = \sum_{q=-1,3} u_0^q \cos(\frac{q\pi z}{h}) \sin(\omega t - kx_{1/2})$. As these modes do not propagate at the same speed, ($c_1 \approx 2.2 \text{ ms}^{-1}$, $c_2 \approx 1.1 \text{ ms}^{-1}$, $c_3 \approx 0.7 \text{ ms}^{-1}$), and as the duration of the reference simulation is limited in order to avoid wave reflections on the downstream boundary, highest (slowest) modes are excited sooner. The amplitudes are defined as:

$$ 
\begin{aligned}
u_0^3 &= u_0 \\
u_0^2 &= 0 \text{ if } t < 8 \text{ days, } u_0^2 = u_0 \text{ if } t \ge 8 \text{ days} \\
u_0^1 &= 0 \text{ if } t < 16 \text{ days, } u_0^1 = u_0 \text{ if } t \ge 16 \text{ days}
\end{aligned} 
\quad (15) 
$$

so that the three modes simultaneously reach after 24 days the middle of the reference (large) domain (Fig. 10), which also corresponds to the downstream boundary of the small domain. Then, the simulation lasts 8 more days, the additional time required for the fastest mode (Mode 1) to cross the second half of the large domain and reach the downstream boundary. The amplitude $u_0$ is set to $0.01 \text{ ms}^{-1}$. The residual energy obtained with boundary conditions following the normal mode approach, (13), is now compared to the residual energy obtained with the global approach (9). Concerning the latter, an ensemble of runs has been performed with $c$, the phase speed appearing in (9), ranging from 0 to $3 \text{ ms}^{-1}$. The residual energy has been computed for the global residual current but also for its modal components. As expected, we can see in Fig. 11 that the residual energy of the separated modes is minimum when the celerity used in Eq. (9) is equal to the theoretical phase speed of the corresponding mode. The global residual energy is not satisfying because none of the single celerity values really suits the three modes simultaneously. The global minimum, obtained somewhere between $c_1$ and $c_3$ ($c = 1.65 \text{ ms}^{-1}$ in our case), represents 5% of $E_0$, the residual energy obtained with fully reflective conditions (that is with $c = 0$). Use of boundary conditions based on the normal mode approach leads to a residual energy which is two orders of magnitude smaller, namely about $5 \times 10^{-4} \times E_0$. This result is obviously very promising but not as small as in the experiment shown in Section 3.3, for which the residual energy is about

<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Amplitudes (in ms$^{-1}$) $\dot{u}_i(x, t)$, of the normal mode decomposition, $u' = u_i(z)\dot{u}_i(x, t)$, over the first half of the large domain, during the first 24 days. a: Mode 1, b: Mode 2, c: Mode 3.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Fig. 11. Residual kinetic energy as a function of the phase speed (in ms⁻¹). Curves are normalized by the residual kinetic energy obtained with a zero phase speed. Solid line, with full current; *, using current Mode 1 component only; x, with Mode 2 component only; +, with Mode 3 component only.</figcaption></figure>

$5 \times 10^{-6} \times E_0$. This can be explained by the fact that the modal decomposition is not perfect and leads to a non-zero residual current, $\epsilon$, which is not radiated by boundary conditions (13).

### 4.4. Effect of the bottom slope

The hypothesis that each mode propagates independently from the others with a given phase speed, requires a constant bathymetry. If not, bottom slope is likely to induce a coupling effect between modes. A probable consequence is a decrease of radiation properties of boundary condition (13). We now present a numerical experiment, similar to the previous one, but including a gradient of bathymetry in the vicinity of the open boundary of the small domain. The bathymetry is given by:

$$ \begin{aligned} h &= h_0 \quad 0 < x < L/2 - 16\Delta x \\ h &= h_0 - \alpha(x - (L/2 - 16\Delta x)) \quad L/2 - 16\Delta x < x < L/2 \end{aligned} \quad (16) $$

where $x = 0$ and $x = L/2 = 1500 \times 10^3$ m, respectively, correspond to the upstream and downstream boundaries, $\alpha$ is the bottom slope at the downstream boundary (other parameters are given in Table 1). Bathymetry is thus constant ($h = h_0 = 5000$ m) in the upstream part of the domain and varies linearly in the last 48 km of the downstream part of the domain. The bathymetry of the large domain, used for the simulation of reference, is the same as the one used for the small domain for $0 < x < L/2$. The domain is flat for $x > L/2$, i.e.:

$$ h = h_0 - \alpha 16 \Delta x \quad x > L/2 \quad (17) $$

At the upstream boundary, the first 3 baroclinic modes are excited according to the procedure described in Section 4.2. An ensemble of simulations is performed in order to evidence the sensitivity of the radiation properties of (13) to the bottom slope $\alpha$. For a given value of $\alpha$, three runs are carried out. The first run is the reference simulation over the large domain. The second run is computed over the small domain, boundary condition (13) being applied at the downstream boundary. The kinetic energy of the reflected current is computed according to (6). The latter is normalized by $E_0$, the kinetic energy of the reflected current obtained with a fully reflecting condition applied at the downstream boundary of the small domain (the third run), $p'_M = u'_{M+1/2} = 0$.

<figure><img src="image_13.png" /><figcaption>Fig. 12. Residual kinetic energy as a function of the bottom slope $\alpha$ used in (16) and (17).</figcaption></figure>

Fig. 12 shows the evolution of the normalized energy of the reflected current as a function of $\alpha$. As expected, the smallest value, $5 \times 10^{-4} \times E_0$, is obtained for $\alpha = 0$. The kinetic energy of the reflected signal increases with $\alpha$ but remains smaller than $10^{-2} \times E_0$ as long as $\alpha < 0.05$. In other words, the radiation properties of (13) are satisfactory for a wide range of topographic slopes. For larger values of $\alpha$, possibly encountered in coastal zones (for instance over shelf breaks), the radiation properties rapidly worsen.

When $\alpha > 0.07$, results (not presented) are more difficult to analyse, as a new regime of propagation is evidenced by the reference (large domain) simulation, waves being partly reflected by the bottom slope, as it approaches the critical value. Thus, the energy of the waves reaching the top of the slope significantly decreases as $\alpha$ gets close to 0.1 and the computation of the kinetic energy of waves reflected by the open boundary condition of the small domain at $x = L/2$ becomes meaningless.

## 5. Cases of multi-dimensional propagation

### 5.1. Basic equation

The schemes presented in the previous subsections are based on the assumption that the wave propagation is horizontal and locally perpendicular to the open boundaries. This is obviously not a realistic hypothesis, which has potentially negative effects on the efficiency of the open boundary conditions, since ignoring the angle of incidence at which waves cross the open boundary is similar to using a wrong apparent phase velocity in the direction normal to the boundary. As far as Sommerfeld type conditions are concerned, the Orlanski scheme seems well adapted to the computation of the normal phase velocity but we saw in Section 3 that the accuracy of the Orlanski scheme is in fact strongly limited. We also note the existence of Orlanski type methods (Raymond and Kuo, 1984; Marchesiello et al., 2001) apparently better suited to multi-dimensional wave propagation.

We now test a multi-dimensional version of the PRM scheme proposed in the previous subsection. First of all, we note that the multi-dimensional problem is simplified by the normal mode approach of Section 4. Indeed, when decomposing the solution in normal modes propagating horizontally in a two-dimensional space, the problem of the vertical propagation naturally disappears. Thus, let us consider a wave that propagates in an arbitrary horizontal direction. We consider the open boundary that is perpendicular to the Ox axes, which points outward. The pressure
anomaly is given by $p' = p_0 \sin(kx + ly - \omega t)$ and since only waves leaving the domain are considered, we have $k > 0$. We continue to neglect Coriolis accelerations and, thus, current components are simply deduced from the momentum Eq. (1) which is now associated with a momentum equation for the y component of the current, namely $\partial v'/\partial t = -(1/\rho_0)\partial p'/\partial y$. We obtain (Gill, 1982, p. 263, Eq. 8.5.3) $(u', v') = (\frac{kp'}{\rho_0 \omega}, \frac{p'}{\rho_0 \omega})$ so that we have $p'^2 = (\rho_0 c U)^2$, where $c = \sqrt{\omega^2/(k^2 + l^2)}$ is the horizontal phase velocity and $U = \sqrt{u'^2 + v'^2}$ is the modulus of the current anomaly. At this stage, we may simply deduce that $|p'| = \rho_0 c U$ but we also know that $u'$ and $p'$ share the same sign (since $k > 0$) and consequently the ambiguity on the pressure sign is easily removed if the current anomaly is known. We finally obtain a boundary condition for the pressure anomaly, $p' = \rho_0 c \operatorname{sign}(u') \times U$, where $\operatorname{sign}(u')$ is defined as follow: $\operatorname{sign}(u') = 1$ if $u' > 0$, $\operatorname{sign}(u') = -1$ if $u' \le 0$. On the second boundary perpendicular to the Ox axes, but with the latter pointing inward, outgoing waves are characterised by $k < 0$ and the boundary condition is given by $p' = -\rho_0 c \operatorname{sign}(u') \times U$. In a similar way, the pressure anomalies on boundaries perpendicular to the Oy axes, are computed according to $p' = \pm \rho_0 c \operatorname{sign}(v') \times U$.

## 5.2. Discretization

Let us consider the boundary normal to the Ox direction on which a boundary condition for $p_M^t$ is computed. A value for the Ox current component, $u'$, at pressure grid point $(M, j)$, is given by the extrapolation of the inner current values, as in Eq. (9b). A similar extrapolation of the tangential component, $v'$, but taking into account the distribution of the state variables on the C-grid (Fig. 1), leads to: $0.5 \times (2 \times v_{M-1,j-1/2}^{t-\Delta t} - v_{M-2,j-1/2}^{t-2\Delta t} + 2 \times v_{M-1,j+1/2}^{t-\Delta t} - v_{M-2,j+1/2}^{t-2\Delta t})$.

Concerning the boundary condition suggested in Section 5.1, we note the possibility of a discontinuity, in extreme cases for which the direction of propagation is parallel to the open boundary. In such situations, small fluctuations of a quasi-zero normal current are likely to abruptly change the sign of the pressure anomaly. In order to counter this possible source of numerical instabilities, we constrain the tangential current to vanish ($^t$) when the normal current becomes close to zero. This is equivalent to prescribing a zero pressure anomaly in the case of a parallel propagation. As far as the two boundaries, normal to the Ox axes, are concerned, we ensure that the tangential current $v$ is small enough compared to the normal current. Practically, we impose $v^2 < \beta u^2$ ($^t$), where $\beta$ is a positive number. The tangential current is not modified as long as the angle of propagation with the normal direction, $\theta = \tan^{-1}(v/u)$, does not exceed the limit angle $\theta_{max} = \tan^{-1}(\sqrt{\beta})$, (i.e. $-\theta_{max} < 0 < \theta_{max}$) and, conversely, $v$ is progressively attenuated as $\theta$ approaches $\pi/2$ or $-\pi/2$. A "controlled incidence angle" scheme is finally expressed as:

$$p' = \pm \rho_0 c \operatorname{sign}(u') \times (u'^2 + \min(\beta u'^2, v'^2))^{1/2} \quad (18)$$

where $\min(\beta u'^2, v'^2) = v^2$ if $\beta u'^2 \ge v^2$ and $\min(\beta u'^2, v'^2) = \beta u'^2$ if $\beta u'^2 < v^2$. The larger $\beta$ is, the less the tangential current is constrained. In that case the boundary condition aims at permitting a wide range of propagation directions, but the risks: of developing the aforementioned instabilities are high. In contrast, Eq. (18) tends to be equivalent to the basic scheme (9) when $\beta$ is small; the scheme is more stable but radiation properties are partially lost if the direction of the wave propagation is sensibly different from

<table><caption>**Table 2**
Model parameters for 3D Oxyz simulations, small and large domains.</caption><tr><td>Δx × Δy</td><td>3000 × 3000 m</td></tr><tr><td>Horizontal extent</td><td>900 × 900 km (large domain) 300 × 300 km (small domain)</td></tr><tr><td>Duration of simulation</td><td>80 h<sup>a</sup></td></tr><tr><td>Other parameters</td><td>As in Table 1</td></tr></table>

<sup>a</sup> Time required for waves to reach the downstream boundary of the large domain (Fig. 13a).

the Ox direction. If $\beta = 1$ waves are correctly radiated as long as their angle of propagation with the normal direction does not exceed $45^\circ$. Analogously, the boundary conditions at the boundaries normal to the Oy axes are given by:

$$p' = \pm \rho_0 c \operatorname{sign}(v') \times (v'^2 + \min(u'^2, \beta v'^2))^{1/2} \quad (19)$$

## 5.3. Numerical tests

The numerical experiment is similar to the one of Section 3 but we now use a three-dimensional domain in order to simulate internal waves crossing the outgoing boundaries with an oblique angle of incidence. A non rotating regime is still assumed. The spatial extent of the new domains and the duration of simulations are indicated in Table 2.

Among the four boundaries, two of them (Ox and Oy axes pointing locally outward) use a radiation condition while the two others are fixed to the value of the entering waves. As in Section 3, the latter, corresponding to baroclinic propagation Mode 1, are given by:

$$\begin{align}
u' &= R u_0 \cos \left(\frac{-\pi z}{H}\right) \sin(\omega t - k_x x - k_y y) \nonumber \\
v' &= R v_0 \cos \left(\frac{-\pi z}{H}\right) \sin(\omega t - k_x x - k_y y) \tag{20}
\end{align}$$

with $(k_x, k_y) = \frac{\omega}{c}(\cos(\theta), \sin(\theta))$, $(u_0, v_0) = A(\cos(\theta), \sin(\theta))$ where $\theta$ is the angle made by the wave propagation with the Ox axes. The period is fixed to 12 h, the phase speed is $2.2 \text{ ms}^{-1}$, current amplitude is $A = 0.01 \text{ m}$ and $R$ is a function representing the propagation of the frontal wave: $R = 1$ if $(k_x x - k_y y)/(k_x^2 + k_y^2)^{1/2} \le ct$ and $R = 0$ if $(k_x x - k_y y)/(k_x^2 + k_y^2)^{1/2} > ct$. Lastly, $(t, x, y)$, respectively, measures time from the beginning of the simulation and the distances along Ox and Oy axes, whose origin corresponds to the upstream corner of the propagation (Fig. 13a). The angle of propagation $\theta$ is fixed to $45^\circ$. As for the two-dimensional experiment of Section 3 we use two domains: on the one hand, a large domain providing a reference solution, shown by Fig. 13a and b, that is not affected by the possible wave reflections over the downstream boundaries and on the other hand, a small domain revealing the radiation properties of the boundary conditions (19) and (18).

The numerical experiment consists in measuring the energy of the residual current (defined as in Section 3) for an ensemble of simulations, each simulation differing from the others by the value of the $\beta$ parameter used in (18) and (19). Practically, $\beta$ ranges from 0 to $\infty$ (i.e. $\theta_{max} = \tan^{-1}(\sqrt{\beta})$ ranges from 0 to $90^\circ$). The residual energy, shown by Fig. 14, is normalized by $E_0$, the residual energy obtained for the fully reflective downstream boundary condition $p' = 0$. The less favourable case corresponds to $\beta = 0$, with a residual energy of about $0.015 \times E_0$. As mentioned previously, this particular case makes the boundary condition (18) and (19) equivalent to the unidirectional scheme (9). The residual energy diminishes, as $\beta$ increases, and falls down to a minimum value of the order of $2.9 \times 10^{-5} \times E_0$ for $\beta \ge 1$, that is $\tan^{-1}(\sqrt{\beta}) > 45^\circ$. This result is consistent with the $45^\circ$ angle at which modelled waves actually cross the open boundaries. This minimum residual energy can be compared to the numerical test with a normal boundary propagation, presented in Section 3.3 (we obtained $5 \times 10^{-6} \times E_0$). The per-

<footnote>* The value used by the boundary condition, not the one used by the inner equations.</footnote>
<figure><img src="image_14.png" /><figcaption>**Fig. 13.** (a) Surface currents (ms⁻¹) after 80 h computed with the large domain. Solid lines indicate the spatial extent of the small domain. The dashed line indicates the location of the vertical section used to compute (b). (b) Vertical section (location shown in (a)) of currents (ms⁻¹) after 80 h.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Residual kinetic energy for a multi-dimensional propagation case (angle of incidence equal to 45°), as a function of the limit angle $\vartheta_{max} = \tan^{-1}(\sqrt{\beta})$, in degree, where $\beta$ is the parameter appearing in (18) and (19). The curve is normalized by $E_0$, the residual kinetic energy obtained for the fully reflective condition $p' = 0$.</figcaption></figure>

formances are comparable, but as expected, the normal boundary propagation case is better since it is not penalized by additional computations induced by the tangential boundary current.

The case of the parallel propagation can be viewed as an asymptotic regime, beyond which the waves can be considered as entering

the domain and consequently the boundary conditions considered as active conditions. In this latter case, the boundary condition performance largely depends on a correct definition of the reference state, ($u_{ref}$, $p_{ref}$), mentioned in Section 2.1.

## 6. The reference state

### 6.1. Low-frequency signal

In the previous sections, the reference state was presented as a way to distinguish low- and high-frequency motion. We assumed that the high frequencies are composed of (incident) gravity waves, whose free propagation out of the domain is the primary purpose of the PRM boundary conditions. However, the efficiency of the boundary conditions does not only depend on the radiation properties of the scheme proposed in the previous sections, but also on whether the low-frequency component of the model solution is faithfully represented by the reference state. In an embedded modelling approach, a simulation over a larger domain can be performed in order to provide the external forcing terms at the open boundaries of the nested model. In the case of a one-way nesting method, the solution in the larger domain is computed first, enabling an a posteriori calculation of its low-frequency component in order to define ($u_{ref}$, $p_{ref}$). An example of such an approach is the Northern Aegean Sea model of Kourafalou and Tsiaras (2006) forced by the daily averaged outputs of a larger and coarser model of the Eastern Mediterranean Sea (Sofianos et al., 2006) in the frame of the MFSPP project (Pinardi et al., 2003). In this kind of configuration, the radiation properties of the boundary conditions depend, among other factors, on the consistency between the
low-frequency component of the model solution and its reference counterpart. However, inconsistencies may indeed occur, notably in case of outgoing circulation. This can be avoided by using a nudging layer inside which the model solution is slowly relaxed toward the reference field. The time scale associated with the restoring force should be large enough to avoid interference with gravity waves and spurious frontal patterns near the open boundaries (Marchesiello et al., 2001).

## 6.2. High-frequency signal

Conversely, we can also consider the case where the reference field is used to introduce incident waves into the domain. As a matter of fact, if radiation schemes are primarily designed to let disturbances out rather than imposing incoming waves, a major interest of Flather type schemes is also to provide an efficient way to force regional or coastal models with external fields, the latter being for instance obtained from a larger scale model (Oey and Chen, 1992; Nycander and Döös, 2003). Concerning this particular problem, Sommerfeld type conditions are generally not very convincing. Different schemes, based on the sign of the Orlanski's phase speed (8), have been proposed to distinguish outgoing and incoming processes in order to replace, in this latter case, the wave equation by a stronger constraint, for instance a Dirichlet condition (Miyakoda and Rosati, 1977). Discontinuities of regimes, mentioned in Section 3.3, can be avoided by combining a wave equation with a relaxation term toward the external field (Marchesiello et al., 2001) but this is a weak constraint that does not necessarily prevent discrepancies between the model fields and their external counterparts (Blumberg and Kantha, 1985) nor guarantees the consistency of the boundary conditions with external fields (Marsaleix et al., 2006).

The following experiment tests the ability of the PRM condition to let the outgoing waves radiate out and, simultaneously, the incoming waves enter into the domain. As in Section 3, we use a two-dimensional Oxy numerical domain with a constant depth of 5000 m (see also Table 1). The boundary conditions are computed using the multimodal scheme (13). The experiment is summarized by Fig. 15. Incoming waves are excited at the two open boundaries. In order to distinguish waves propagating in the positive x direction from those propagating in the opposite direction, the former are excited according to baroclinic Mode 2, whereas the latter are given by baroclinic Mode 1. Thus, at the $x = L/2$ boundary the reference field is given by $u_{ref} = u_0 \cos(\frac{\pi x}{h}) \sin(\omega t)$ where, as in previous sections, the amplitude of the horizontal current is set to $u_0 = 0.01 \text{ ms}^{-1}$ and the period, $2\pi\omega^{-1}$, equals 12 h. Assuming

that incoming waves travel in the negative x direction at the phase speed of Mode 1, $c_1 = 2.2 \text{ ms}^{-1}$, the reference pressure is simply $p_{ref} = P_b(z) - \rho_0 c_1 u_{ref}$, with $p_b(z)$ the background pressure field at rest. At the $x = 0$ boundary, incoming waves travelling in the positive x direction at the phase speed of Mode 2, $c_2 = 1.1 \text{ ms}^{-1}$, are given by $u_{ref} = u_0 \cos(\frac{2\pi x}{h}) \sin(\omega t)$ and $p_{ref} = P_b(z) + \rho_0 c_2 u_{ref}$. Because Mode 1 travels faster than Mode 2, the latter is excited first (8 days before Mode 1) so that the two modes reach the outgoing boundaries after 16 days. The simulation is extended 8 more days, during which boundary conditions have, at the same time, to radiate out outgoing waves and excite incoming waves. The interest of this experiment is to show that the PRM conditions can simultaneously act as passive and active conditions (see Marchesiello et al., 2001, for a detailed study of these two concepts), in contrast with Orlanski type conditions which require the choice of one of the two regimes.

The results are compared with a reference solution. The latter is provided by two runs, separately computing the propagation of the two modes over a numerical domain that is large enough to avoid any spurious reflection at the downstream boundary (Fig. 15). These two reference runs impose the incoming waves at their upstream boundary through a Dirichlet condition.

The full velocity field computed over the small domain is projected over its two first baroclinic components. Each of these two components is then compared to the corresponding reference solution. As in the previous sections, errors are given by the kinetic energy of the residual current computed according to (6). Here, the errors are $1 \times 10^{-4} \times E_0^1$ and $2 \times 10^{-3} \times E_0^2$ for Modes 1 and 2, respectively, where, similarly to Section 3.1, $E_0^1$ and $E_0^2$ are the kinetic energy computed over the second half of the large domain, for Modes 1 and 2, respectively, (Fig. 15). These small errors indicate that the PRM condition performs as well as the Dirichlet condition as far as the forcing of incoming waves is concerned and, at the same time, maintains its good radiation properties.

# 7. The rotating case

## 7.1. Limitations of the rotating case

The boundary conditions proposed in the previous sections are based on a rather restrictive hypothesis. While the hydrostatic regime is consistent with the basic equations of our model, omission of the Coriolis effects is more problematic. At frequencies close to the inertial frequency, the wave propagation features are sensibly different from what we hitherto assumed. The wave propagation is notably partially vertical, while its angle with the vertical axes

<figure><img src="image_16.png" /><figcaption>**Fig. 15.** Middle panel, small domain using PRM boundary conditions, forced by a baroclinic Mode 1 wave at $x = L/2$ and a baroclinic Mode 2 incoming wave at $x = 0$, as obtained from the reference fields. Top panel, large domain providing the reference solution for the baroclinic Mode 2 wave. Bottom panel, large domain providing the reference solution for the baroclinic Mode 1 wave.</figcaption></figure>
depends on the frequency. The apparent horizontal phase speeds are higher than those predicted by the non-rotating theory based on Eq. (11). Most of all, their dependency on the frequency limits the utility of the normal mode approach in correctly identifying the constant phase velocity for each baroclinic mode. As far as the multi-dimensional approach of Section 5 is concerned, we assumed that the anomalies of the current were parallel to the wave propagation, whereas a normal current component is now expected in the rotating regime. The limitations of the boundary conditions introduced in Section 5 are actually well evidenced by the relations of polarization in the rotating regime (Gill, 1982, p. 263, Eq. 8.5.3):

$$ \tilde{u} = \frac{k\omega + ilf}{\omega^2 - f^2} \frac{\tilde{p}}{\rho_0}, \quad \tilde{v} = \frac{l\omega - ikf}{\omega^2 - f^2} \frac{\tilde{p}}{\rho_0} \qquad (21) $$

These relations are consistent with, on the one hand, a pressure anomaly of the form $p' = p_0\tilde{p}(z) \exp(i(kx + ly - \omega t))$ and on the other hand, the basic equations given in the previous sections taking now the Coriolis force into account:

$$ \frac{\partial u'}{\partial t} - fv' = - \frac{1}{\rho_0} \frac{\partial p'}{\partial x}, \quad \frac{\partial v'}{\partial t} + fu' = - \frac{1}{\rho_0} \frac{\partial p'}{\partial y} \qquad (22) $$

The physical solution is given by the real part of Eq. (21), that is:

$$ \begin{aligned} \tilde{p}(x,y,t) &= p_0 \cos(kx + ly - \omega t) \\ \tilde{u} &= \frac{p_0}{\rho_0} \left( \frac{k\omega}{\omega^2 - f^2} \cos(kx + ly - \omega t) - \frac{lf}{\omega^2 - f^2} \sin(kx + ly - \omega t) \right) \quad (23) \\ \tilde{v} &= \frac{p_0}{\rho_0} \left( \frac{l\omega}{\omega^2 - f^2} \cos(kx + ly - \omega t) + \frac{kf}{\omega^2 - f^2} \sin(kx + ly - \omega t) \right) \end{aligned} $$

from which it can be deduced that:

$$ \tilde{p}^2 = \rho_0^2 \frac{\omega^2 - f^2}{k^2 + l^2} (\tilde{u}^2 + \tilde{v}^2) - p_0^2 \frac{f^2}{\omega^2 - f^2} \qquad (24) $$

If we introduce the horizontal phase speed $c = \frac{\omega}{\sqrt{k^2 + l^2}}$ and the constant:

$$ \gamma = \left( \frac{\omega f p_0}{\rho_0 c (\omega^2 - f^2)} \right)^2, \qquad (25) $$

Eq. (23) is also equivalent to $\tilde{p}^2 = \rho_0^2 c^2 (1 - \frac{f^2}{\omega^2}) (\tilde{u}^2 + \tilde{v}^2 - \gamma)$. In cases where $\omega^2 \ll N^2$ the horizontal phase speed $c$ is related to the celerity $c_e$, computed from Eq. (11) for the non-rotating regime, according to (Gill, 1982, p. 299):

$$ c = c_e (1 - \frac{f^2}{\omega^2})^{-1/2} \qquad (26) $$

Eq. (24) is finally equivalent to $\tilde{p}^2 = \rho_0^2 c_e^2 (\tilde{u}^2 + \tilde{v}^2 - \gamma)$, in other words we have:

$$ |\tilde{p}| = \rho_0 c_e (\tilde{u}^2 + \tilde{v}^2 - \gamma)^{1/2} \qquad (27) $$

Eq. (27) is relatively similar to boundary conditions proposed in Section 5. An advantage of (27) is that it relies on $c_e$, a celerity that does not depend on frequency (even though $c$ does), but we must, however, note two difficulties.

The first one comes from the appearance of $\gamma$, a constant that is obviously not easy to determine.

Secondly, we note that it is now more complicated to clarify the ambiguity concerning the sign of $\tilde{p}$ at the left hand side of (27). In the non-rotating regime, the sign of the pressure anomaly was simply deduced from the sign of the boundary-normal current. In the rotating case, the Eq. (23) show that this option is not valid anymore, except if the propagation remains normal to the boundary.

We propose to overcome these two problems by defining new variables:

$$ u^* = u' - f \int_{t_0}^{t} v' dt', \qquad (28a) $$

$$ v^* = v' + f \int_{t_0}^{t} u' dt' \qquad (28b) $$

If $(u^*, v^*)$ are introduced into the basic Eq. (22) we obtain:

$$ \frac{\partial u^*}{\partial t} = - \frac{1}{\rho_0} \frac{\partial p'}{\partial x} \qquad (29a) $$

$$ \frac{\partial v^*}{\partial t} = - \frac{1}{\rho_0} \frac{\partial p'}{\partial y} \qquad (29b) $$

The polarization relations (23) become:

$$ \begin{aligned} \tilde{p}(x, y, t) &= p_0 \cos(kx + ly - \omega t) \\ \tilde{u}^* &= \frac{p_0}{\rho_0} \frac{k}{\omega} \cos(kx + ly - \omega t) \qquad (30) \\ \tilde{v}^* &= \frac{p_0}{\rho_0} \frac{1}{\omega} \cos(kx + ly - \omega t) \end{aligned} $$

They can be used to propose a set of open boundary conditions. As far as we consider the open boundary that is perpendicular to the Ox axes, which points outward, a similar reasoning to the one developed in Section 5.1 (namely that the outgoing waves are characterised by $k > 0$) leads to:

$$ p' = +\rho_0 c \operatorname{sign}(u^*) \sqrt{u^2 + v^2} \qquad (31) $$

At the other boundary (Ox pointing inward) the boundary condition is given by $p' = -\rho_0 c \operatorname{sign}(u^*) \sqrt{u^2 + v^2}$. Analogously, boundary conditions at the boundaries perpendicular to Oy are given by $p' = \pm \rho_0 c \operatorname{sign}(v^*) \sqrt{u^2 + v^2}$.

The interest of the boundary condition (31) is its compatibility with multidirectional propagation (as explained in Section 5) but it has the inconvenience of depending on $c$ which according to (26) depends on the a priori unknown frequency. We note however that the latter may be known in advance in some cases (see for instance the simulation of M2 internal waves by Pichon and Correard, 2006).

## 7.2. Determination of the phase speed

The uncertainty in $c$ is the main shortcoming of the boundary condition (31). At this stage, one could be tempted to estimate the phase speed by means of the Orlanski method, as in Eq. (8). However we have seen in Section 3.2 that this method did not provide convincing results. One reason is the possibility of a numerical instability when the denominator of (8) becomes too small. A second reason is that the Orlanski method is not adapted to the even-tuality of reflected waves. Indeed, if we assume that the velocity field $u$ consists of an incident wave $u_1 \cos(kx - \omega t)$ and a reflected wave $u_2 \cos(kx + \omega t)$, then the phase speed deduced from $-\frac{\partial u}{\partial t} = -\frac{\partial v}{\partial x}$ leads to $\frac{\omega}{k} \frac{u_1 - u_2}{u_1 + u_2}$, that is not equal to $\omega/k$, unless the reflected wave is zero (i.e. $u_2 = 0$). In contrast, the general principle of the Orlanski method (i.e. using the inner solution to identify some parameters of the wave propagation) does seem suited to the identification of the wave frequency. Indeed, considering again the case of the combination of an incident wave and a reflected wave, $u = u_1 \cos(kx - \omega t) + u_2 \cos(kx + \omega t)$, we see that $-\frac{\partial^2 u}{\partial u^2} = -\frac{\partial^2 v}{\partial v^2}$ equals $\omega^2$ whatever $u_1$ or $u_2$. We note that this calculation can become inaccurate when $u$ approaches zero, but because theoretically $u$ and $v$ cannot both vanish at the same time, the frequency can be alternatively obtained from $-\frac{\partial^2 v}{\partial v^2}$. In practice, to eliminate the inaccuracy when either $u$ or $v$ approaches zero we propose computing $\omega$ using a weighted average of these two options, namely $\omega^2 = -(W_1 \frac{\partial^2 u}{\partial u^2} + W_2 \frac{\partial^2 v}{\partial v^2}) (W_1 + W_2)^{-1}$, where the weights, $W_1$ and $W_2$, are respectively $u^2$ and $v^2$. We obtain:

$$ \omega^2 = - \frac{u \partial^2 u / \partial t^2 + v \partial^2 v / \partial t^2}{u^2 + v^2} \qquad (32) $$
The resulting frequency can now be used to compute the phase speed from Eq. (26).

### 7.3. Discretization

We note that the boundary condition (31) is a generalization of the multi-directional condition given in Section 5. They actually become equivalent if the effects of rotation are ignored. In this particular case, according to (28), we simply have $u^* = u'$ and $v^* = v'$. The options of discretization proposed in Sections 2.2, 3.4 and 5.2 can be applied, including the stabilising procedure in case of nearly tangential propagation as stated in (18) and (19). When effects of rotation are not neglected, the time integral terms appearing at the RHS of (28), $S_v = -f \int_{t0}^t v'dt'$ and $S_u = f \int_{t0}^t u'dt'$, have to be taken into consideration. The latter can be computed using a first order scheme such as:

$$S_u^t = S_u^{t-\Delta t} + f u^{t\Delta t} \quad (33a)$$

$$S_v^t = S_v^{t-\Delta t} - f v^{t\Delta t} \quad (33b)$$

At this stage, we have tested two options. The first one is a straightforward discretization of (28) and (33), namely:

$$u^* = u' + S_v \quad (34a)$$

$$v^* = v' + S_u \quad (34b)$$

In the second option (note that we are still considering the case of the boundary perpendicular to $Ox$), we keep (34a) but we constrain $v^*$ to satisfy the tangential momentum Eq. (29b). Indeed, because of truncation errors, (34b) and (33a) may be numerically inconsistent with (29b). This requirement can be recovered by replacing (34b) by $v^{*t} = -\frac{1}{\rho_0} \int_{t0}^t \frac{\partial p'}{\partial y} dt' + v^{*t0}$ discretized according to:

$$v^{*t} = v^{*t-\Delta t} - \frac{1}{\rho_0} \frac{\partial p'}{\partial y} \Delta t \quad (35)$$

where the tangential pressure gradient, $\frac{\partial p'}{\partial y}$, has been computed on inner grid points contiguous to the 'open' boundary, i.e. at $i = M - 1$ and at $i = 2$ on Fig. 1.

### 7.4. Numerical tests

The boundary condition (31) is tested in the case of multi-directional wave propagation at latitude $45^\circ$N. The basic parameters of these simulations are summarized in Table 3.

Waves corresponding to the baroclinic Mode 1 are excited with a period of 12 h, by a generating force acting in the centre of the domain. This generating force, introduced in the momentum equations of the model, is derived from a potential function, $\Psi$, defined as:

$$\Psi = 0.5 \cos \left( \frac{\pi z}{h} \right) \sin(\omega t) e^{-\frac{d}{2\lambda}}$$

where $d$ is the distance from the central point of the domain and $\Delta x$ the horizontal size of a grid cell. We note that $e^{-\frac{d}{2\lambda}}$, and consequently $\Psi$, becomes negligible a few grid points away from the centre of the domain. As in the previous sections, two domains are

<table><caption>**Table 3**
Model parameters for 3D $Oxyz$ simulations in the rotating case.</caption><tr><td>Δx × Δy</td><td>3000 × 3000 m</td></tr><tr><td>Horizontal extent</td><td>1500 × 1500 km(large domain) 300 × 300 km(small domain)</td></tr><tr><td>Duration of simulation</td><td>4 days<sup>a</sup></td></tr><tr><td>Coriolis parameter</td><td>f = 1.028 × 10<sup>-4</sup> s<sup>-1</sup></td></tr><tr><td>Other parameters</td><td>As in Table 1</td></tr></table>

used: a large domain which provides a reference solution, and a small domain.

Fig. 16a and b show the reference solution after 4 days, which is the time required for waves to reach the open boundaries of the large domain. Fig. 16c shows the corresponding solution computed with the small domain and PRM conditions including the rotating effect. Fig. 16b and c show a good agreement and, as a matter of fact, the kinetic energy of the current reflected on the open boundaries of the small domain, computed as in the previous sections, is particularly small (about $2 \times 10^{-3} \times E_0$). We note that Fig. 16c has been obtained with the second option, using (35), as described in Section 7.3. The solution obtained with the first option, using (34b), developed unrealistic small-scale patterns near the boundaries and is consequently not presented. Fig. 16d was obtained with the non-rotating boundary condition described in Section 5. Its relative disagreement with Fig. 16b shows that the omission of the rotating effect significantly damages the radiation properties of the PRM condition. The Sommerfeld condition is neither convincing even if, as in the case of Fig. 16e, the phase speed used in (7) is set to the phase speed of the baroclinic Mode 1 at the 12 h period, i.e. $c = c_e(1 - f^2/\omega^2)^{-1/2} \approx 3.1$ ms<sup>-1</sup>.

## 8. Consistency with boundary conditions for temperature and salinity

In this study, the boundary conditions for pressure have been treated independently from the tracer boundary conditions. We have shown, within the limits of the numerical tests presented in this paper, that this option seems relevant. Consistency between the tracers and the pressure conditions can, however, easily be ensured. In other words, it is possible to ensure that the pressure computed by the boundary conditions, and the pressure that would be derived from the boundary temperature and boundary salinity, through the equation of state, are identical. This can be done by using the hydrostatic equilibrium to relate the pressure and density anomalies, $\rho' = -\frac{1}{g} \bar{p}(x, y, t) \frac{\partial \bar{p}}{\partial z}$, and then the two additional constraints converting the latter into temperature and salinity perturbations. These two constraints may consist of, on one hand, a linear state equation of the type $\rho' = \alpha T' + \beta S'$ (where $T'$ and $S'$ are temperature and salinity anomalies) and on the other hand, on an assumption that the tracer variations are due to the vertical displacements, $h_w$, of the reference field, namely $(T', S') = h_w \left( \frac{\partial T_{\text{ref}}}{\partial z}, \frac{\partial S_{\text{ref}}}{\partial z} \right)$. If the density anomaly is known, the vertical displacements are given by $h_w = \frac{-\rho'}{\alpha \frac{\partial T_{\text{ref}}}{\partial z} + \beta \frac{\partial S_{\text{ref}}}{\partial z}}$, from which the temperature and salinity anomalies are easily deduced.

## 9. Conclusions

We have presented an open boundary condition scheme based on polarization relations for internal waves. In simple cases of non-rotating regime and boundary-normal propagation, the so called PRM scheme can be seen as a three-dimensional version of the well-known Flather barotropic conditions. The discretization, based on considerations of numerical stability and energy conservation, relates the open boundary pressure to inner domain currents. An additional boundary condition is required in order to compute advective and diffusive momentum fluxes. The latter has little influence on radiation properties as long as the current is weak enough. Numerical tests, reproducing a classical baroclinic Mode 1 propagation of internal waves in a continuously stratified ocean, show that the PRM radiation properties are similar to those of usual Sommerfeld conditions if not better, particularly when the latter estimate phase speed according to the Orlanski's method. The PRM conditions also have the advantage, over the wave

<footnote>a Time required for waves to reach the open boundaries of the large domain.</footnote>
<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /><figcaption>**Fig. 16.** Surface $Ox$ velocity component after 4 days. (a) Reference solution over the large domain. The black square indicates the open boundaries of the small domain. (b) Same as the previous figure but enlarged over the small domain area in order to facilitate the comparisons with the next figures. (c) Small domain with PRM conditions of Section 7, (d) small domain with PRM conditions of Section 5, (e) small domain with Sommerfeld conditions.</figcaption></figure>

equation type conditions, of being easily compatible with multi-directional wave propagation and of providing an efficient forcing scheme when incoming waves are introduced as a reference state of the boundary conditions.

In the case of the non-rotating regime, the dispersive nature of internal waves is well treated by the normal mode approach that

permits each mode to radiate separately with a proper phase speed. The PRM condition is theoretically limited to waves of small amplitude but our numerical tests showed that it performs reasonably well when amplitudes are comparable to those observed in the Bay of Biscay, a region well known for the strength of internal waves (Gerkema et al., 2004). This potentially weak point can
anyway be reduced by ensuring, when possible, that open boundaries are not too close to the area of wave generation (Chini and Leibovich, 2003). The bottom slope induces a coupling effect between modes that possibly limits the scope of the modal approach. However, our numerical tests show that radiation properties remain correct for a wide range of bathymetric situations and we note that strong bathymetry variations are anyway expected to reflect internal waves.

The modal approach leads to a Sturm-Liouville type problem. Its solution is possibly prohibitive. As far as our numerical tests are concerned, its computing cost could be comparable to that of the simulation itself. But on the other hand the variations of the background density field are in many cases slow enough so that it is not necessary to recalculate the vertical modes and their corresponding phase speed at each time step and consequently the computing cost of the modal approach can be significantly reduced, until being small compare to that of the full model solution.

The case of the rotating regime is more complicated, but it is possible to retain the formulation derived from the equations of the non-rotating regime, provided that a set of pseudo-velocity variables is suitably defined and used instead of the current field. The uncertainty on the wave phase speed is obviously a weak point of this scheme but it can be overcome if the frequency of waves can be identified.

The numerical experiments presented in this study are encouraging. However, the relevance of the PRM conditions in more general cases, e.g. bores, fronts or large scale flows such as Rossby waves, still needs to be demonstrated.

## Acknowledgements

This study was funded by the European MFSTEP Project (EU Contract EVK3-CT-2002-00075), the European INSEA Project (SST4-CT-2005-012336), the "LEFE-IDAO-2006" and "OPA Symphonique" French projects. The authors thank Jérôme Chanut (Mercator-Océan) for helpful discussions, the reviewers, who provided many suggestions and comments that improved the manuscript, Cyril Nguyen and the Laboratoire d'Aérologie computer team, Serge Prieur, Laurent Cabanas, Jérémy Leclercq, Didier Gazen, and Juan Escobar for their support.

## References

Arakawa, A., Lamb, V.R., 1977. Computational design of the basic dynamical processes of the UCLA general circulation model. Methods in Computational Physics 17, 174-267.

Barnier, B., Marchesiello, P., De Miranda, A.P., Molines, J.-M., Coulibaly, M., 1998. A sigma-coordinate primitive equation model for studying the circulation in the South Atlantic. Part I: Model configuration with error estimates. Deep Sea Research Part I 45, 543-572.

Bennett, A.F., Chua, B.S., 1994. Open-ocean modelling as an inverse problem: the primitive equation. Monthly Weather Review 122, 1326-1336.

Blayo, E., Debreu, L., 2005. Revisiting open boundary conditions from the point of view of characteristic variables. Ocean Modelling 9, 231-252.

Blumberg, A.F., Kantha, L.H., 1985. Open boundary conditions for circulation models. Journal of Hydraulic Engineering 11, 237-255.

Blumberg, A.F., Mellor, G.L., 1987. A description of a three-dimensional coastal circulation model, three-dimensional coastal ocean models. In: Heaps, N. (Ed.), Coastal Estuarine Science, vol. 4. AGU, Washington, DC, pp. 1-16.

Bougeault, P., 1983. A non-reflective upper boundary condition for limited-height hydrostatic models. Monthly Weather Review 111, 420-429.

Carter, G.S., Merrifield, M.A., 2007. Open boundary conditions for regional tidal simulations. Ocean Modelling 18, 194-209.

Chini, G.P., Leibovich, S., 2003. An analysis of the Klemp and Durran radiation boundary condition as applied to dissipative internal waves. Journal of Physical Oceanography 33, 2394-2407.

Engquist, B., Mayda, A., 1977. Absorbing boundary conditions for numerical solution of waves. Mathematics of Computation 31, 629-651.

Flather, R.A., 1976. A tidal model of the northwest european continental shelf. Memoires. Société Royale des Sciences de Liège 6 (10), 141-164.

Gaspar, P., Gregoris, Y., Lefevre, J.M., 1990. A simple eddy kinetic energy model for simulations of the oceanic vertical mixing: tests at station Papa and long-term upper ocean study site. Journal of Geophysical Research 95, 16,179-16,193.

Gerkema, T., Lam, F.P.A., Maas, L.R.M., 2004. Internal tides in the bay of Biscay: conversion rates and seasonal effects. Deep-Sea Research Part II 51, 2995-3008.

Gill, A.E., 1982. Atmosphere-Ocean Dynamics. Academic Press, 662 pp.

Higdon, R.L., 1994. Radiation boundary conditions for dispersive waves. SIAM Journal on Numerical Analysis 31, 64-100.

Jensen, T.G., 1993. Equatorial variability and resonance in a wind-driven Indian ocean model. Journal of Geophysical Research 98, 22533-22552.

Jensen, T.G., 1998. Open boundary conditions in stratified ocean models. Journal of Marine Systems 16, 297-322.

Johns, B., Sinha, P.C., Dube, S.K., Mohanty, U.C., Rao, A.D., 1983. Simulation of storm surges using a three-dimensional numerical model: an application to the 1977 Andhra Cyclone. Quarterly Journal of the Royal Meteorological Society 109, 211-224.

Klemp, J.B., Durran, D.R., 1983. An upper boundary condition permitting internal gravity waves radiation in numerical mesoscale models. Monthly Weather Review 111, 430-444.

Kourafalou, V.H., Oey, L.-Y., Wang, J.D., Lee, T.N., 1996. The fate of river discharge on the continental shelf, 1. Modelling the river plume and the inner shelf current. Journal of Geophysical Research 101, 3415-3434.

Kourafalou, V.H., Tsiaras, K.P., 2006. A nested circulation model for the North Aegean Sea. Ocean Science Discussions 3, 343-372.

Lavelle, J.W., Thacker, W.C., 2008. A pretty good sponge: dealing with open boundaries in limited-area ocean models. Ocean Modelling 20, 270-292.

Marchesiello, P., McWilliams, J.C., Shchepetkin, A., 2001. Open boundary conditions for long-term integration of regional oceanic models. Ocean Modelling 3, 1-20.

Martinsen, E.A., Engedahl, H., 1987. Implementation and testing of lateral boundary schemes as an open boundary condition in a barotropic ocean model. Coastal Engineering 11, 603-627.

Marsaleix, P., Auclair, F., Estournel, C., 2006. Considerations on open boundary conditions for regional and coastal ocean models. Journal of Atmospheric and Oceanic Technology 23, 1604-1613.

Marsaleix, P., Auclair, F., Floor, J.W., Herrmann, M.J., Estournel, C., Pairaud, I., Ulses, C., 2008. Energy conservation issues in sigma-coordinate free-surface ocean models. Ocean Modelling 20, 61-89.

Miyakoda, K., Rosati, A., 1977. One-way nested grid models: The interface conditions and the numerical accuracy. Monthly Weather Review 105, 1092-1107.

Morel, Y., Baraille, R., Pichon, A., 2008. Time splitting and linear stability of the slow part of the barotropic component. Ocean Modelling 23, 73-81.

Nycander, J., Döös, K., 2003. Open boundary conditions for barotropic waves. Journal of Geophysical Research 108 (C5), 3168. doi:10.1029/2002JC001529.

Nycander, J., McC. Hogg A., Leela, M. Frankcombe, L.M., 2008. Open boundary conditions for nonlinear channel flow. Ocean Modelling 24, 108-121.

Oey, L.Y., Chen, P., 1992. A model simulation of circulation in the northeast atlantic shelves ans seas. Journal of Geophysical Research 97, 20,087-20,115.

Orlanski, I., 1976. A simple boundary condition for unbounded hyperbolic flows. Journal of Computational Physics 21, 251-269.

Palma, E.D., Matano, R.P., 1998. On the implementation of passive open boundary conditions for a general circulation model: the barotropic mode. Journal of Geophysical Research 103, 1319-1341.

Pichon, A., Correard, S., 2006. Internal tides modelling in the bay of Biscay. Comparisons with observations. Scientia Marina 70S1, 65-88.

Pinardi, N., Allen, I., Demirov, E., de Mey, P., Korres, G., Lascaratos, A., Le Traon, P.-Y., Maillard, C., Manzella, G., Tziavos, C., 2003. The Mediterranean Ocean forecasting system: first phase of implementation (1998-2001). Annales Geophysicae 21, 3-20.

Raymond, W.H., Kuo, H.L., 1984. A radiation boundary condition for multidimensional flows. Quarterly Journal of the Royal Meteorological Society 110, 535-551.

Shchepetkin, A.F., McWilliams, J.C., 2005. Regional Ocean Model System: a split- explicit ocean model with a free-surface and topography-following vertical coordinate. Ocean Modelling 9, 347-404.

Shulman, I., 1997. local data assimilation in specification of open boundary conditions. Journal of Atmospheric and Oceanic Technology 14, 1409-1419.

Shulman, I., Lewis, J.K., Blumberg, A.F., Nicholas Kim, B., 1998. Optimized boundary conditions and data assimilation with application to the M2 tide in the Yellow Sea. Journal of Atmospheric and Oceanic Technology 15, 1066-1071.

Sofianos, S.S., Skliris, N., Mantziafou, A., Lascaratos, A., Zodiatis, G., Lardner, R., Hayes, D., Georgiou, G., 2006. Nesting operational forecasting models in the Eastern Mediterranean: active and slave mode. Ocean Science Discussions 3, 1225-1254.

Trefethen, L.N., 1982. Group velocity in finite difference schemes. SIAM Review 24, 113-132.

Treguier, A.M., Barnier, B., de Miranda, A.P., Molines, J.M., Grima, N., Imbard, M., Messenger, C., Reynaud, T., Michel, S., 2001. An eddy-permitting model of the Atlantic circulation: evaluating open boundary conditions. Journal of Geophysical Research 106, 22115-22130.