About Interface Conditions for Coupling Hydrostatic and Nonhydrostatic Navier-Stokes Flows

. In this work we are interested in the search of interface conditions to couple hydrostatic and nonhydrostatic ocean models. To this aim, we consider simpliﬁed systems and use a time discretization to handle linear equations. We recall the links between the two models (with the particular role of the aspect ratio δ = H/L (cid:28) 1) and introduce an iterative method based on the Schwarz algorithm (widely used in domain decomposition methods). The convergence of this method depends strongly on the choice of interface conditions: this is why we look for exact absorbing conditions and their approximations in order to provide tractable and eﬃcient coupling algorithms.

(Communicated by the associate editor name) Abstract. In this work we are interested in the search of interface conditions to couple hydrostatic and nonhydrostatic ocean models. To this aim, we consider simplified systems and use a time discretization to handle linear equations. We recall the links between the two models (with the particular role of the aspect ratio δ = H/L 1) and introduce an iterative method based on the Schwarz algorithm (widely used in domain decomposition methods). The convergence of this method depends strongly on the choice of interface conditions: this is why we look for exact absorbing conditions and their approximations in order to provide tractable and efficient coupling algorithms.

1.
Introduction. Almost all ocean circulation models, either at large scale or at regional scale, are based on the so called primitive equations, which make use of the hydrostatic approximation. From a dynamical point of view, this approximation is justified by the fact that oceanic flows are generally characterized by large differences between horizontal and vertical scales and by a strong vertical stratification that avoids vertical mixing. From a computational point of view, the hydrostatic approximation decreases the computational burden by one order of magnitude w.r.t. solving the nonhydrostatic equations. However continuous improvement in numerical modeling and in computing resources leads to more and more sophisticated ocean modeling systems, which aims at representing the full ocean physics. A natural idea is thus to build systems that couple local nonhydrostatic models to larger scale hydrostatic ones. Since the pioneering work of Marshall et al. [13,14], that aimed at locally adapting the physics of a unique model rather than at coupling two different models, some few works addressed this question. Let us cite for instance the effort by Fringer et al. [6,7] to couple the Navier-Stokes SUNTANS model to the regional hydrostatic oceanic model ROMS, or the work by Gallacher et al. [8] nesting a nonhydrostatic model in a hydrostatic one. Let us also mention the work of Audusse et al. [1] that tackles domain decomposition issues using Schwarz methods for the hydrostatic equations. However the coupling procedure in these works is essentially ad hoc and does not rely on a sound mathematical basis. This comes from the fact that such a coupling is quite delicate from a mathematical point of view, due to the different nature of hydrostatic and nonhydrostatic Navier-Stokes equations (where the vertical velocity is either a diagnostic or a prognostic variable). To the best of our knowledge, there is still almost no work addressing this problem. In this paper we propose to couple such systems through a Schwarz iterative algorithm and we seek relevant interface conditions that would both lead to a well posed coupled problem and minimize the computational cost.
2. Navier-Stokes equations. We consider in the following a shallow fluid domain Ω, meaning that its horizontal length scale L is much larger than its vertical one H. Thus the aspect ratio δ = H/L satisfies δ 1.

2.1.
Hydrostatic and nonhydrostatic systems. The usual 3D Navier-Stokes system reads: where v = (u, v, w) T is the velocity vector, p is the pressure, ∇ and div are the gradient and divergence operators, ∆ h = ∂ 2 x + ∂ 2 y is the horizontal Laplacian operator, and ν h , ν v are the horizontal and vertical diffusion coefficients (with ν v = δ 2 ν h ). This system must of course be complemented with initial and boundary conditions. In order to write this system in a dimensionless form, let us introduce the following dimensionless variables (denoted by a prime symbol): where U is a typical scale for the horizontal velocity. Making these changes in (1) and (2) leads to the dimensionless system (see e.g. [2,12] for further details): where prime symbols have been dropped for the sake of readability and Re = LU/ν h is the Reynolds number. The hydrostatic approximation consists in considering the dominant term only in (5) since δ 1, which leads to the hydrostatic Navier-Stokes equations (HNS): div v = 0 (10)

COUPLING HYDROSTATIC AND NONHYDROSTATIC NAVIER-STOKES FLOWS 3
Allowing for a variable density and considering the gravity force leads to the socalled primitive equations, which are widely used for the simulation of large scale ocean circulation (see e.g. [11,18,19] and the numerous references herein). Naturally, before coupling such systems of equations, one should check that each of the two systems is well-posed. We know that it is not (yet) the case for the 3D Navier-Stokes equations, but as far as the primitive equations are concerned, we can rely on the recent works of [4] or [10] (see also [18]).

2.2.
Transmission conditions. Since we are interested by coupling the nonhydrostatic Navier-Stokes (NS) system (3)-(6) with its hydrostatic (HNS) approximation (7)-(10), let us compute the corresponding so-called natural transmission conditions, i.e. the physical quantities that are naturally conserved through any fluid boundary 1 . These quantities are obtained by writing the variational formulation of each system and correspond to the argument of the boundary integrals. They will be used in Section 4 to provide relevant interface conditions for the coupling problem.
Let denote by n the outward normal vector to any fluid boundary. Simple calculations (see also e.g. [17]) show that those naturally transmitted quantities for (NS) are equal to where v δ = (u, v, δ 2 w) T . Similar calculations conducted with the hydrostatic system (7)-(10) leads to the transmitted This corresponds to making δ equal to 0 in (11), which is obviously consistent with the fact that (HNS) is obtained as a limit of (NS) for δ → 0 (see [2]). Moreover, considering that the fluid boundary is vertical (which is actually almost always the case), n is horizontal, which means that the third component of the transmitted quantities is equal to 0 for (HNS). In other words, (HNS) transmits one scalar quantity less than (NS), and thus it requires one boundary condition less than (NS). We will be back later on this important point.
3. The Schwarz coupling method. As indicated in Section 1, our aim is to couple two existing (NS) and (HNS) numerical models in a way that requires as few modifications as possible to the models. That is why the Schwarz algorithm, developed in the context of domain decomposition methods, is appropriate (see [9] for a historical review of this method).
3.1. Schwarz algorithm. Let us first introduce this method in a very general way. We are interested in solving a coupled problem of the form: where Ω 1 and Ω 2 are two adjacent open domains, which interface is Γ =Ω 1 ∩Ω 2 . The external boundary of Ω i (i = 1, 2) is denoted ∂Ω ext i = ∂Ω i \ Γ, and L i and B i denote partial differential operators representing the model equations and the boundary conditions respectively. This coupled problem must be complemented with compatibility conditions of the form where F i (i = 1, 2) also denote partial differential operators. Note that since we are exclusively interested in this paper in coupling two different sets of equations, we consider that Ω 1 and Ω 2 do not overlap (otherwise the solution in Ω 1 ∩ Ω 2 should exactly satisfy both systems, which generally does not make sense). Note also that time evolution can be easily included in this general formulation, by simply adding initial conditions and replacing the space domain Ω i by a space-time domain In order to solve for (12)-(13), the Schwarz algorithm reads: where m is the (Schwarz) iteration index. Interface operators C i and D i (i = 1, 2) are not specified at this stage. They must of course lead to a well posed problem on each domain, and must also ensure that the algorithm converges. The link between the converged solutions u 1 and u 2 of (14) at the interface Γ obviously reads: while we know from the original coupled problem that it also reads (13). This means that (15) must imply (13).
Remark 1. If the compatibility conditions (13) are not sufficient to ensure wellposedness of the coupled system (which will be our case for the hydrostatic / nonhydrostatic coupling, see below), it is clear that they cannot be equivalent to conditions (15) that must provide a well-posed (and convergent) algorithm. Consequently conditions (15) imply (13) but are not always equivalent.
Note that (14) is the so-called multiplicative (or sequential) version of the algorithm, for which the computation on Ω 1 must be performed before the computation on Ω 2 . Replacing D 1 u m 1 by D 1 u m−1 1 allows to compute both systems in parallel, and is called the additive (or parallel) version.

Convergence issues.
In order to study the convergence of the Schwarz algorithm, let us introduce the errors e m i = u m i − u i . Assuming that all operators are linear (or approximating non linear operators by their linearized version) and that (15) is satisfied, we can derive the algorithm satisfied by the errors: Under this form it is clear that the convergence speed of the algorithm relies on the choice of the interface operators C i and D i (i = 1, 2). We can even notice that if there exists some operator C 2 such that C 2 e 1 2 = 0 and/or some operator D 1 such that D 1 e 1 1 = 0, then the algorithm will converge in only two iterations. Such operators are called transparent or perfectly absorbing interface operators ( [5,15,16]). Note however that they are generally non local in time nor in space, which means that they cannot be used in actual computations. Moreover there is no reason why absorbing operators corresponding to (16) should satisfy the compatibility conditions (13). Therefore our goal in the following will be to look for interface operators C i and D i that both satisfy the compatibility conditions (13) (in order to get the solution of the original coupled problem) and are local approximations of the perfectly absorbing interface operators (in order to ensure a fast convergence of the coupling algorithm).

Simplified systems of equations.
For the sake of simplicity, we will restrict in the following to a 2D x-z version of the NS system (3) where ∆t is the time step and where the superscript n relates to the approximate solution at the n th time step. Thus one has to solve a Stokes problem (corresponding to the L 1 operator in section 3) at each time step: Its hydrostatic counterpart in Ω 2 (corresponding to the L 2 operator in section 3) reads: In summary we want to couple the two systems (20)-(22) and (23)-(25), complemented of course with boundary conditions on ∂Ω ext i . These two sets of equations correspond to the coupled system written in generic form (12) in section 3. Moreover their compatibility conditions, corresponding to (13) in section 3, can be directly deduced from (11). Since we couple the hydrostatic and nonhydrostatic models, the interface must obviously stand in a hydrostatic region (where both models are relevant). Consequently the compatibility conditions read: We have now to determine interface conditions on Γ for the Schwarz coupling algorithm, i.e. to choose operators corresponding to C i and D i in (14) that will make the corresponding boundary value problem well-posed. At convergence, these conditions must imply the compatibility conditions (26) and (27) (and possibly more, see Remark 1 above). Following the method explained at the end of Section 3.2, we will first determine absorbing conditions for the two systems (20)-(22) and (23)-(25).

4.2.
Expression of the errors in the Schwarz algorithm. The first step consists in determining approximate analytical expressions for the errors e m 1 and e m 2 . Regarding the Stokes system (20)-(22), a very similar calculation was already performed in our paper [3], except that it was not conducted in a shallow domain (i.e. there were no δ in the equations). The adaptation of this previous work being quite straightforward, we will only give here a rapid sketch of the derivation, the reader being referred to [3] for further details.
The same kind of computation in the hydrostatic case (x > 0) leads to: so that we finally have: 4.3. Exact absorbing conditions. We now look for exact absorbing conditions on Γ, as defined in section 3.2. As we did in [3], we search them as perturbations of the natural transmission conditions: for the nonhydrostatic equations, and for the hydrostatic equations, where the operators S and S H are to be chosen appropriately. Note once again the fact that 2D NS leads to two scalar equations while 2D HNS leads to only one. For the nonhydrostatic side Ω 1 , we obtain from (39) after a Fourier transform: whereŜ is defined byŜv = Sv. In order to obtain the exact absorbing boundary conditions, the operator S must be such that the right-hand-side of (41) is zero. Replacing the errorê 2 by its expression (38) and making x = 0, one can easily check that this is achieved witĥ i.e.Ŝ = λ Re Id, where we recall that λ = Re ∆t + k 2 .
For the hydrostatic side Ω 2 , we proceed similarly starting from (40): In order to cancel the right-hand-side of (43), we use the expression ofê 1 given by (34). ThusŜ H = (a b) is such that: Finally, the exact absorbing boundary conditions are provided by equations (39) and (40) with operators S and S H defined as the inverse Fourier transforms of Remark 2. In the hydrostatic limit δ −→ 0, one can check that so that in this limit,Ŝ H = −λ/Re (1 0) is coherent with the definition ofŜ = λ Re Id.

4.4.
Coupling algorithm. We summarize hereafter the iterative algorithm to be run until convergence: the operators S and S H being defined by their Fourier transform in equations (46). Unfortunately (and this is a traditional issue with absorbing boundary conditions), these operators are non-local (see the definition of λ) and we must approximate them by local operators to get a tractable algorithm. As we did in [3], we propose to introduce the parameter ε = k 2 ∆t/Re so that λ = Re ∆t Assuming that ε is small we can introduce two approximations for λ: Order 0: Order 1: These approximations provide local approximations for the operator S (and also for S H since the aspect ratio δ is small). For λ = λ 0 , S is a simple Dirichlet operator, and for λ = λ 1 it reads Re ∆t Id − ∆t 2 Re ∂ zz .
These computations provide ansatz for the search of accurate interface operators to be implemented in the above coupling algorithm. These ansatz can be used to provide the type of boundary operators, with scalar coefficients that will need to be optimized, as we did in [3], keeping in mind that the chosen operators must absolutely ensure that the converged solution satisfies the compatibility conditions (26) and (27). Interface boundary conditions could be search with with coefficients α, α H and β to be optimized. As a first guess (see discussion above), one could set β = 0 and α H = α.
This work can be extended in several directions: • accounting for the free surface in geophysical flows is a challenging question, • adding the third dimension along the y axis would not change the modeling principles above, but of course this would complicate the calculations (in particular an interface can seldom be defined as x = cst), • numerical validation for coupled systems is much more challenging than with traditional DDM problems. Indeed, coupled problems generally do not have reference solutions. In the present case, the coupled reference solution is not expected to be the fully non-hydrostatic solution, although not too far from it. Testing and validating our coupling algorithms is thus a demanding work, that will require a collaboration with physicists.
Finally a future objective is now to address global-in-time Schwarz methods in order to improve the cost/accuracy ratio of this coupling algorithm. This is let to future works.