E Audusse

B. MERLET P Dreyfuss

OPTIMIZED SCHWARZ WAVEFORM RELAXATION FOR PRIMITIVE EQUATIONS OF THE OCEAN

Keywords: Domain Decomposition, Schwarz Waveform Relaxation Algorithm, Fluid Mechanics, Primitive Equations, Finite Volume Methods AMS subject classifications. 65M55, 76D05, 76M12

In this article we are interested in the derivation of efficient domain decomposition methods for the viscous primitive equations of the ocean. We consider the rotating 3d incompressible hydrostatic Navier-Stokes equations with free surface. Performing an asymptotic analysis of the system with respect to the Rossby number, we compute an approximated Dirichlet to Neumann operator and build an optimized Schwarz waveform relaxation algorithm. We establish the wellposedness of this algorithm and present some numerical results to illustrate the method.

Introduction.

A precise knowledge of ocean parameters (velocity, temperature...) is an essential tool to obtain climate and meteorological previsions. This task is nowadays of major importance and the need of global or regional simulations of the evolution of the ocean is strong. Moreover the large size of global simulations and the interaction between global and regional models require the introduction of efficient domain decomposition methods.

The evolution of the ocean is commonly modelized by the use of the viscous primitive equations. This system is deduced from the full three dimensional incompressible Navier-Stokes equations with free surface with the use of the hydrostatic approximation and of the Boussinesq hypothesis. It is implemented in all the major softwares that are concerned with global or/and regional simulations of ocean and/or atmosphere (we refer for example to NEMO [START_REF] Madec | OPA 8.1 Ocean General Circulation Model reference manual, Note du Pole de modélisation[END_REF], MOM [START_REF] Pacanowski | Manual[END_REF] or HYCOM for global models and ROMS [START_REF] Arango | ROMS : A Regional Ocean Modeling System[END_REF] or MARS for regional models). The primitive equations have been studied for twenty years and important theoretical results are now available [START_REF] Lions | New formulations of the primitive equations of the atmosphere and applications[END_REF][START_REF] Temam | Some mathematical problems in geophysical fluid dynamics[END_REF][START_REF] Cao | Global well-posedness of the three-dimensional viscous primitive equations of large scale ocean and atmosphere dynamics[END_REF]. The numerical treatment of this system has been also strongly investigated [START_REF] Temam | Computational methods for the oceans and the atmosphere[END_REF]. But the key point here is to simulate global circulation on the earth for long time and/or with small space discretization. This type of computations can not be performed on a single computer in realistic CPU time and need to be parallelized. The problem is then to allow the different subdomains to interact in an efficient way. Another type of applications that are commonly investigated in the oceaonographic and/or meteorological community is to couple global and regional models in order to obtain precise regional previsions. The problem is also to construct an efficient interaction between the two models. In [START_REF] Cailleau | Comparison of different numerical methods used to handle the open boundary of a regional ocean circulation model of the Bay of Biscay[END_REF] the authors exhibit that most of the existing algorithms are not able to compute this kind of problem in an efficient way. We propose in this article to investigate these still open questions in the context of a quite recent performing domain decomposition method : the Schwarz waveform relaxation type algorithms.

The development of domain decomposition techniques have known a great development for the last decades and our purpose is not to make an exhaustive presentation of these methods. We refer the reader to [START_REF] Quarteroni | Domain Decomposition Methods for PDEs[END_REF][START_REF] Toselli | Domain decomposition methods -Algorithms and theory[END_REF] for a general presentation and we restrict ourselves to the description of Schwarz waveform relaxation method. It is a relatively new domain decomposition technique. It has been developed for the last decade and has been successfully applied to different types of equations. This type of algorithms is the result of the interaction between classical Schwarz domain decomposition techniques and waveform relaxation algorithms. Its great interest is to be explicitly designed for evolution equations and to allow different strategies for the space time discretization in each subdomain. Moreover we can even consider different models in each subdomain without modifying the architecture of the interaction. The heart of the classical Schwarz method is to solve the problem on the whole domain thanks to an iterative procedure where a problem is solved on each subdomain by the use of boundary conditions that contain the information coming from the neighboring subdomains. It comes from the early work of Schwarz [START_REF] Schwarz | Über einen Grenzübergang durch alternierendes Verfahren[END_REF] where this idea was introduced to prove the well-posedness of a Poisson problem in some nontrivial domains. This method is designed for stationary problems and presents two main drawbacks : it needs an overlapping between subdomains and it converges slowly [START_REF] Lions | On the Schwarz alternating method[END_REF]. In the last decade, some works have been devoted to cure these disagreements [START_REF] Lions | On the Schwarz alternating method II: a variant for nonoverlapping subdomains[END_REF]. We refer to [START_REF] Gander | Optimized Schwarz methods[END_REF] for a complete presentation. The extension to time evolution problems was performed at the end of the nineties by Gander [START_REF] Gander | Overlapping Schwarz for parabolic problems[END_REF][START_REF] Gander | A waveform relaxation algorithm with overlapping splitting for reaction diffusion equations[END_REF] and Giladi & Keller [START_REF] Giladi | Space time domain decomposition for parabolic problems[END_REF] and was denoted Schwarz waveform relaxation algorithms. The authors mixed the classical Schwarz approach with waveform relaxation techniques developed in the context of the solutions of large system of ordinary differential equations [START_REF] Lelarasmee | The waveform relaxation method for time-domain analysis of large scale integrated circuits[END_REF][START_REF] Jeltsch | Waveform relaxation with overlapping splittings[END_REF]. The exchanged quantities were of Dirichlet type. Optimized Schwarz waveform relaxation methods were developed with the introduction of more sophisticated information to compute the interaction between the subdomains. These optimized algorithms were based on previous works [START_REF] Engquist | Absorbing boundary conditions for the numerical simulation of waves[END_REF][START_REF] Halpern | Artificial boundary conditions for the advection diffusion equations[END_REF][START_REF] Halpern | Artificial boundary conditions for incompletely parabolic perturbations of hyperbolic systems[END_REF] about the derivation of absorbing boundary conditions respectively for hyperbolic, elliptic and incompletely parabolic equations. The same ideas were used to derive efficient transmission conditions between the subdomains : since the exact transparent conditions can not be implemented in general (it may lead to non-local pseudo-differential operators), the derivation of some approximate conditions is performed. These conditions can be optimized with respect to some free parameters which justifies the name of the method. The optimized Schwarz waveform relaxation method was first applied to the wave equation [START_REF] Gander | Optimal Schwarz waveform relaxation for the one dimensional wave equation[END_REF] and then to the advection-diffusion equation with constant or variable coefficients [START_REF] Martin | An optimized Schwarz waveform relaxation method for unsteady convection diffusion equation[END_REF]. A recent paper [START_REF] Gander | Optimized Schwarz Waveform Relaxation for Advection Reaction Diffusion Problems[END_REF] gives the complete solution of the one dimensional optimization problem for constant coefficients equations. More recently the method has been extended to the linearized viscous shallow water equations without advection term by V. Martin [START_REF] Martin | A Schwarz Waveform Relaxation Method for the Viscous Shallow Water Equations[END_REF]. Here we are interested in the application of the method to the system of Primitive Equations of the ocean. It leads to non-trivial new problems (new transmission conditions, well-posedness of the problem, convergence of the algorithm...) that we address in this article.

The outline of the paper is the following : in Section 2 we write the equations and we precise the asymptotic regime that we consider. In Section 3 we derive an approximated Dirichlet to Neumann operator, and define the associated Schwarz waveform relaxation algorithm. In Section 4 we define a weak formulation of the problem on the whole domain and prove that it is well-posed in the natural functional spaces. In Section 5 we introduce a weak formulation for the Schwarz waveform relaxation algorithm and prove that each sub-problem solved in the algorithm is well-posed. Finally we present some numerical results in Sections 6.

2. The set of equations. We first write the primitive equations of the ocean. Then we present the simplified system from which we are able to derive efficient transmission conditions.

2.1. The primitive equations of the ocean. We consider the primitive equations of the ocean on the domain (x, y, z, t) ∈ R×R×[-H(x, y), ζ(x, y, t)]×R + where -H(x, y) denotes the topography of the ocean and ζ(x, y, t) denotes the altitude of the free surface of the ocean. The primitive equations are commonly written [START_REF] Cushman-Roisin | Introduction to Geophysical Fluid Dynamics[END_REF]

∂ t U h + U h • ∇ h U h -ν∆U h + 2 ρ 0 Ω ∧ U h + 1 ρ 0 ∇ h p = 0, (2.1)
∇ h • U h + ∂ z w = 0, (2.2) ∂ z p = -ρg, (2.3) ρ = ρ(z, T, S), (2.4) ∂ t T + U 0 • ∇T -ν T ∆T = Q T , (2.5) ∂ t S + U 0 • ∇S -ν S ∆S = Q S , (2.6)
where the unknowns are the 3d-velocity (U h , w) = (u, v, w), the pressure p, the density ρ, the temperature T and the salinity S. The parameters are the gravity g, the eddy viscosity ν, the eddy diffusion coefficients for the tracers ν T and ν S and the earth rotation vector Ω. The source terms Q T and Q S for the temperature and salinity model the influence of the sun, rivers and atmosphere for these tracers. Note that we consider here the classical but non-symmetric viscosity tensor

σ =   σ xx σ xy σ xz σ yx σ yy σ yz σ zx σ zy σ zz   = ν   ∂ x u ∂ y u ∂ z u ∂ x v ∂ y v ∂ z v ∂ x w ∂ y w ∂ z w   .
Other form of the viscosity tensor can be found in [START_REF] Gerbeau | Derivation of viscous Saint Venant system for laminar shallow water; numerical simulation[END_REF]. Note also that it is possible to consider different viscosity coefficients in the horizontal and vertical directions [START_REF] Lucas | New developments and cosine effect in the viscous shallow water and quasi geostrophic equations[END_REF]. These equations are supplemented by initial and boundary conditions. At initial time, we impose

U h (•, 0) = U h,i in Ω, ζ(•, 0) = ζ i in ω,
where the subscript letters i means "initial". At the bottom of the ocean we impose a non-penetration condition and a friction law of Robin type (α b > 0)

U h (-H) • ∇ h (H) -w(-H) = 0, ∂ n U t (-H) + α b U t (-H) = 0, (2.7)
where U t stands for the tangential velocity and n denotes the outward normal vector to the bottom of the ocean. The free surface is transported by a kinematic boundary condition

∂ t ζ + U h (ζ) • ∇ h ζ -w(ζ) = 0. (2.8)
The equilibrium of the stresses at the free surface implies

[σ -(p -p a )Id] • 1 1 + (∂ x ζ) 2 + (∂ y ζ) 2   ∂ x ζ ∂ y ζ 1   = 0, (2.9)
where p a (x, y, t) denotes the atmospheric pressure.

Free surface

Ocean bottom

H h z = 0

2.2.

A linearized hydrostatic model. In order to derive simple and efficient transmission conditions for the Schwarz waveform relaxation method we make some assumptions on this set of equations.

First we neglect the influence of the tracers (temperature and salinity) on the density. Thus we suppose that the density is constant (we assume ρ 0 = 1) and we do not solve the equations on the tracers (2.5)-(2.6). Note that these equations are classical advection-diffusion equations for which the optimized transmission conditions are well known [START_REF] Gander | Optimized Schwarz Waveform Relaxation for Advection Reaction Diffusion Problems[END_REF][START_REF] Martin | An optimized Schwarz waveform relaxation method for unsteady convection diffusion equation[END_REF]. Then we use the divergence-free condition (2.2) and the non-penetration condition (2.7) to write the vertical velocity w as a function of the horizontal velocity U h and we use the hydrostatic assumption (2.3) to write the pressure p as a function of the water height ζ. The remaining unknowns in the system are the horizontal velocity U h and the water height ζ. The set of equations (2.1)-(2.3) stands

∂ t U h + U h • ∇ h U h -ν∆U h + f CU h + g∇ h ζ = 0, ∂ t ζ + ∇ h • ζ -H U h dz = 0, with C = 0 -1 1 0 , and f := 2 Ω • e z .
The first equation is written on the initial domain

R x ×R y ×[-H(x, y), ζ(x, y, t)] z ×R + t while the second one is written on R x × R y × R + t .
We consider for simplicity a flat bottom and a constant atmospheric pressure. Then we linearize the problem around a constant state which corresponds to a horizontal velocity U 0 = (u 0 , v 0 ) and a horizontal free surface located at z = 0. It follows that the water height ζ is a small perturbation. In the sequel U h denotes the perturbation on the horizontal velocity. The linearized problem stands

∂ t U h + U 0 • ∇ h U h -ν∆U h + f CU h + g∇ h ζ = 0, (2.10) ∂ t ζ + H∇ h • U h + U 0 • ∇ h ζ = 0, (2.11)
where

U h = u v := 1 H 0 -H U h dz,
denotes the mean horizontal velocity of the flow. This mean velocity is called barotropic velocity by the oceanographic community while the deviation U h -U h is called baroclinic velocity [START_REF] Cushman-Roisin | Introduction to Geophysical Fluid Dynamics[END_REF].

Note that the first equation (2.10) is now written in the fixed domain

R x × R y × [-H, 0] z × R + t .
The associated boundary conditions are

∂ z U h (z = 0) = 0, {∂ z U h + α b U h } (z = -H) = 0, (2.12)
where the boundary condition at z = 0 is deduced from the equilibrium of the stresses at the free surface (2.9). Indeed, with the help of the linearization procedure, we first deduce that at first order ∂ z U h (z = ζ) ≃ 0 and then ∂ z U h (z = 0) ≃ 0 since we assume that ζ is small. In order to derive the transmission conditions we assume α b = 0 in the sequel. However in the definition of the Schwarz waveform relaxation algorithm and in the numerical simulations, the condition α b > 0 will be supported.

2.3. Dimensionless system. We choose characteristic horizontal and vertical lengths (denoted L and H respectively) and velocity U of the problem. We introduce the dimensionless quantities

(x, y) = L(x, ỹ), t = (L/U ) t, ζ = H ζ, z = H z, U h = U U h , U 0 = U U 0 .
The spatial domains of computation are Ω = R x × R y × (-1, 0) z for the momentum equation and ω = R x × R y for the continuity equation. We study both equations in the time interval [0, T ], where T > 0 is fixed. Dropping the "∼" for a better readability, the system in dimensionless variables stands

∂ t U h + U 0 • ∇ h U h - 1 Re ∆ h U h - 1 Re ′ ∂ 2 z U h + 1 ε CU h + 1 F r 2 ∇ h ζ = 0, (2.13)
∂ z U h (x, y, 0, t) = ∂ z U h (x, y, -1, t) = 0, (2.14) U h (•, 0) = U h,i , (2.15) ∂ t ζ + U 0 • ∇ h ζ + ∇ h • U h = 0, (2.16) ζ(•, 0) = ζ i .
(2.17)

We have introduced the characteristic quantities,

-ε = U/(f L) the Rossby number, -Re = U L/ν the horizontal Reynolds number, -Re ′ = H 2 /L 2 Re the vertical Reynolds number, -F r = U/ √ gH the Froude number.
We choose to exhibit the Rossby number as a small parameter since we are interested in long-time oceanographic circulation for which the Rossby number is typically of magnitude 10 -2 . The values of Reynolds and Froude numbers vary with respect to the turbulent processes and to the depth of the area that is considered respectively.

3. The optimized Schwarz waveform relaxation algorithm. We are now interested in finding efficient transmission conditions for equations (2.13)-(2.17). We first present the Schwarz waveform relaxation method. Then we derive the relations satisfied by the optimal transmission conditions. Since we are not able to solve analytically these equations, we perform an asymptotic analysis with respect to the Rossby number ε in order to derive some approximated transmission conditions. Finally we present the related optimized Schwarz waveform relaxation algorithm.

3.1. The Schwarz waveform relaxation method. The heart of the method is the following. We first divide the computational domain into an arbitrary number of subdomains. Then we solve each sub-problem independently for the whole time interval. The interactions between neighboring subdomains are entirely contained in the boundary conditions. An iterative procedure is considered until a prescribed precision is reached. The advantages of the method are clear : the parallelization is almost optimal: at each step the sub-problems are solved independently, so the space-time discretization strategies (or even the models...) can be chosen independently on each subdomain. Moreover at the end of each step only a small amount of informations are exchanged. The main drawback is related to the needed number of iterations : the method is efficient if it converges quickly (in two or three iterations typically). This requirement needs the derivation of efficient transmission conditions. In the sequel we consider for simplicity two subdomains but the method extends to an arbitrary number of subdomains.

We begin with some notations. First we introduce the left and right spatial subdomains Ω -and Ω + defined by:

Ω -:= (-∞, 0) x × R y × (-1, 0) z , Ω + := (0, +∞) x × R y × (-1, 0) z ,
and their interface

Γ = {0} x × R y × (-1, 0) z ≃ R y × (-1, 0) z .
We also introduce the domains ω ± := ±(0, +∞) x × R y for the unknowns that do not depend on the z variable, and their interface γ

:= {0} x × R y ≃ R y .
Let D be some spatial open domain and T > 0 be a given real number. Then we will write D T to denote the cylindrical domain D T := D × (0, T ). We denote by PE the set of equations (2.13), (2.14) and (2.16) and X := (U h , ζ) stands for the solution of this system with associated initial data

X i := (U h,i , ζ i ).
Then the Schwarz waveform relaxation algorithm is defined as follows:

           PE(X n+1 -) = 0 on Ω - T , X n+1 -(•, 0) = X i on Ω -, B -X n+1 - = B -X n + on Γ T ,            PE(X n+1 + ) = 0 on Ω + T X n+1 + (•, 0) = X i on Ω + , B + X n+1 + = B + X n -on Γ T , (3.1
) where the operators B ± contain the transmission conditions. In the classical Schwarz waveform relaxation algorithm [START_REF] Gander | Overlapping Schwarz for parabolic problems[END_REF][START_REF] Daoud | Overlapping Schwarz waveform relaxation for convection reaction diffusion problems[END_REF], the transmitted quantities are of Dirichlet type and the operators B ± are thus chosen to be the identity operator. Note that in this case an overlap is needed in the definition of the subdomains.

In the sequel we are interested in deriving more efficient transmission conditions. In order to reach such a goal we will first describe the general method to obtain optimal transmission conditions. The transmission conditions are said to be optimal if the algorithm converges in two iterations to the solution of the initial problem. These optimal transmission conditions involve the Dirichlet to Neumann operator associated to PE on the subdomains Ω ± T . Here we will see that we are not able to obtain an explicit formulation for these optimal conditions. Anyway these optimal boundary conditions are not local and consequently too expensive to be useful from a numerical point of view. Recent methods have been developed recently in order to approximate these optimal conditions by analytical or numerical means -see the review paper [START_REF] Gander | Optimized Schwarz methods[END_REF] for elliptic problems and [START_REF] Gander | Optimized Schwarz Waveform Relaxation for Advection Reaction Diffusion Problems[END_REF] for parabolic evolution equations. Here we will perform an asymptotic analysis of the system with respect to the Rossby number ε in order to deduce a set of approximated and efficient transmission conditions. This strategy has been initiated in [START_REF] Martin | A Schwarz Waveform Relaxation Method for the Viscous Shallow Water Equations[END_REF] for the shallow water equation without advection term. In our case, it turns out that these approximate transmission conditions lie in a two parameter family of boundary conditions. In Section 6 we optimize numerically the transmission conditions in this two parameter family.

Let us first describe in a formal setting the ideal case of optimal transmission conditions for the Schwarz waveform relaxation algorithm (3.1). We consider the case u 0 > 0. The case u 0 < 0 is deduced by applying the symmetry "x ′ = -x". Integrating the linearized Primitive equations PE on a subdomain, we see that the flux of the unknown (U h , ζ) through the interface Γ is given by

1 Re ∂ x U h -u 0 U h - 1 F r 2
DN U h -X b := 1 Re ∂ x U h -u 0 U h - 1 F r 2 ζ 0 c |ΓT , where X = (U h , ζ) solves    PE(X) = 0 on Ω + T X(•, 0) = 0 on Ω + , X = X b on Γ T .
Symmetrically, consider a Dirichlet data U h,b , we set

DN U h + X b DN ζ + X b :=   1 Re ∂ x U h -u 0 U h - 1 F r 2 ζ 0 u 0 ζ + u   |ΓT , where X = (U h , ζ) solves    PE(X) = 0 on Ω - T , X(•, 0) = 0 on Ω -, U h = U h on Γ T .
Notice that since we consider the case u 0 > 0, the continuity equation (2.16) the boundary condition is relevant only in the subdomain Ω + T . This is why in the later case we do not have to prescribe a boundary condition for ζ.

Once these Dirichlet to Neumann operators are defined we can introduce the optimal transmission conditions

B -X = 1 Re ∂ x U h -u 0 U h - 1 F r 2 ζ 0 -DN U h -X , (3.2)
B + X =     - 1 Re ∂ x U h + u 0 U h + 1 F r 2 ζ 0 -DN U h + X u 0 ζ + u -DN ζ + X     , (3.3)
Proposition 3.1. With this particular choice of transmission operators B ± , the algorithm (3.1) converges in two iterations.

Proof. By linearity, we may assume that the exact solution is 0 (X i ≡ 0). At the initial step the solutions on each subdomain do not satisfy any particular property. But the first iterate solves the primitive equations with vanishing initial data. It follows from the very definition of the operators DN ± that in the definition of the second iterate, the right hand sides of the transmission conditions vanish for both sub-problems. We deduce that this second iterate vanish: the algorithm converges in two steps.

The operators (3.2)(3.3) being non-local pseudo-differential operator, they are not well suited for numerical implementation. Our strategy is to approximate these operators by numerically cheap operators. Of course the two-step convergence property will be lost. The quality of the approxiamation will be measured through the convergence rate of the algorithm. From the structure of (3.2)(3.3), we choose to write B ± as perturbations of the natural operators transmitted through the interface:

B -X = 1 Re ∂ x U h -u 0 U h - 1 F r 2 ζ 0 -S U h -X , (3.4)
B + X =     - 1 Re ∂ x U h + u 0 U h + 1 F r 2 ζ 0 -S U h + X u 0 ζ + u -S ζ + X     , (3.5)
where S U h ± and S ζ + are pseudo-differential operators that will approximate the Dirichlet to Neumann operators. Let us finally remark that the differences in the expression of the two transmission operators B ± are due to the sign u 0 > 0. Since B -contains the information that is transmitted from Ω + T to Ω - T it is constructed on three boundary values (velocities and water height) but it has to transmit only two boundary conditions for momentum equations (2.13). On the contrary B + is constructed on two boundary values (velocities) but has to send three boundary conditions (for momentum and continuity equations). In the next subsections we will identify optimal and approximated transmission operators. To carry out the computation of the Dirichlet to Neumann operators we perform Fourier-Laplace transforms.

3.2. Laplace-Fourier transform of the primitive equations. We perform on the set of primitive equations (2.13)-(2.16) a Fourier transform in the y variable and a Laplace transform in time. The dual variables are respectively denoted η ∈ R and s = σ + iτ ∈ C. The real part σ is assumed to be strictly positive. We obtain in each subdomain the same set of differential equations

s + u 0 ∂ x + iηv 0 - 1 Re ∂ 2 x + 1 Re η 2 - 1 Re ′ ∂ 2 z + 1 ε B Ûh + 1 F r 2 ∂ x iη ζ = 0, {s + u 0 ∂ x + iηv 0 } ζ + ∂ x û + iη v = 0.
In the z direction we introduce the eigenmodes of the operator -∂ 2 z on (-1, 0) with homogeneous Neumann boundary conditions (2.14) e n (z) := α n cos(µ n z) with µ n := nπ; α 0 := 1 and

α n := √ 2 if n > 0.
Then we search for the solution on the form

Ûh (x, z) = ∞ n=0 Û n h (x)e n (z).
Note that we obviously obtain Û h = Û 0 h . It means that the first vertical mode Û 0 h represents the barotropic velocity while the sum of the other ones denotes the baroclinic deviation. The barotropic mode is coupled with the water height and it is the solution of the following system of three ordinary differential equations,

- 1 Re ∂ 2 x Û 0 h + u 0 ∂ x Û 0 h + s + iηv 0 + 1 Re η 2 + 1 ε B Û 0 h + 1 F r 2 ∂ x ζ iη ζ = 0, (3.6) u 0 ∂ x ζ + (s + iηv 0 ) ζ + ∂ x û0 + iηv 0 = 0.(3.7)
This last system is exactly the Laplace-Fourier transform of the so-called linearized viscous shallow water equations [START_REF] De Saint-Venant | Théorie du mouvement non-permanent des eaux, avec application aux crues des rivières et à l'introduction des marées dans leur lit (in french)[END_REF].

For the other vertical modes we have a set of two coupled reaction advection diffusion equations,

- 1 Re ∂ 2 x Û n h + u 0 ∂ x Û n h + s + iηv 0 + 1 Re η 2 + 1 Re ′ µ 2 n + 1 ε B Û n h = 0. (3.8)
3.3. Optimal transmission conditions for the baroclinic modes. The derivation of optimal transmission conditions for an advection diffusion equation was performed in [START_REF] Martin | An optimized Schwarz waveform relaxation method for unsteady convection diffusion equation[END_REF]. Here we are interested in the set of coupled reaction advection diffusion equations (3.8). The baroclinic modes are not coupled with the evolution of the water height. Hence for these modes the transmission operators have two components and will be searched on the form

B n ± = ∓ 1 Re ∂ x U h ± u 0 U h -S u,n ± . (3.9)
We search for the solution of system (3.8) as a sum of exponentials x → e λx . Plugging this ansatz in the system, we obtain that e λx solves (3.8) if and only if λ is a root of the determinant of the matrix

M n (λ) :=      - λ 2 Re + u 0 λ + s + η 2 Re + µ 2 n Re ′ + iηv 0 - 1 ε 1 ε - λ 2 Re + u 0 λ + s + η 2 Re + µ 2 n Re ′ + iηv 0      .
This determinant is a polynomial of degree four in λ and we can compute its four roots

λ n,+ ± := Re 2 u 0 + ∆ n ± , λ n,- ± := Re 2 u 0 -∆ n ± , (3.10)
where

∆ n ± := u 2 0 + 4 Re η 2 Re + µ 2 n Re ′ + s + iηv 0 ± i ε . (3.11)
Every solution λ n,± ± is associated with a one dimensional kernel generated by the vector φ n,± ± defined by

φ n,± + = 1 -i , φ n,± - = 1 i .
Since the solutions must vanish at infinity we search for solutions in Ω -on the form

Û n h,-(x) = α n,+ + e λ n,+ + x φ n,+ + + α n,+ -e λ n,+ -x φ n,+ - = Φ n,+ • exp xΛ n,+ • α n,+ , (3.12)
where

Φ n,± := 1 1 -i i , Λ n,± := λ n,± + 0 0 λ n,± - , α n,± := α n,± + α n,± - . (3.13)
In Ω + we search for the solution on the form

Û n h,+ (x) = α n,- + e λ n,- + x φ n,- + + α n,- -e λ n,- -x φ n,- - = Φ n,-• exp xΛ n,-• α n,-. (3.14)
It follows from relations (3.12) and (3.14) that

∂ x Û n h,∓ (x) = Φ n,± • Λ n,± • exp xΛ n,± • α n,± = Φ n,± • Λ n,± • Φ n,± -1 • Û n h,∓ . (3.15)
We can now define the operator S u,n ± in (3.9) in order to derive an optimal algorithm. This is done through its Laplace-Fourier symbol:

Ŝu,n ± := ∓ 1 Re Φ n,± Λ n,± Φ n,± -1 ± u 0 Id. (3.16)
3.4. Approximate transmission conditions for baroclinic modes. Since we want to construct an efficient but simple Schwarz waveform relaxation algorithm we will derive approximated transmission conditions by considering an asymptotic analysis of the results of the previous subsection. The definition (3.11) of ∆ n ± leads to the expansion Ŝn

± = Ŝu,n ±, app + O( √ ε) with Ŝu,n ±, app := 1 2       ±u 0 - 2 Re 1 √ ε 2 Re 1 √ ε - 2 Re 1 √ ε ±u 0 - 2 Re 1 √ ε       . (3.17)
Note that the approximated operator (3.17) does not depend on n. Consequently the related approximated transmission operators (3.9) can be applied to the whole baroclinic velocity, i.e. to the sum of the baroclinic modes.

3.5. Approximate transmission conditions for the barotropic mode. The derivation of optimal transmission conditions for the linearized viscous shallow water equations without advection term was performed in [START_REF] Martin | A Schwarz Waveform Relaxation Method for the Viscous Shallow Water Equations[END_REF]. Here we are interested in the linearized viscous shallow water equations (3.6)-(3.7). The transmission operators will be searched on the form (3.4)-(3.5).

As for the baroclinic modes we search for the solution of system (3.6)-(3.7) as a sum of exponentials e λx . Here λ has to be a root of the determinant of the matrix

M 0 (λ) defined by          - λ 2 Re + u 0 λ + s + η 2 Re + iηv 0 - 1 ε λ/F r 2 1 ε - λ 2 Re + u 0 λ + s + η 2 Re + iηv 0 iη F r 2 λ iη s + u 0 λ + iηv 0          .
This determinant is a polynomial of degree five which does not admit a trivial decomposition. Hence it is not possible to derive an explicit formula for the solutions of (3.6)-(3.7). Consequently, we are not able to obtain an explicit form for the optimal transmission conditions for the barotropic mode, even in Fourier-Laplace variables.

In order to derive approximated transmission conditions we use the fact that the Rossby number is a small parameter to compute approximated values of the roots of the determinant of M 0 (λ). The related approximated transmission conditions will be coherent with the results of the previous subsection for the baroclinic modes. Since u 0 is positive we first notice that three roots (3.18)-(3.19) have a negative real part and two roots (3.20) have a positive real part. The negative roots will be denoted λ 0,- ± and λ 0 0 . The positive ones will be denoted λ 0,+ ± . The notations for the related quantities that we introduce later are coherent with the previous ones (3.13). As above, we search for the solution in Ω -on the form

X-(x) = α 0,+ + e λ 0,+ + x φ 0,+ + + α 0,+ -e λ 0,+ -x φ 0,+ - =: Φ 0,+ • exp xΛ 0,+ • α 0,+ .
In Ω + we search for the solution on the form

X+ (x) = α 0,- + e λ 0,- + x φ 0,- + + α 0,- -e λ 0,- -x φ 0,- -+ α 0 0 e λ 0 0 x φ 0 0 =: Φ 0,-• exp xΛ 0,-• α 0,-.
We compute the following approximations for the roots of the determinant of M 0 (λ):

λ 0 0 = - s + iηv 0 u 0 + O(ε 2 ), (3.18) λ 0,- ± = - √ ±iRe √ ε + Re u 0 2 - Re 4F r 2 u 0 + O( √ ε), (3.19)
λ 0,+ ± = √ ±iRe √ ε + Re u 0 2 - Re 4F r 2 u 0 + O( √ ε). (3.20)
The associated kernel is always one dimensional and spanned by:

Φ 0 0 =   0 0 1   + O(ε 2 ), Φ 0,- ± =              u 0 ± i √ 2Re 4F r 2 √ ε ±iu 0 ∓ i √ 2Re 4F r 2 √ ε -1 - ± √ 2 4 i √ Re u 0 F r 2 + √ 2 2 1 ± i √ Re ((±u 0 + iv 0 )η + s) √ ε              + O(ε), Φ 0,- ± =              u 0 - i √ 2Re 4F r 2 √ ε -iu 0 + i √ 2Re 4F r 2 √ ε -1 - - √ 2 4 i √ Re u 0 F r 2 + √ 2 2 1 -i √ Re ((-u 0 + iv 0 )η + s) √ ε              + O(ε).
As in the baroclinic modes case, we compute the approximated transmission operators in Laplace-Fourier variables by Ŝu,0

-, app = 1 Re Φ 0,-Λ 0,-Φ 0,--1 2,3 - u 0 0 -1 Re 0 u 0 0 , Ŝu,0 +, app = - 1 Re Φ 0,+ Λ 0,+ Φ 0,+ -1 + u 0 Id,
where M 2,3 denotes the first 2 × 3 matrix extracted from the 3 × 3 matrix M . It leads to the following Laplace-Fourier symbols Ŝu,0

-, app = 1 2       - √ 2 √ Reε -u 0 - 1 F r 2 u 0 √ 2 √ Reε + 1 2F r 2 u 0 -2 F r 2 - √ 2 √ Reε + 1 2F r 2 u 0 - √ 2 √ Reε -u 0 0       + O( √ ε), (3.21)
and

Ŝ0 +, app = 1 2           - √ 2 √ Reε + u 0 - 1 F r 2 u 0 √ 2 √ Reε - 1 2F r 2 u 0 - √ 2 √ Reε - 1 2F r 2 u 0 - √ 2 √ Reε + u 0 0 0           + O( √ ε).(3.22)
By using relations (3.17), (3.21) and (3.22), we notice that

Ŝu,0 ±, app 2,2 = Ŝu,n ±, app +     - 1 2F r 2 u 0 ∓ 1 4F r 2 u 0 ∓ 1 4F r 2 u 0 0     .
It follows that a part of the transmission conditions will be applied to the whole velocity (sum of baroclinic and barotropic modes) while a second part will be applied only to the barotropic mode. The first part corresponds to the operator Ŝu,n ±, app . The second one corresponds to the remaining terms in the operator Ŝu,0 ±, app . 3.6. The optimized Schwarz waveform relaxation algorithm. Thanks to the computed approximated operators (3.17), (3.21) and (3.22) we can now derive an approximated Schwarz waveform relaxation algorithm for the linearized primitive equations (2.13)-(2.17).

Since the computed operators (3.17), (3.21) and (3.22) do not depend neither on the Fourier variable η nor on the Laplace variable s the related operators in the real space are identical to their Laplace-Fourier symbols. It follows that the approximated transmission operators B ± (3.4)-(3.5) have the following form

B -X =         1 Re ∂ x u + √ 2 2 √ Re ε - u 0 2 u - √ 2v 2 √ Re ε + u -v/2 2F r 2 u 0 1 Re ∂ x v + √ 2 2 √ Re ε - u 0 2 v + √ 2u 2 √ Re ε - u 4F r 2 u 0         , (3.23)
B + X =             - 1 Re ∂ x u + ζ F r 2 + √ 2 2 √ Re ε + u 0 2 u - √ 2v 2 √ Re ε + u + v/2 2F r 2 u 0 - 1 Re ∂ x v + √ 2 2 √ Re ε + u 0 2 v + √ 2u 2 √ Re ε + u 4F r 2 u 0 u 0 ζ + u             ,(3.24)
for which we recall that u and v represent the mean-values with respect to the z variable of the velocities u and v. Note that by replacing the first component (B + X) 1 by the linear combination

(B + X) 1 -1/(F r 2 u 0 )(B + X) 3 ,
we replace (3.24) by the equivalent transmission conditions

B ∼ + X =             - 1 Re ∂ x u + √ 2 2 √ Re ε + u 0 2 u - √ 2v 2 √ Re ε - u -v/2 2F r 2 u 0 - 1 Re ∂ x v + √ 2 2 √ Re ε + u 0 2 v + √ 2u 2 √ Re ε + u 4F r 2 u 0 u 0 ζ + u             . (3.25)
In the sequel we use (3.25) rather than (3.24) and we drop the superscripts "∼".

Next, we remark that the transmission conditions (3.23)(3.25) are a particular case of the generalized transmission conditions

B -X = 1 Re ∂ x U h - u 0 2 U h + α √ ε AU h + βBU h , (3.26)
B + X =   - 1 Re ∂ x U h + u 0 2 U h + α √ ε AU h -βBU h u 0 ζ + u   (3.27)
where

A := 1 -1 1 1 , B := 1 -1/2 -1/2 0 . (3.28)
The original transmission operators (3.23)-(3.25) correspond to the choice

α = 1 √ 2Re , β = 1 2F r 2 u 0 . (3.29)
Notice that B -X and (B + X) [START_REF] Arakawa | Computational design of the basic dynamical processes of the UCLA general circulation model[END_REF][START_REF] Arango | ROMS : A Regional Ocean Modeling System[END_REF] do not depend on the water height ζ, so we may rewrite

B -X = B U h -U h and B + X = t (B U h + U h , B ζ + X) as B U h ± U h := ∓ 1 Re ∂ x U h ± u 0 2 U h + α √ ε AU h ∓ βBU h B ζ + X := u 0 ζ + u. (3.30)
Let us emphasize the identity:

B U h + U h + B U h -U h = 2 α √ ε AU h . (3.31)
This relation will be useful both for defining a weak formulation of the algorithm in Section 5 and for the numerical implementation of this algorithm in Section 6.

Finally the Schwarz waveform relaxation algorithm (3.1) writes

           PE(X n+1 -) = 0 on Ω - T , X n+1 -(•, 0) = X i on Ω -, B U h -U n+1 h,-= B U h -U n h,+ on Γ T ,                    PE(X n+1 + ) = 0 on Ω + T , X n+1 + (•, 0) = X i on Ω + , B U h + U n+1 h,+ = B U h + U n h,-on Γ T , B ζ + X n+1 + = B ζ + X n - on γ T . (3.32)
where the operators B U h ± , B ζ + are defined by equalities (3.28)(3.30) and where α and β are free parameters. These generalized transmission conditions can now be optimized with respect to the two parameters α and β. In the case of a one dimensional reaction advection diffusion equation this optimization problem has been solved analytically (see [START_REF] Gander | Optimized Schwarz Waveform Relaxation for Advection Reaction Diffusion Problems[END_REF]). Here, we will present a numerical procedure in Section 6.

4.

Well-posedness of the linearized Primitive Equations. In the previous sections we have performed formal computations on the linearized Primitive Equations leading to the construction of the Schwarz waveform relaxation algorithm (3.32). The aim of this section is to be more precise: we will define a weak formulation of the system (2.13)-(2.17) and then prove that this system is well-posed in the natural spaces associated to this weak formulation.

From now on we relax the boundary condition on the bottom, i.e. we assume α b ≥ 0 instead of α b = 0. Moreover, in order to prepare the study of the well posedness of the algorithm (3.1) in the next section, we consider non-homogeneous right-hand sides Y = (F 1 , F 2 , f ) = Y (x, y, z, t). The system of linearized primitive equations PE(X) = Y writes

∂ t + U 0 • ∇ h - 1 Re ∆ h - 1 Re ′ ∂ 2 z + 1 ε C U h + 1 F r 2 ∇ h ζ = F in Ω T , (4.1)
∂ z U h (x, y, 0, t) = 0 on ω T , (4.2) -∂ z U h (x, y, -1, t) + α b U h (x, y, -1, t) = 0 on ω T , (4.3) {∂ t + U 0 • ∇ h } ζ + ∇ h • U h = f in ω T . (4.4)
We supplement this system with the initial conditions

U h (•, 0) = U h,i , in Ω, (4.5) ζ(•, 0) = ζ i in ω. (4.6)
Note that if we consider that the water height ζ + is given, the system (4.1)-(4.3), (4.5) with unknown U h is a classical linear parabolic problem. On the other hand if we consider that the mean horizontal velocities U h are given then ζ solves the linear transport problem with source term (4.4), (4.6). We will proceed as follows: first we recall the classical weak formulations both for the parabolic problem (with prescribed water height) and for the transport equations (with prescribed velocity). These two problems define two maps S 1 : ζ → U h and S 2 : U h → ζ. Finally we define the weak solutions of the Primitive Equations to be the fixed points of the map τ : (U h , ζ) → (S 1 (ζ), S 2 (U h )) and conclude by proving the existence of a unique fixed point.

Let us first introduce some functional spaces and some notations. We will work with initial data and right hand sides satisfying

U h,i ∈ H := L 2 (Ω, R 2 ), ζ i ∈ L 2 (ω), F ∈ L 2 (0, T ; V ′ ), f ∈ L 2 (0, T ; L 2 (ω)),
where

V ′ is the topological dual of V := H 1 (Ω, R 2 ).
The weak solutions will satisfy

U h ∈ C ([0, T ], H) ∩ L 2 (0, T ; V), ζ ∈ C [0, T ], L 2 (ω) ∩ C(R x , L 2 (R y × (0, T ))) .
We will need the following bilinear forms:

a(U, V ) := 1 Re (∇ h U, ∇ h V ) ΩT + 1 Re ′ (∂ z U, ∂ z V ) ΩT + α b Re ′ (U, V ) ω-1,T + 1 ε (CU, V ) ΩT + (U 0 • ∇U, V ) ΩT , (4.7)
c(ζ, V ) := 1 F r 2 (ζe x , ∂ x V ) ΩT . (4.8)
where ω -1 := R x × R y × {-1} z and (U, V ) Σ denotes the L 2 scalar product on Σ.

Assuming that we have a strong solution, and taking the scalar product of equation (4.1) with V ∈ D(Ω × (0, T ), R 2 ), we obtain (after integrating by parts) the following weak formulation for the equations governing the horizontal velocities:

∀ V ∈ D(Ω × (0, T ), R 2 ), (∂ t U h , V ) ΩT + a(U h , V ) = c(ζ, V ) + F, V . (4.9)
We now state

Definition 4.1. Let F ∈ L 2 (0, T ; V ′ ) and ζ ∈ L 2 (ω T ), we say that U h ∈ L 2 (0, T ; V) is a weak solution of the system (4.1) if (4.9) holds. Proposition 4.2. Let U h,i ∈ L 2 (Ω), F ∈ L 2 (0, T ; V ′ ) and ζ ∈ L 2 (ω T ), there exists a unique weak solution U h ∈ C([0, T ]; H) ∩ L 2 (0, T ; V) of (4.1)
satisfying the initial condition (4.5). Moreover, we have the energy inequality

1 2 U h 2 Ω (t) + t 0 1 Re ∇ h U h 2 Ω (s) + 1 Re ′ ∂ z U h 2 Ω (s) + α b Re ′ U h 2 ω-1 (s) ds ≤ 1 2 U h,i 2
Ω + t 0 { F , U h (s) + (∂ x u, ζ) ω (s)} ds. (4.10)
Proof. The method is classical and we only sketch the proof. We obtain the existence of a solution satisfying (4.10) by the Galerkin method. Let (E m ) be an increasing sequence of finite dimensional sub-spaces of V such that ∪E m is dense in V. For every m, there exists U m ∈ C ∞ ([0, T ], E m ) such that (4.9) holds for every V ∈ D(0, T ; E n ) -we only have to solve a finite system of linear ordinary differential equations. Using U m × 1 [0,t] as a test function in (4.9), we conclude that U m satisfies (4.10). So using the Cauchy Schwarz inequality and the Grönwall Lemma, we see that the sequence (U m ) is uniformly bounded in L 2 (0, T ; V). Now from the weak formulation, we deduce that (∂ t U m ) is bounded in L 2 (0, T ; V ′ ), thus by Aubin-Lions Lemma, (U m ) is compact in C([0, T ], H). Extracting a subsequence we obtain a solution satisfying (4.10). Regularizing in time and using the weak formulation, we see that any solution satisfies (4.10) and uniqueness follows by the energy method.

Let us turn our attention to the equations (4.4), (4.6) governing the evolution of the water height ζ. This is a linear transport equation with constant coefficients and a source term. Assuming that ζ is a strong solution, multiplying (4.4) by a test function χ ∈ D(ω × [0, T )), integrating on ω T , integrating by parts in space and time and then using the initial condition (4.6), we obtain

∀χ ∈ D(ω × [0, T )), -(ζ, {∂ t + U 0 • ∇ h } χ) ωT = (ζ i , χ(•, 0)) ω + (f -∇ h • U h , χ) ωT . (4.11) Definition 4.3. Let f ∈ L 2 (ω t ), U h ∈ L 2 (0, T ; V) and ζ i ∈ L 2 (ω), we say that ζ ∈ L 2 (ω T
) is a weak solution of the system (4.4), (4.6) if (4.11) holds.

Remark that the test function does not necessarily vanish at time 0 and that the initial data is prescribed by the weak formulation.

Proposition 4.4. Let f ∈ L 2 (ω t ), U h ∈ L 2 (0, T ; V) and ζ i ∈ L 2 (Ω).
There exists a unique weak solution ζ ∈ L 2 (ω T ) of (4.4), (4.6). Moreover this solution is given by the characteristic formula:

ζ(x, y, t) = ζ i (x -u 0 t, y -v 0 t) + t 0 (f -∇ h • U h )(x -u 0 s, y -v 0 s, t -s)ds. (4.12)
This solution lies in C [0, T ]; L 2 (ω) ∩ C R x ; L 2 (R y × (0, T )) and satisfies the following estimates for every t ∈ [0, T ] and every x ∈ R,

ζ(•, •, t) ω ≤ ζ i ω + t 0 f -∇ h • U h ω (s)ds, (4.13)
ζ(x, •, •) γt ≤ 1 u 0 ζ i ω + t 0 f -∇ h • U h ω (s)ds . (4.14)
Proof. First, notice that the estimates (4.13) (4.14) are direct consequences of the characteristic formula (4.12). Next, remark that if the data ∇ h U h , f and ζ i are sufficiently smooth then the function ζ given by the formula (4.12) solves (4.11). Hence we obtain the existence of a solution of (4.11) by density. For the uniqueness, by linearity we may assume that the data ∇ h U h , f and ζ i vanish. Then let ψ ∈ D(ω) and ρ ∈ D([0, T )) and define the test function χ by χ(x, y, t) := -ψ(x -u 0 t, y -v 0 t) T t ρ(s)ds, so that:

∂ t χ + U 0 • ∇χ = ψ(x -u 0 t, y -v 0 t)ρ(t), and (4.11) yields 0 = ωT ζ(x, t)ρ(t)ψ(x -u 0 t, y -v 0 t) = ωT ζ(x + u 0 t, y + v 0 t, t)ρ(t)ψ(x, y).
Since this is true for every (ψ, ρ) ∈ D(ω) × D([0, T )), we have ζ ≡ 0 on ω T . Finally, we define the notion of weak solution for the linearized primitive equations.

Definition 4.5. Let Y = (F, f ) ∈ L 2 (Ω T , V ′ ) × L 2 (ω T ) and X i = (U h,i , ζ i ) ∈ L 2 (Ω) × L 2 (ω). We say that X = (U h , ζ) ∈ C(0, T ; H) × L 2 (ω T )
(•, 0) = U h,i . Theorem 4.1. Let Y = (F, f ) ∈ L 2 (Ω T , V ′ ) × L 2 (ω T ) and X i = (U h,i , ζ i ) ∈ L 2 (Ω) × L 2 (ω). There exists a unique weak solution X = (U h , ζ) ∈ (C(0, T ; H) ∩ L 2 (0, T ; V)) × L 2 (ω T ) of (4.1)-(4.6).
Proof. The right hand side Y and the initial data X i being fixed, Proposition 4.2 and Proposition 4.4 define two maps

S 1 : L 2 (ω T ) → C(0, T ; H) ∩ L 2 (0, T ; V), ζ → U h , and
S 2 : L 2 (0, T, V) → C [0, T ]; L 2 (ω) ∩ C R x ; L 2 (R × (0, T )) , U h → ζ.
Denoting by T the affine mapping (U h , ζ) → (S 1 (ζ), S 2 (U h )), the application X is a weak solution of (4.1)-(4.6) if and only if it is a fixed point of T in

E T := L 2 (0, T ; V) × C([0, T ], L 2 (ω)). Let X 1 , X 2 ∈ E T and let (U h , ζ) := X 1 -X 2 and ( Ũh , ζ) := T (X 1 ) -T (X 2 )
, by linearity, using (4.10), we get for 0

≤ t ≤ T , 1 2 Ũh 2 Ω (t) + t 0 1 Re ∇ h Ũh 2 Ω (s) + 1 Re ′ ∂ z Ũh 2 Ω (s) ds ≤ t 0 (∂ x ũ, ζ) ω (s)ds ≤ t 0 ∇ Ũh 2 Ω (s) ds 1/2 t 0 ζ 2 ω (s) ds 1/2
.

By Young inequality, we may absorb the term in ∇ Ũh in the left hand side and get:

Ũh 2 Ω (t) + t 0 ∇ Ũh 2 Ω (s) ds ≤ κt sup s∈[0,t] { ζ 2 ω (s)} for 0 ≤ t ≤ T, (4.15)
for some κ > 0. Now (4.13) and the Cauchy Schwarz inequality yield

ζ 2 ω (t) ≤ t t 0 ∇U h 2 Ω (s)ds, for 0 ≤ t ≤ T. (4.16)
Finally, inequalities (4.15) (4.16) imply that, for T ′ ∈ (0, T ] small enough, the mapping T is strictly contracting in E T ′ yielding the existence of a unique fixed point of T in E T ′ . Repeating the argument on the intervals [T ′ , 2T ′ ], [2T ′ , 3T ′ ], ... we obtain the result on [0, T ].

5.

Weak formulation and well-posedness of the Schwarz waveform relaxation algorithm. We study in this section the well-posedness of the algorithm (3.32). First, we will define weak formulations for the two sub-problems and prove that they are well-posed. We will pay a particular attention to the weak form of the transmission conditions. In particular we will establish that the solutions X n+1 ± of the n th step of the algorithm (3.32) are in the right spaces, allowing the construction of the transmission conditions for the next step.

As in the previous section, we also consider non-homogeneous right-hand sides Y = (F, f ). Every step of the algorithm may be split in the two following sub-problems.

First in the domain {x < 0}, we search for a solution X n+1 -:= X -= (U h,-, ζ -) solving the initial and boundary value parabolic problem,

                             . ∂ t + U 0 •∇ h - 1 Re ∆ h - 1 Re ′ ∂ 2 z + 1 ε C U h,-+ 1 F r 2 ∇ h ζ -= F in Ω - T , -∂ z U h,-(x, y, -1, t) + α b U h,-(x, y, -1, t) = 0 , ∂ z U h,-(x, y, 0, t) = 0 on ω - T , B U h -U h,-= B U h -U n h,+ on Γ T , U n+1 h,-(•, 0) = U h,i in Ω -, (5.1
) and the transport problem,

{∂ t + U 0 • ∇ h } ζ -+ ∇ h • U h,-= f in ω - T , ζ n+1 -(•, 0) = ζ i in ω -. (5.2)
In the right subdomain {x > 0} we search for a solution X n+1

+ := X + = (U h,+ , ζ + ) solving the initial and boundary value parabolic problem,                              . ∂ t + U 0 •∇ h - 1 Re ∆ h - 1 Re ′ ∂ 2 z + 1 ε C U h,+ + 1 F r 2 ∇ h ζ + = F in Ω + T , -∂ z U h,+ (x, y, -1, t) + α b U h,+ (x, y, -1, t) = 0 , ∂ z U h,+ (x, y, 0, t) = 0 on ω + T , B U h + U h,+ = B U h + U n h,- on Γ T , U h,+ (•, 0) = U h,i in Ω + , ( 5

.3) and the transport problem with entering characteristics on the boundary γ

T ,    {∂ t + U 0 • ∇ h } ζ -+ ∇ h • U h,+ = f in ω + T , B ζ + X + = B ζ + X n - on γ T , ζ + (•, 0) = ζ i in ω + .
(5.4)

To prove that these two sub-problems are well-posed, we proceed as in Section 4. First we study the parabolic problems with prescribed water heights: we introduce a weak formulation for these problems and prove that they are well-posed. Then we study the transport equations, introduce their weak formulations and establish their well-posedness. Finally, the solutions of the coupled parabolic-transport problems are obtained via a fixed point method.

As in Section 4, the initial data

X i (U h,i , ζ i ) satisfy U h,i ∈ H, ζ i ∈ L 2 (ω). We choose right hand sides Y = (F, f ) in L 2 (0, T ; H) × L 2 (ω t )
. (In section 4, we only assumed F ∈ L 2 (0, T ; V ′ ), but here this choice would cause difficulties at the interface). We will search for weak solutions

X ± = (U h,± , ζ ± ) in the spaces, U h,± ∈ C [0, T ], H ± ∩ L 2 (0, T ; V ± ), (5.5) ζ ± ∈ C [0, T ], L 2 (ω ± ) ∩ C(R ±,x , L 2 (R y × (0, T ) t )) , (5.6)
with

H ± := L 2 (Ω ± , R 2 ) and V ± := H 1 (Ω ± , R 2 ).
5.1. The parabolic problems. Let us define the weak-formulation for the parabolic problems (5.1) and (5.3). First we introduce the bilinear forms a ± and c ± :

a ± (U, V ) := 1 Re (∇ h U, ∇ h V ) Ω ± T + 1 Re ′ (∂ z U, ∂ z V ) Ω ± T + α b Re ′ (U, V ) ω ± -1,T + 1 ε (CU, V ) Ω ± T + (U 0 • ∇U, V ) Ω ± T , (5.7
)

c ± (ζ, V ) = 1 F r 2 ζe x , ∂ x V ω ± T ± 1 F r 2 ζe x , V γT , (5.8)
where

ω ± -1 := R ± x × R y × {-1} z .
Next, taking the scalar product of the first equation of (5.1) or (5.3) with some test map V ∈ D(Ω ± × (0, T ), R 2 ), we obtain:

(∂ t U h,± , V ) Ω ± ×(0,T ) + a ± (U h,± , V ) = c ± (ζ ± , V ) ∓ 1 Re (∂ x U h,± , V ) Γ + (F, V ) Ω ± T .
Then, using the transmission conditions to express ∂ x U h,± on Γ, we get

(∂ t U h,± , V ) Ω ± ×(0,T ) + a ± (U h,± , V ) + b ± (U h,± , V ) = c ± (ζ ± , V ) + B U h ± U n h,∓ , V Γ + (F, V ) Ω ± T . with b ± (U, V ) := ± u 0 2 (U, V ) Γ + α √ ε (AU, V ) Γ ∓ β(BU , V ) Γ .
(5.9)

We are still not satisfied with this weak formulation. Indeed, the knowledge of

∂ x U n h,∓
on the boundary Γ × (0, T ) is needed for defining the term (B U h ± U n h,∓ , V ) Γ in the right hand side. Unfortunately, (5.5) only gives: ∂ x U n h,∓ ∈ L 2 (Ω ± × (0, T )) which is not sufficient to define a trace. To overcome this difficulty, we use relation (3.31) to define recursively the terms (B

U h ± U n h,∓ , V ) Γ . Indeed, for strong solutions, we have on Γ T B U h ∓ U h,± (3.31) = -B U h ± U h,± + 2 α √ ε AU h,± (3.32) = -B U h ± U n h,∓ + 2 α √ ε AU h,± . Thus, identifying B U h ± U n h,∓ with a distribution B n ± ∈ L 2 (0, T ; W ′ )
, where W denotes the space H 1/2 (Γ, R 2 ) ; we obtain a weak formulation of the algorithm for the horizontal velocities: For the first step, we choose

B 0 ± ∈ L 2 (0, T ; W ′ ) (5.10)
Then for n ≥ 0, the horizontal velocity is defined by U n+1 h,± = U h,± where U h,± solves

∀ V ∈ D(Ω ± × (0, T ), R 2 ), (∂ t U h,± , V ) Ω ± ×(0,T ) + a ± (U h,± , V ) + b ± (U h,± , V ) = c ± (ζ ± , V ) + B n ± , V ΓT + (F, V ) Ω ± T , (5.11)
where a ± , c ± , and b ± are defined in (5.7)-(5.9). Once U n+1 h,± is known, we can define the boundary conditions for the next step in the opposite domain by

B n+1 ∓ := -B n ± + 2 α √ ε AU n+1 h,± |ΓT .
(5.12)

Notice that assuming that the maps U n+1 h,± satisfy (5.5) then their traces on Γ T are well defined in L 2 (0, T ; W) ⊂ L 2 (0, T ; W ′ ). Consequently, the transmission conditions B n+1 ∓ defined recursively by (5.12) stay in the space L 2 (0, T ; W ′ ). (5.6). Then there exists a unique U n+1 h,± = U h,± with regularity (5.5) satisfying (5.11) and the initial condition U n+1 h,± (0) ≡ U h,i on Ω ± . Moreover, we have the energy inequality

Proposition 5.2. Let U h,i ∈ H, F ∈ L 2 (0, T ; H), B n ± ∈ L 2 (0, T ; W ′ ± ) and ζ ± (= ζ n+1 ± ) satisfying
1 2 U h,± 2
Ω ± (t) + α √ ε ± u 0 /2 U h,± 2 Γt ∓ β u ± 2 γt -(u ± , v ± ) γt + t 0 1 Re ∇ h U h,± 2
Ω ± (s) + 1 Re ′ ∂ z U h,± 2
Ω ± (s) + α b Re ′ U h,± 2 ω ± -1 (s) ds ≤ 1 2 U h,i 2
Ω ± + (F , U h,± ) Ω ± t + B n ± , U h,± Γt + t 0 (∂ x u ± , ζ ± ) ω ± (s) ± (u ± , ζ ± ) γ ± (s) ds. (5.13)
Proof. We proceed as in the proof of Proposition 4.2: we apply the Galerkin method. Here we only check that the a priori inequality (5.13) is sufficient for applying this method. In order to bound the quadratic terms in the left hand side of (5.13) and the last term in the right hand side, we will use the inequality

U 2 Γ ≤ 2 U Ω ± ∂ x U Ω± , valid for U ∈ V ± . (To prove it, write |U (0, y, z)| 2 = 2 0 -∞ (∂ x U • U )(x ′ , y, z) dx ′ inte- grate on R y × (-1, 0) z
and use the Cauchy-Schwarz inequality). From this inequality, the Cauchy-Schwarz inequality, the Young inequality and the fact that the trace on Γ defines a continuous embedding Π : V ± → W, we see that (5.13) implies

U h,± 2 Ω ± (t) + t 0 ∇U h,± 2 Ω ± (s) ds -κ t 0 U h,± 2 Ω ± (s) ds ≤ κ U h,i 2 Ω ± + F 2 Ω ± t + t 0 B n ± 2 W ′ (s) ds + ζ ± 2 Ωt + ζ ± 2 γt (5.14)
for some κ > 0. Taking a Galerkin sequence (U m ) associated to (5.11), the elements of this sequence satisfy (5.13) and then inequality (5.14) and the Grönwall Lemma imply that this sequence is bounded in L 2 (0, T ; V). Extracting a subsequence (as in Proposition 4.2) we obtain a solution of (5.11). Then using the weak formulation satisfied by U m we see that (∂ t U m ) is bounded in L 2 (0, T ; V ′ ) and from Aubin-Lions Lemma (see e.g. [START_REF] Showalter | Monotone operators in Banach space and nonlinear partial differential equations[END_REF]), the sequence (U m ) is compact L 2 (0, T ; H s (Ω ± , R 2 )) for s < 1. Thus we may let m tend to ∞ in the quadratic boundary terms in the left hand side of (5.13).

The uniqueness follows from (5.14) and Grönwall Lemma.

5.2. The transport equations. We now consider that the velocities U h,± = U n+1 h,± are known and study the transport problems (5.2) and (5.4). We begin with the domain {x < 0}. Proceeding exactly as in Section 4, we obtain that a strong solution of Problem (5.2) satisfies

∀χ ∈ D(ω -× [0, T )), -(ζ -, {∂ t + U 0 • ∇ h } χ) ω - T = (ζ i , χ(•, 0)) ω -+ (f -∇ h • U h,-, χ) ω - T . (5.15) Definition 5.3. Let f ∈ L 2 (ω t ), U h,-(= U n+1 h,-) ∈ L 2 (0, T ; V -) and ζ i ∈ L 2 (ω). We say that ζ -∈ L 2 (ω - T ) is a weak solution of Problem (5.
2) if (5.15) holds. The following result is proved exactly as Proposition 4.4

Proposition 5.4. Let f ∈ L 2 (ω t ), U n+1 h,-∈ L 2 (0, T ; V -) and ζ i ∈ L 2 (ω). There exists a unique weak solution ζ n+1 - = ζ -∈ L 2 (ω - T ) of (5.2)
. Moreover this solution is explicitly given by the formula:

ζ -(x, y, t) = ζ i (x -u 0 t, y -v 0 t) + t 0 (f -∇ h • U h,-)(x -u 0 s, y -v 0 s, t -s)ds. (5.16) It lies in C [0, T ]; L 2 (ω -) ∩ C (-∞, 0] x ; L 2 (
R y × (0, T )) and satisfies the following estimates for every t ∈ [0, T ] and every x ≤ 0,

ζ -(•, t) ω -≤ ζ i ω -+ t 0 f -∇ h • U h,-ω -(s)ds, (5.17)
ζ -(x, •) γt ≤ 1 u 0 ζ i ω -+ t 0 f -∇ h • U h,-ω -(s)ds . (5.18)
Once the solutions of (5.1)-(5.2) are known it is possible to define the transmission conditions on the water-height for the next step (see (3.30))

B ζ + X n+1 - := u 0 ζ n+1 -(0, •)) + u n+1 -(0, •). (5.19)
In the domain x > 0, the situation is slightly different since there are ingoing characteristics on γ T . So we choose test functions that do not necessarily vanish on the boundary and use the transmission condition to prescribe the value of the solution on γ T . Finally, a solution of (5.4) satisfies

∀χ ∈ D(ω + × [0, T )), -(ζ + , {∂ t + U 0 • ∇ h } χ) ω + T = (ζ i , χ(•, 0)) ω + + (ζ b , χ(0, •)) Rt + (f -∇ h • U h,+ , χ) ω + T , (5.20)
where the boundary value ζ b is defined on γ T by

ζ b := 1 u 0 B ζ + X n --u + .
(5.21)

Definition 5.5. Let f ∈ L 2 (ω t ), U h,+ (= U n+1 h,+ ) ∈ L 2 (0, T ; V + ), ζ i ∈ L 2 (ω).
Assuming that ζ b defined by (5.21) belongs to L 2 (γ T ) , we say that ζ n+1 T ) of (5.4). Moreover it is given by the characteristic formula:

+ = ζ + ∈ L 2 (ω + T )
ζ + (x, y, t) = ζ i (x -u 0 t, y -v 0 t) + t 0 (f -∇ h • U h,+ )(x -u 0 s, y -v 0 s, t -s)ds if x > u 0 t, and ζ + (x, y, t) = ζ b y - v 0 u 0 x, t - x u 0 + x u 0 0 (f -∇ h • U h,+ )(x -u 0 s, y -v 0 s, t -s)ds, with ζ b given by (5.21), if x ≤ u 0 t. The solution belongs to C [0, T ]; L 2 (ω + ) ∩ C [0, +∞) x ; L 2 (R y × (0, T ))
and satisfies the following estimates for every t ∈ [0, T ] and every x ≥ 0,

ζ + (•, t) ω + ≤ ζ i ω + + u 0 ζ b γt + t 0 f -∇ h • U h,+ ω + (s)ds, (5.22)
ζ + (x, •) γt ≤ 1 u 0 ζ i ω + + u 0 ζ b γt + t 0 f -∇ h • U h,+ ω + (s)ds . (5.23)
5.3. Well-posedness of the algorithm. First we define a weak formulation for the left and right sub-problems at step n of the algorithm. (5.4) if it has regularity (5.5)-(5.6) and if U h,+ , (respectively ζ + ) is a weak solution of (5.3)(respectively (5.4)).

Definition 5.7. Let Y = (F, f ) ∈ L 2 (Ω T ) × L 2 (ω T ), let X i = (U h,i , ζ i ) ∈ L 2 (Ω) × L 2 (ω). For n ≥ 0. • Let B n -∈ L 2 (0, T ; W ′ ). Then X n+1 - = (U h,-, ζ -) is a weak solution of Prob- lem (5.1), (5.2) if it has regularity (5.5)-(5.6) and if U h,-(respectively ζ -) is a weak solution of (5.1)(respectively (5.2)). • Let B n + ∈ L 2 (0, T ; W ′ ) and B ζ + X n -∈ L 2 (γ T ). Then X n+1 + = (U h,+ , ζ + ) is a weak solution of Problem (5.3),
Then we give a weak formulation for the complete algorithm.

Definition 5.8. The weak formulation of Algorithm (3.32) is defined by

• Choose B ζ + X 0 -∈ L 2 (γ T ) and B 0 ± ∈ L 2 (0, T ; W ′ ). Then, for n ≥ 0,
• Find X n+1 -weak solution of (5.1)-(5.2) and X n+1 + weak solution of (5.3)-(5.4).

• Define the transmission conditions for step n + 1 by (5.12) and (5.19).

Theorem 5.1. With the hypotheses of Definition 5.7, there exists a unique weak solution X n+1 -(respectively X n+1 + ) of Problem (5.1),(5.2) (respectively (5.3),(5.4)). Proof. We only prove the result for the left sub-problem, the other one being similar. As in the proof of Theorem 4.1, we use a fixed point method. Let us introduce the spaces

E 1 T := C [0, T ], H - ∩ L 2 (0, T ; V -), E 2 T := C [0, T ], L 2 (ω -) ∩ C((-∞, 0] x , L 2 (R y × (0, T ) t )) .
Proposition 5.2 (respectively Proposition 5.4) defines an affine mapping S 1 :

E 2 T → E 1 T , ζ -→ U h,-(respectively S 2 : E 1 T → E 2 T , U h,-→ ζ -). Setting E - T := E 1 T ×E 2 T -, an application X n+1
is a weak solution of Problem (5.1),(5.2) if and only if it is a fixed point in E - T of the mapping

T -: (U h,-, ζ -) → (S - 1 (ζ -), S - 2 (U h,-)).
We now show that T -has a unique fixed point. Let

X 1 , X 2 ∈ E - T and let (U h,-, ζ -) := X 1 -X 2 and ( Ũh , ζ) := T -(X 1 ) -T -(X 2 ). By linearity (5.14) yields: for 0 ≤ t ≤ T , Ũh,- 2 Ω -(t) + t 0 ∇ Ũh,- 2 Ω -(s) ds -κ t 0 Ũh,- 2 Ω -(s) ds ≤ κ ζ - 2 Ωt + ζ - 2 γt .
And from Grönwall lemma, we obtain for 0 ≤ t ≤ T , Ũh,-

2 Ω (t) + t 0 ∇ Ũh,- 2 Ω (s) ds ≤ κe κT t sup [0,t] ζ -(•, s) 2 ω + sup R- ζ -(w, •) 2 γt . (5.24)
Now from (5.17) and (5.18), we get

ζ- 2 ω -(t) + u 0 ζ-(x, •) γt ≤ t ∇U h,- 2 Ω - t for 0 ≤ t ≤ T. (5.25)
Finally, we endow

E - t with the norm (U h,-, ζ -) E - t := sup [0,t] U h,-(•, s) 2 Ω -+ ∇U h,- 2 Ω - t + sup [0,t] ζ -(•, s) 2 ω + 2κe κT sup R- ζ -(x, •) 2 Rt 1/2
.

With this norm (5.24) (5.25) imply that for T ′ ∈ (0, T ] small enough, T -is contracting in E T ′ . This yields the existence of a unique fixed point of T -in E - T ′ . We obtain the result on [0, T ] by continuation.

Finally, we can state Theorem 5.2. The algorithm (5.8) is well-defined.

Proof. We only have to check that for each step the hypotheses of Theorem 5.1 are satisfied. The solutions X n+1

± , build at step n have regularity (5.5), (5.6). We easily deduce that B n+1 ± defined by (5.12) belongs to L 2 (0, T ; W ′ ) and B ζ + X n+1 -defined by (5.19) belongs to L 2 (γ T ). Thus the hypotheses of Theorem 5.1 hold for step n + 1.

Remark 5.9. Although we do not exhibit a proof here, we are able to establish the convergence of the algorithm in some cases. More precisely, if the matrices A and B defined by (3.28) are replaced by diagonal matrices à and B, à being positive definite and B being non negative, then the algorithm converges. The proof relies on the energy method developed for the Shallow water equations without advection term in [START_REF] Martin | A Schwarz Waveform Relaxation Method for the Viscous Shallow Water Equations[END_REF]. Modifying slightly the proof, we can allow à and B to have non vanishing skew-symmetric off-diagonal parts. This generalization still does not cover the situation (3.28) because B has a symmetric non vanishing off-diagonal part. Nevertheless, numerical evidences of the convergence of the algorithm are given in the next section.

6. Numerical results.

6.1. Numerical scheme in the subdomains. For the numerical applications we consider for simplicity a 2 dimensional domain and the related two dimensional (x, z) version of the primitive equations (2.13)-(2.17). Note that the transmission conditions (3.26)-(3.27) are independent of the transverse y-variable and are not affected by this simplification.

In this subsection we do not deal with the boundary conditions. Hence the processes are the same in both subdomains Ω ± and we restrict ourselves to the subdomain Ω + . We first describe the space discretization of the subdomains. We consider a regular cartesian grid of nx × nz points and we apply a finite volume method. We introduce the horizontal space step ∆x and the vertical space step ∆z. For Euler or Navier-Stokes type problems it is well known that a good way to recover some numerical stability is to compute velocities and pressure on different cells (see for instance [START_REF] Arakawa | Computational design of the basic dynamical processes of the UCLA general circulation model[END_REF] and the publications devoted to the so-called C-grids). Here we only deal with the horizontal velocity and the water height (depending only depending on x and t) plays the role of the pressure. We thus have to introduce two types of finite volume meshes -see Figure 6.1. The first one is a 2d finite volume mesh and is related to the computation of the velocities. For i = 0...n x -1 and j = 0...n z we denote I = i + jn x . The cells of this first mesh will be denoted C + I = X + I + (-∆x/2, ∆x/2) × (-∆z/2, ∆z/2). where the points X + I stand for X + I = (0, -H) + (i∆x, j∆z) (they are represented by a black circle in Figure 6.1). The second grid is a 1d finite volume mesh devoted to the computation of the water height. The cells of this second mesh will be denoted c i+1/2 = x i+1/2 + (-∆x/2, ∆x/2). where the points x i+1/2 stand for x i+1/2 = (i + 1/2)∆x (they are represented by a circle with a number inside in Figure 6.1).

Let us now consider the discretization of the equations. Let us start with momentum equation (2.13). We integrate it on the time-space cell [t k , t k+1 ] × C I . We compute the interface fluxes at time t k+1/2 by classical centered formulas. We recover the well-known Crank-Nicolson scheme. It is known to be second order accurate and conditionally stable in the L ∞ norm under a CFL type condition on the time step ∆t k = t k+1 -t k . This strategy is applied for all the velocity nodes such that the neighboring nodes are included inside the considered subdomain. The related discrete x u I,k+1 -

1 Re ′ D 2 z u I,k+1 - 1 ε v I,k+1 + 1 F r 2 ∆x D x1 ζ j,k+1 = u I,k - ∆t 2 u 0 D x u I,k - 1 Re D 2 x u I,k - 1 Re ′ D 2 z u I,k - 1 ε v I,k + 1 F r 2 D x1 ζ j,k . (6.1) v I,k+1 + ∆t 2 u 0 D x v I,k+1 - 1 Re D 2 x v I,k+1 - 1 Re ′ D 2 z v I,k+1 + 1 ε u I,k+1 = v I,k - ∆t 2 u 0 D x v I,k - 1 Re D 2 x v I,k - 1 Re ′ D 2 z v I,k + 1 ε u I,k . (6.2)
where Let us now consider the mass equation (2.16). We integrate it on time space cells [t k , t k+1 ] × c i+1/2 -except for i = 0 where we need to use the transmission conditions. We compute the interface fluxes by using explicit upwind formulas. The resulting scheme is known to be first order and also conditionally stable under a CFL type condition. The related formula stands

D x u I,k = (u I+1,k -u I-1,k )/∆x
ζ i+1/2,k+1 = 1 - ∆t ∆x u 0 ζ i+1/2,k + ∆t ∆x u 0 ζ i-1/2,k - ∆t∆z ∆x u +,n+1 i,k -u +,n+1 i-1,k where u +,n+1 i,k = u j,k 2 + nz-1 j=1 u jnx+j,k + u nznx+j,k 2
denotes the discrete mean velocity of the flow along the vertical direction.

6.2. Numerical discretization near the boundaries. We now have to explain how we compute the numerical solution when one of the interfaces of the cell belongs to the physical boundaries of the domain or to the fictitious one that is related to the domain decomposition method. For all cases we choose to work in the same finite volume framework that we use in the interior of the subdomains.

For the physical boundary conditions (2.14) we use the ghost cells method. This method consists in introducing a fictitious cell along the boundary and then using the same finite volume strategy as in the interior domain. For no slip conditions (2.14) we choose the values of the unknowns in the fictitious cell to be equal to their values in the neighboring interior cell.

Let us now focus on the numerical treatment on the cells that are connected with the interface Γ = ∂Ω + ∩ ∂Ω --see Fig. 6.1. Here we will use a discrete version of the transmission conditions (3.26)- (3.27). This discrete information will be the only data that will be transmitted from a subdomain to the other one. Let us first consider the mass equation (2.16). We integrate it on the cell c 3/2 to obtain

∆x ζ +,n+1 1/2,k+1 -ζ +,n+1 -1/2,k + ∆t u 0 ζ +,n+1 1/2,k -ζ +,n+1 -1/2,k + u +,n+1 1,k -u +,n+1 0,k = 0.
where ζ +,n+1 1/2,k+1 denotes the water height computed in cell c + 1/2 at time t k+1 and for iteration n + 1 of the algorithm. The quantities ζ +,n+1 -1/2,k and u +,n+1 0,k have to be considered as unknown quantities since the corresponding cells are not included in Ω + . We will use the transmission conditions (3.32) to evaluate them. Hence we obtain thanks to (3.30)

ζ +,n+1 1/2,k+1 = 1 - u 0 ∆t ∆x ζ +,n+1 1/2,k - ∆t ∆x u +,n+1 1,k + ∆t ∆x B ζ,n +,k
where B ζ,n +,k has been computed in Ω -during the previous Schwarz iteration and is given by

B ζ,n +,k = u 0 ζ -,n nx,k + u -,n nx,k
The basic idea is the same for the momentum equation (2.13). Here we integrate the equation on the semi-cell C+ I = X + I + (0, ∆x/2) × (-∆z/2, ∆z/2). for I = jn x with j = 0, ..., n z . We obtain ∆x∆z 2 (u +,n+1 I,k+1 -u +,n+1

I,k ) +∆t∆z u 0 2 u +,n+1 I,k (r) + u +,n+1 I,k+1 (r) 2 -u +,n+1 I,k+1/2 (l) - 1 Re D x u +,n+1 I,k (r) + D x u +,n+1 I,k+1 (r) 2 -∂ x u +,n+1 I,k+1/2 (l) + 1 F r 2 ζ +,n+1 1/2,k + ζ +,n+1 1/2,k+1 2 -ζ +,n+1 -1/2,k+1/2 - ∆t∆x∆z 2 1 Re ′ D 2 z u +,n+1 I,k + D 2 z u +,n+1 I,k+1 2 - ∆t∆x∆z 2 1 ε v +,n+1 I,k + v +,n+1 I,k+1 2 = 0
where l (respectively r) denotes quantities that are evaluated on the left (respectively right) boundary of the cell C+ I . Hence quantities u +,n+1 I,k+1/2 (l), ∂ x u +,n+1 I,k+1/2 (l) and ζ +,n+1 -1/2,k+1/2 are unknown quantities since they involve quantities that are computed outside the domain Ω + . Here also we use transmission conditions (3.32) and we obtain thanks to (3.30)

u +,n+1 I,k+1 ∆t + 1 2 u 0 u +,n+1 I,k+1 (r) ∆x - 1 Re 2D x u +,n+1 I,k+1 (r) ∆x - 1 Re ′ D 2 z u +,n+1 I,k+1 - 1 ε v +,n+1 I,k+1 - 1 ∆x - α √ ε u +,n+1 I,k+1 + α √ ε v +,n+1 I,k+1 -( 1 F r 2 u 0 -β)ū +,n+1 1,k+1 - β 2 v+,n+1 1,k+1 = - 2 ∆x B u,n +,j,k + 1 ∆x - α √ ε u +,n+1 I,k + α √ ε v +,n+1 I,k -( 1 F r 2 u 0 -β)ū +,n+1 j,k - β 2 v+,n+1 j,k - 1 2 u 0 u +,n+1 I,k (r) ∆x - 1 Re 2D x u +,n+1 I,k (r) ∆x - 1 Re ′ D 2 z u +,n+1 I,k - 1 ε v +,n+1 I,k + u +,n+1 I,k ∆t - 1 F r 2 ζ +,n+1 1/2,k + ζ +,n+1 1/2,k+1 ∆x (6.3)
where B u,n +,j,k has been computed in Ω -during the previous Schwarz iteration and is deduced from relation (3.31)

B u,n +,j,k = B u,n-1 -,j,k + 2 α √ ε u -,n jnx,k + u -,n jnx,k+1 2 -2 α √ ε v -,n jnx,k + v -,n jnx,k+1 2
Note that B u,n-1 -,j,k is known since it has been computed in Ω + at iteration n-1 and has been transmitted to the domain Ω -before iteration n. Same type of computations for the transverse component of the velocity lead to the following scheme v +,n+1

I,k+1 ∆t + 1 2 u 0 v +,n+1 I,k+1 (∆x) ∆x - 1 Re 2D + x v +,n+1 I,k+1 ∆x - 1 Re ′ D 2 z v +,n+1 I,k+1 + 1 ε u 2,n+1 I,k+1 (6.4)
- 1 ∆x - α √ ε u +,n+1 I,k+1 - α √ ε v +,n+1 I,k+1 - β 2 ū+,n+1 1,k+1 = - 2 ∆x B v,n +,j,k + 1 ∆x - α √ ε u +,n+1 I,k - α √ ε v +,n+1 I,k - β 2 ū2,n+1 j,k + v +,n+1 I,k ∆t - 1 2 u 0 v +,n+1 I,k (r) ∆x - 1 Re 2D x v +,n+1 I,k (r) ∆x - 1 Re ′ D 2 z v +,n+1 I,k + 1 ε u +,n+1 I,k
where B v,n +,j,k has been computed in Ω -during the previous Schwarz iteration and is deduced from relation (3.31)

B v,n +,j,k = B v,n-1 -,j,k + 2 α √ ε u -,n jnx,k + u -,n jnx,k+1 2 + 2 α √ ε v -,n jnx ,k + v -,n jnx,k+1 2
The derivation of the discrete boundary condition in Ω -is based on the same type of computations. Note that only the components of the velocity are concerned by the transmission problem in Ω -.

6.3. Numerical optimization of the transmission conditions. In this section we are interested in the optimization the transmission conditions (3.30) with respect to the free parameters α and β. To optimize the conditions means that we choose parameters α and β such that the Schwarz waveform relaxation algorithm (3.32) reaches a given error for as small as possible number of iterations. The analytical solution of this problem is quite complex in the considered framework and we only present here a numerical strategy to reach the optimum. In the simpler case of a 1D advection diffusion equation a complete solution of the related optimization problem is given in [START_REF] Gander | Optimized Schwarz Waveform Relaxation for Advection Reaction Diffusion Problems[END_REF].

We consider a test case for which all the initial data (velocities and perturbation of the water height) are taken equal to zero. We initialize the algorithm (3.1) with random boundary conditions on the interface and we study the convergence of the solution towards the analytical ones. This test is quite classical to study the convergence of a domain decomposition algorithm. It is interesting since the initial quantities do contain all frequencies. In all the computations the physical parameters Re and F r are taken equal to one but the Rossby number ε remains free. For a given value of ε we apply the transmission conditions (3.30) for several values of the parameters α and β and we compare the L 2 error between the computed and the analytical solutions after a given number of iterations. It allows us to find an optimal pair (α opt , β opt ) that minimizes this error. This first study exhibit that the influence of the parameter β is quite small. In the following this parameter will be kept equal to its theoretical value (3.29). In a second step we study the dependency of the optimal parameter α opt with respect to the Rossby number ε. The results are presented in Fig. 6.2 for different values of ε. We found that this optimized parameter does depend on ε in a nontrivial way. We now present the evolution of the error on the computed solution as a function of the number of iterations of the Schwarz waveform relaxation algorithm (3.1) in both cases α = α opt and α = α T ay . in Fig. 6.3 we present the results for two different values of the Rossby number ε : ε = 10 -3 and ε = 10 -2 . The curves (Log of the error) all look like straight lines, at least after a sufficiently large number of iterations. The method appears to be more efficient when the Rossby number is smaller since the error decreases much faster in the case ε = 10 -3 -Fig. 6.3 on the left. This result is consistent with the previous theoretical study that is based on an asymptotic analysis in ε. We also observe that for a given value of ε the curves look similar for both optimized and Taylor approximation parameters even if the error decreases faster for the optimal value α opt . Moreover let us observe that to reach an error of 10 -4 (that is enough for the applicability of the Schwarz waveform relaxation algorithm) both algorithms (with optimized or Taylor approximation parameter) need a very close number of iterations. We compute the same test with Rossby number ε = 10 -2 but with a sinusoidal initial guess (instead of the random ones) for the transmission conditions. We consider two different sinusoids with one or ten periods in the space-time considered interval and we use Taylor approximation parameters α T ay and β T ay . In Fig. 6.4 the results appears to be much better for the low frequency sinusoid as for high frequency one. The results for the high frequency sinusoid look similar to the results that were obtained with the random initial guess. It follows that the method is particularly well adapted to low frequency signals : the relative error is smaller than 10 -4 after only two iterations. and with optimal parameter for a low frequency signal (down) and for a high frequency signal (middle) and for a random signal (up) 6.4. Numerical application. In this section we consider the case of a flow with a constant positive background velocity u 0 = 1.m/s and an initial local decreasing step on the water height. We choose the Rossby number ε equal to 10 -3 . We choose nx = 40, nz = 10 and nt = 40 in order to ensure the CFL condition. We present the initial solution and the solution computed at final time T = 1.3s after 20 iterations by the proposed Schwarz waveform relaxation algorithm in Fig. 6.6. The 2d horizontal velocity vector field (u, v) is presented in the 2d vertical domain (in the (x, z) plane) which is occupied by the flow. A horizontal vector denotes a velocity which is collinear to the x-direction and a vertical one denotes a velocity which is collinear to the ydirection. Since we consider the linearized version of the equations the step just moves without deformation from the left to the right of the domain. Since the Coriolis effect is dominant we observe the formation of a transverse jet which moves with the step. Another consequence of the Coriolis effect is the formation of a stationary eddy at the initial location of the step. We now compare the solution that is computed on the whole domain with the solution that is obtained by considering the presented domain decomposition strategy. In Fig. 6.7 we present the evolution of the relative error between the two solutions versus the number of considered iterations. It exhibits the fast convergence of the algorithm for such a case. After two iterations the relative error is around 10 -6 and it reaches the factor 10 -10 after eight iterations.

Conclusion.

We presented in this article a new domain decomposition method for the viscous primitive equations. It involves a Schwarz waveform relaxation type algorithm with approximated transmission conditions for which we proved well-posedness. We presented a numerical optimization of the transmission conditions and we study the speed of convergence of the algorithm for several test cases. Academic numerical In forthcoming papers we plan to prove the convergence of the algorithm and we want to present oceanographic configurations and to increase the efficiency of the algorithm by deriving more complex transmission conditions based on another asymptotic regime that corresponds to quasi-geostrophic flows.

Fig. 2 . 1 .

 21 Fig. 2.1. Schematic representation of the ocean

  Using this flux as a Neumann operator, we define the Dirichlet to Neumann operators as follows. Consider a Dirichlet data X b = (U h,b , ζ b ), we set

Definition 5 . 1 .

 51 Assuming that the functions ζ ± = ζ n+1 ± are known, the weak formulation of the parabolic part (5.1) and (5.3) of the algorithm (3.32) are defined as follows:

  is a weak solution of Problem (5.4) if (5.20) holds. Using the characteristic method, we have Proposition 5.6. Let f , U h,+ (= U n+1 h,+ ), ζ i and ζ b be as in Definition 5.5. There exists a unique weak solution ζ n+1 + ∈ L 2 (ω +

Fig. 6 . 1 .

 61 Fig. 6.1. Space discretization of Ω +

  denotes a classical approximation of the first derivative in space in horizontal direction,D 2 x u I,k = (u I+1,k -2u I,k + u I-1,k )/(2∆x) and D 2 z u I,k = (u I+nx,k -2u I,k + u I-nx,k)/(2∆x) denote classical approximations of second derivatives in space in horizontal and vertical directions, respectively.

Fig. 6 . 2 .

 62 Fig. 6.2. Quotient αopt/α T ay between the numerically optimized parameter and the Taylor approximation parameter as a function of the Rossby number ε (in Log scale)

Fig. 6 . 3 .

 63 Fig. 6.3. Log of the error on the computed solution as a function of the number of Schwarz iterations for Rossby number ε = 10 -3 (left) and ε = 10 -2 (right) and for a random initial guess using the Taylor approximation parameter α T ay (up) and the optimized one αopt (down)

Fig. 6 . 4 .

 64 Fig.6.4. Log of the error on the computed solution as a function of the number of Schwarz iterations for Rossby number ε = 10 -2 and with optimal parameter for a low frequency signal (down) and for a high frequency signal (middle) and for a random signal(up)

Fig. 6 . 5 .

 65 Fig. 6.5. Water height and velocity field at initial time

Fig. 6 . 6 .

 66 Fig. 6.6. Water height and velocity field at final time

Fig. 6 . 7 .

 67 Fig. 6.7. Log of the relative error on the solution computed by using the Schwarz waveform relaxation algorithm versus number of iterations

Acknowledgements. The authors thank L. Halpern and V. Martin for fruitful discussions and helpful comments. This work was partially supported by ANR program COMMA (http://www-lmc.imag.fr/COMMA/).


domain: info scco

text:
Baptiste Rohou-Claquin

Malika Auvray

Jean-Loïc Le Carrou

David Gueorguiev

Perception of guitar strings on a flat visuo-haptic display ⋆

Keywords: Vibrotactile rendering, Guitar string discrimination, Visuohaptic perception

There is a rapid growth of interest in digital musical instruments due to advantages they present over their physical counterparts such as portability and novel experiences. However, these instruments still rarely provide haptic feedback, which results in weakened user experience. This study investigates whether haptic renderings of three guitar strings are discriminated by touch and whether they provide a realistic sensation of plucking a guitar string. Specifically, five methods to record and replay the vibrations induced by the plucking of a guitar string were tested. Three procedures relied on recordings made on the right hand's index finger, which underwent three types of filtering. Another one was recorded on the left hand kept motionless on the guitar's fingerboard. The fifth method recorded the vibrations directly on the neck of the guitar. The results revealed that good discrimination of the rendered string occurred only for recordings made on the motionless left hand's index finger or on the guitar. Participants also rated these conditions as being the most realistic, regardless of whether they were acquainted with musical instruments. Overall, this study shows that distinct guitar strings can be rendered by vibrotactile feedback and that tactile noise due to the interaction impairs their recognition and perceived realism.

Introduction

The field of musical instruments becomes increasingly digitalized, creating a category of instruments called digital musical instruments (DMI). This new category of instruments brings several advantages in comparison to their physical counterparts by overcoming the physical limitations imposed on tangible instruments [START_REF] Smith | Physical Audio Signal Processing[END_REF]. However, despite these advantages [START_REF] Serafin | Virtual reality musical instruments: State of the art, design principles, and future directions[END_REF][START_REF] Frid | Accessible digital musical instruments-a review of musical interfaces in inclusive music practice[END_REF], DMIs currently suffer from a large drawback compared to real instruments, which is the absence of haptic feedback. When a musician plays an instrument, a loop is created in the sense ⋆ Supported by UFR d'Ingénierie de Sorbonne Université and ANR Maptics of touch. The musician gives energy and perceives the instrument's response, creating a rich haptic exchange. The absence of compelling haptic feedback does not allow for this loop to be obtained in an optimal manner, thus diminishing the musician's ability to control the device and express itself. In addition, the performance and fidelity of the notes are reduced, and the listeners are also impacted by a less expressive performance [START_REF] Papetti | Musical Haptics[END_REF]. Indeed, the sense of touch is an essential sense for interacting with world around us. It enables us to explore the surrounding environment [START_REF] Lederman | Hand movements: A window into haptic object recognition[END_REF] and exhibits an astonishing sensitivity to subtle vibrations up to 1000 Hz [START_REF] Bolanowski | Four channels mediate the mechanical aspects of touch[END_REF] as well as small transient changes in the contact forces [START_REF] Gueorguiev | The tactile perception of transient changes in friction[END_REF].

Various studies have already investigated the impact of adding a haptic feedback on musical instruments [START_REF] Merchel | Tactile music instrument recognition for audio mixers[END_REF][START_REF] Turchet | Touching the audience: musical haptic wearables for augmented and participatory live music performances[END_REF][START_REF] Altınsoy | Touchscreens and Musical Interaction[END_REF][START_REF] Passalenti | No strings attached: Force and vibrotactile feedback in a virtual guitar simulation[END_REF] and conclude that such feedback can be helpful for musical perception. The addition of a haptic feedback is useful to improve the quality of play of the musician, but also, it increases the ease of learning, and the pleasure of practicing. Moreover, it facilitates access to people suffering from visual, auditory or deaf blindness impairments [START_REF] Turchet | Musical haptic wearables for synchronisation of visually-impaired performers: A co-design approach[END_REF][START_REF] Trivedi | Wearable musical haptic sleeves for people with hearing impairment[END_REF][START_REF] Nanayakkara | Enhancing musical experience for the hearing-impaired using visual and haptic displays[END_REF]. In order to feel these digital instruments, it is necessary to be able to perceive their textures and vibrations and numerous studies have shown that it is possible to discriminate textures, edges, shapes or vibrations that are generated by haptic devices [START_REF] Bau | Teslatouch: Electrovibration for touch surfaces[END_REF][START_REF] Fiedler | A Novel Texture Rendering Approach for Electrostatic Displays[END_REF][START_REF] Gueorguiev | The perception of ultrasonic square reductions of friction with variable sharpness and duration[END_REF]. Moreover, coupling haptic feedback with visual representation improves perceptual capacity [START_REF] Lacey | Visuo-haptic multisensory object recognition, categorization, and representation[END_REF][START_REF] Newell | Visuo-haptic perception of objects and scenes[END_REF], representation in memory [START_REF] Ernst | Merging the senses into a robust percept[END_REF], and spatiotemporal accuracy [START_REF] Altınsoy | Touchscreens and Musical Interaction[END_REF]. This study aims to investigate how different types of haptic rendering of a guitar string, which are coupled with an accurate visual rendering, are discriminated and perceived by users. Specifically, we investigate the capacity of users to recognize each digital string and their rating of the achieved realism in five experimental conditions.

Materials and methods

Data Collection

To be able to compare different methods for haptic rendering of guitar string's vibrations, we measured the vibrations elicited by the string on different locations:

-The tip of the right hand's index finger: The measure is performed while this hand plucks the string (Fig. 1.A). -The tip of the left hand's index finger: This finger pushes on the string at the level of the guitar's fingerboard and stays motionless during the recording (Fig. 1.B). -The guitar fingerboard: The measure was performed with the accelerometer placed directly on the wood of the guitar's neck (Fig. 1.C).

These measuring locations aimed to test whether signal components generated during the movement or during the interaction between the sensor and the nail impact the discrimination capacity of the user and the sensation of realism. The guitar was held in a natural way; the right hand plucks the string, the left hand holds the neck and pushes on the string. The measures were performed with a PCB352A21 accelerometer that has a sensitivity of 10 mV/g, a measure range of ±500 g pk, and a frequency range between 1.0 and 10000 Hz. The data were collected with a sampling frequency of 6250 Hz. For each position, three different notes were recorded corresponding to the vibration of three open strings : High (note E 4 ) at 329 Hz, Medium (note G 3 ) at 196 Hz and Low (note A 2 ) at 110 Hz. These three open strings were found sufficiently spaced in frequency to enable recognition significantly above chance but still challenging enough to prevent perfect performance. We implemented these frequency filters to preferentially activate the tactile channels respectively related to FA1 (<70 Hz) and FA2 (>70 Hz) tactile afferents with the aim to investigate whether filtering low-frequency tactile noise or the high-frequency peaks would impact perception.

Visual interface

In addition to the haptic feedback, a visual feedback of the strings' vibration was implemented on Java Processing (Fig. 3) following the equation of a vibrating string fixed at both ends, which was analytically computed:

y(x, t) = N p=1 A(p) sin pπ L x cos (ω p t) e -0.1ωpt (1)
With :

ω p = pπ L T µ (2)
A(p) = 2a p 2 π 2 1 α (1 -α) sin (pπα) (3)
For each note, Equation 1 was implemented in real time, except for a slight delay generated by the software (≃ 2ms). A represents the amplitude of the string in relation to the y-axis (width of the screen), α represents the ratio of the location where the string is plucked in comparison to the total length. Several mechanical parameters related to the string are also considered. L represents the total length of the string, T is the string tension, ω p is the angular frequency, which coincides with the harmonic frequency (f p ), and µ represents its linear mass. There are also some parameters specific to the simulation such as the sampling frequency, the harmonic number p, and the total number of harmonics N . We computed only the fundamental frequency for the visual rendering because we were limited by the refresh rate of the visual display (p=1, N=1, High=329 Hz, Medium=196 Hz and Low=110 Hz).

Perception study

To conduct the perception study, a vibrotactile actuator MM3C (Tactile Labs, Canada) was used. This actuator is specially designed to cover the frequency range perceived by the human finger. However, it distorts the signal at very low frequencies (< 30 Hz). In order to minimize disturbance during the test, a high-pass finite impulse response filter (FIR) at 30 Hz was used on each sample. No other distortion was observed during reproduction of the signals hence we did not apply additional low-pass filtering. In addition to the actuator, a 7-inch touch screen was used with a display frequency of 47 Hz. It enabled the participants to physically pluck the string, thus adding realism to the action.

Participants had to recognize three types of notes: high (E 4 ), medium (G 3 ) and low (A 2 ). For each of the five types of feedback, a total of 30 samples had to be recognized (10 per note type). The users place the index finger of their non-dominant hand on the actuator, which is located at the top of the screen (Fig. 4). Then, they pluck the string on the tablet with their dominant hand, as if they were playing on a real guitar (Fig. 4). After releasing the string, visual and haptic feedbacks are played simultaneously. Finally, after a few seconds, another screen appears that asks the participants which note was played: high, medium, or low. They choose one of them by pressing the corresponding button. The order of the blocks and of the notes within them are pseudo-randomized. In addition to recognizing the notes, participants are asked to rate the realism of the feedback they had just played in comparison with the impression of playing with a physical guitar. The score is reported using a scroll bar between 0 and 100%. Before each experimental block, which relates to a specific feedback type, participants familiarize themselves with the feedback by playing each note twice with knowledge of its type.

This study took place in an isolated space and participants wore ear muffs and earplugs in order to prevent auditory feedback. The study lasted between 25 and 30 minutes. In total, 21 people participated, 7 women and 14 men (mean age = 31.2, SD = 11.8). Participants also answered a questionnaire about their musical knowledge and experience with guitar playing. They rated their level of musical knowledge on a Likert-scale from 1 to 7 and they additionally reported whether it relates to guitar. 11 participants self-reported as guitarists among the 21 of the study. This question aimed to probe whether their musical knowledge influences their answers to the task. This study was conducted in accordance

Results

To investigate the impact of the feedback type on recognition of the played note and the sensation of realism, we performed one-way ANOVA statistical analyses after verification of the data by a d'Agostino-Pearson normality test. The ANOVA analysis related to recognition of high, medium, and low notes between the five feedbacks (Fig. 5.A) showed a significant effect of the haptic feedback type on the capacity of discrimination (F = 8.713, p = 0.0002). A post-hoc Tukey's multiple comparisons test was performed between individual feedbacks. Five post-hoc comparisons were found significant with p < 0.05: right hand vs. fingerboard, right hand<70 Hz vs. left hand, right hand<70 Hz vs. fingerboard, right hand>70 Hz vs. left hand, and right hand>70 Hz vs. fingerboard. Overall, significant comparisons occurred only between the feedbacks related to the right hand's recordings and those related to the left hand or fingerboard recordings. Another ANOVA was performed on the ratings about the feedback realism (Fig. 5.B). The results from the analyses showed a significant effect of the feedback type on the quality of the rating ( F=7.444, p = 0.0002). Five post-hoc comparisons were also found significant with p < 0.05: right hand vs. left hand, right hand vs. fingerboard, right hand<70 Hz vs. fingerboard, right hand>70 Hz vs. left hand, and right hand>70 Hz vs. fingerboard. Overall, the significant differences observed for the recognition task were similar to those observed for rating the realism. Generally, it appears that haptic feedback generated from recordings on the right hand, which is plucking the string, are less recognized and appreciated compared to the other feedback types.

The statistical test did not show significant differences between the experimental conditions related to the right hand either for note recognition or realism rating. Thus, filtering the signal recorded on the hand that plucks the guitar did not impact perception. The confusion matrices of the three main feedbacks were also computed showing similar trends. The confusion matrix when the feedback stems from recordings on the right hand (Fig. 6.A) has a success rate of about 50%, unlike the other two matrices that have success rates over 70% (Fig. 6.B & C). They also show that confusion mainly occurs between the high and medium notes, whose recognition is worse than for the low note regardless of the feedback type.

Discussion

The results of our study show that haptic feedback that stem from recordings on the hand that is plucking the string is significantly less recognized and appreciated than when the feedback is recorded on the the immobile left hand's index finger or on the guitar itself. Thus, the tactile noise generated by the kinesthetic interaction has a real impact on discrimination capacity, strongly reducing it. Moreover, users rate higher the realism when the recording is not impacted by kinesthetic motion (see Fig. 7). However, even when tactile noise distorts the signal and the signal is further filtered, performance is around 50% hence superior to chance level. This means that the users keep a certain capacity to discriminate visuo-haptic string rendering even if the recording is noisy or part of the frequency spectrum is missing due to filtering. Interestingly, the string vibration recorded directly on the guitar's neck hence without the finger in the recording loop was the highest rated by almost everyone. The absence of perturbation by the guitar plucking movement and of resonance damping by the fingertip have probably played a role in this condition being felt the most compelling. We can make a slight reservation concerning the similarities in the discrimination and rating abilities of the users between the three conditions that were designed by filtering the measurements on the right hand. The objective of the filtering differentiating the three feedbacks was to target the activation frequencies of the tactile channels related to FA1 and FA2 tactile afferents, respectively. However, the resulting signals still encompassed tactile noise around the cut-off frequency, which blurred the differences between the three conditions that used the recording from the right hand. However, this does not affect the conclusion that the feedbacks, in which tactile noise due to kinesthetic interaction is present, are significantly less recognized and are rated as less realistic. In spite of the fact that presentation of the experimental blocks (i.e. feedback types) and the order of the trials were randomized between participants, we still noticed a learning phenomenon during the experiment (Fig. 8). The recognition performance continuously improves with the progression of the blocks (Pearson's Correlation, p = 0.0193). This progression may be due at first to the discovery of the experiment, which might have required some time to master. However, it is astonishing that performance at recognizing the notes continued to improve after the first blocks had passed.

We also asked participants to rate their musical knowledge and whether it relates to guitar. However, we found no indication that musical experience or experience with guitar playing influenced their perception of the rendered guitar strings.

Conclusion

Our study shows that recognition of a digital guitar string varies greatly depending on how the signal is recorded. Users perform better with feedback without the tactile components related to string plucking, thus when sensing solely the instrument's response. Furthermore, the realism of string vibrations recorded on the motionless finger or directly on the guitar is consistently rated higher by participants with a further preference for the direct recording on the guitar.

Fig. 1 .

 1 Fig. 1. Measurement of the string's vibration with an accelerometer placed at several locations. A) On the right hand's index finger. B) On the left hand's index finger. C) On the guitar fingerboard

Fig. 2 .

 2 Fig. 2. Fast Fourier transform of the guitar string vibration signals, which are measured with the accelerometer. A) Right hand's signal. B) Left hand's signal. C) Guitar fingerboard's signal.

Fig. 3 .

 3 Fig. 3. Representation of a guitar string on the visuo-haptic display.

Fig. 4 .

 4 Fig. 4. The experimental set-up and the user's interaction with the simulated string

Fig. 5 .

 5 Fig. 5. A) Percentage of correct answers at the recognition task. The blue bars represent the median across participants and the error bars represent the lower and upper quartiles. B) Same plot for the rating of realism.

Fig. 6 .

 6 Fig. 6. Confusion matrix of the users' correct answers for the different feedback types. A) The right hand's feedback without filtering. B) The left hand's feedback. C) The guitar fingerboard's feedback

Fig. 7 .

 7 Fig.7. User preferences. In blue, the feedbacks that obtained the highest ratings by the users during the experiment. In orange, those whose rating was the 2 nd highest.

Fig. 8 .

 8 Fig. 8. A linear fit is performed on the average of individual performances across the experimental blocks in chronological order (Mean ± SD).