Sophie Thery
email: sophie.thery@univ-grenoble-alpes.fr

Charles Pelletier
email: charles.pelletier@uclouvain.be

Florian Lemarié
email: florian.lemarie@inria.fr

Eric Blayo
email: eric.blayo@univ-grenoble-alpes.fr

Analysis of Schwarz waveform relaxation for the coupled Ekman boundary layer problem with continuously variable coefficients

Keywords: Schwarz waveform relaxation, variable coefficients, multiphysics coupling, coupled Ekman problem, numerical climate modeling

In this paper we present a global-in-time non-overlapping Schwarz method applied to the Ekman boundary layer problem. Such a coupled problem is representative of large-scale atmospheric and oceanic flows in the vicinity of the air-sea interface. Schwarz waveform relaxation (SWR) algorithms provide attractive methods for ensuring a "tight coupling" between the ocean and the atmosphere. However the convergence study of such algorithms in this context raises a number of challenges. Numerous convergence studies of Schwarz methods have been carried out in idealized settings, but the underlying assumptions to make these studies tractable may prohibit them to be directly extended to the complexity of climate models. We illustrate this aspect on the coupled Ekman problem, which includes several essential

1 Introduction

Context and motivations

Owing to advances in computer resources and simulation tools, weather prediction systems now routinely include coupling mechanisms between several Earth system compartments; e.g. atmosphere, waves, ocean, sea-ice [START_REF] Brassington | Progress and challenges in short-to medium-range coupled prediction[END_REF][START_REF] Smith | Impact of coupling with an ice-ocean model on global medium-range NWP forecast skill[END_REF]. For such applications, the coupling methods used between individual components are usually ad-hoc algorithms raising synchrony issues [START_REF] Beljaars | On the numerical stability of surface-atmosphere coupling in weather and climate models[END_REF][START_REF] Lemarié | Analysis of ocean-atmosphere coupling algorithms: Consistency and stability[END_REF]. Following [START_REF] Keyes | Multiphysics simulations: Challenges and opportunities[END_REF], such methods are referred to as loose coupling methods, since it can be shown that the coupling is relaxed to one single iteration of what should be an iterative process. Due to the overwhelming complexity of the coupling mechanisms in the context of weather and climate, the consequences of inaccuracies in coupling algorithms on the stability and regularity of solutions are hard to untangle for realistic applications, unless a strongly coupled solution can be used as a reference. Schwarz waveform relaxation (SWR) algorithms are attractive iterative coupling methods for geophysical applications. Indeed, they are very weakly intrusive, can be enhanced by an appropriate choice of transmission conditions, are naturally able to handle nonconformities in time and space, and can be combined with time parallelization techniques [START_REF] Blayo | Optimized Schwarz waveform relaxation algorithms with nonconforming time discretization for coupling convection-diffusion problems with discontinuous coefficients[END_REF][START_REF] Gander | Optimized Schwarz waveform relaxation methods for advection reaction diffusion problems[END_REF][START_REF] Gander | Parareal Schwarz waveform relaxation methods[END_REF]. Several mathematical studies have investigated applying SWR algorithms to increase the efficiency of domain decomposition methods for increasingly complex problems derived from Earth sciences, e.g. for the so-called shallow-water equations [START_REF] Martin | Schwarz waveform relaxation algorithms for the linear viscous equatorial shallow water equations[END_REF][START_REF] Qaddouri | Optimized Schwarz methods with an overset grid for the shallow-water equations: preliminary results[END_REF], or for the ocean primitive equations [START_REF] Audusse | Optimized Schwarz waveform relaxation for the primitive equations of the ocean[END_REF]. At a more practical level, the dynamical core of at least one Earth system numerical model relies on a SWR algorithm for its domain decomposition method [START_REF] Qaddouri | The Canadian global environmental multiscale model on the Yin-Yang grid system[END_REF], and empirical studies have been undertaken for including Schwarz methods in ocean -atmosphere couplers [START_REF] Marti | A Schwarz iterative method to evaluate ocean-atmosphere coupling schemes. Implementation and diagnostics in IPSL-CM6-SW-VLR[END_REF]. Generally speaking, the theoretical convergence analysis of SWR algorithms is restricted to relatively simple linear model problems, often approximating space-varying quantities as constant coefficients, assuming that scalar variables are independent from each other, so that all the conditions to apply Fourier and Laplace transformations are met. As a matter of fact, coupled problems arising in Earth system modeling often exhibit sharp turbulent boundary layers near the interfaces whose effect is parameterized using a spatially variable viscosity coefficient [START_REF] Monin | Basic laws of turbulent mixing in the surface layer of the atmosphere[END_REF][START_REF] O'brien | A note on the vertical structure of the eddy exchange coefficient in the planetary boundary layer[END_REF]. The objective of this paper is to study the SWR algorithm for a model problem representative of the coupling between the turbulent oceanic and atmospheric surface boundary layers.

The coupled Ekman problem

We specifically study the so-called coupled Ekman layer problem [START_REF] Ekman | On the influence of the Earth's rotation in ocean-currents[END_REF], which has already been investigated by several physicists and mathematicians [e.g. [START_REF] Bye | Inertially coupled Ekman layers[END_REF][START_REF] Grisogono | A generalized Ekman layer profile with gradually varying eddy diffusivities[END_REF][START_REF] Lewis | Time-dependent, coupled, Ekman boundary layer solutions incorporating Stokes drift[END_REF]39] since it is representative of large-scale atmospheric and oceanic flows in the vicinity of the air-sea interface. So far, this type of coupled problems has been rarely studied with a focus on the numerical techniques to solve it, as proposed in the present work. The atmospheric and oceanic flows at large scale are generally represented by the so-called primitive equations, i.e. Navier-Stokes equations written in a rotated frame under the hydrostatic assumption. Assuming that the fluid is vertically homogeneous (i.e. time and space independent density) and that horizontal velocity gradients are small compared to vertical ones, the equation for the horizontal velocity U = (u, v) reduces to the Reynolds averaged primitive equations:

∂ t U + f k × U = -∂ z U w -ρ -1 ∇p + F U , (1)
where k is a vertical unit vector, U w is the vertical component of the Reynolds stress tensor, ρ is the density, p is the corresponding hydrostatic pressure, and F U contains any other forcing terms. f corresponds to the local Coriolis frequency associated with the Earth's rotation, whose values typically range from f = 0 s -1 at the Equator to f ≈ ±1.45 × 10 -4 s -1 at the poles. For the U w term, a common closure assumption (known as the Boussinesq hypothesis) is to consider that turbulence acts as a mixing term, thus leading to U w = -ν(z)∂ z U where ν is a so-called turbulent viscosity, much larger than the each fluid's molecular viscosities, and which varies with z [START_REF] O'brien | A note on the vertical structure of the eddy exchange coefficient in the planetary boundary layer[END_REF]. Introducing the stationary geostrophic velocity U g defined as -f k × U g = ρ -1 ∇p, (1) becomes

∂ t U + f k × U -∂ z (ν(z)∂ z U) = F U + f k × U g . (2)
Assuming that U g is known, the problem is one-dimensional and represents the balance between the Coriolis term and the stress due to small-scale turbulent motion, under the assumption that the large-scale dynamics are in geostrophic balance. This problem typically involves boundary layers (in this context, also known as Ekman layers) which have a typical depth of tens of meters in the ocean and hundreds of meters in the atmosphere.

In the following, the physical domain of interest is Ω =]h 1 , h 2 [ where h 1 < 0 < h 2 , and {z = 0} is the interface between two non-overlapping subdomains. The index j ∈ {1, 2} is used to distinguish subdomains, and for sake of simplicity the interval (a, b) stands for ]a, b[ if a < b and for ]b, a[ if a > b. We note U j = (u j , v j ) the Ekman components of the velocity in subdomain Ω j = (0, h j ), ν j (z) > 0 the turbulent viscosity, and the r.h.s. in [START_REF] Beljaars | On the numerical stability of surface-atmosphere coupling in weather and climate models[END_REF] is grouped as a single forcing term denoted F j .

On the time interval ]0, T ], the coupled Ekman boundary layer problem of interest is:

∂ t U j + f k × U j -∂ z (ν j (z)∂ z U j ) = F j , (z, t) ∈ Ω j ×]0, T ], j ∈ {1, 2}, (3a)
U j (z, 0) = U j,0 (z), z ∈ Ω j , j ∈ {1, 2},

U j (h j , t) = U j,g (t), t ∈ [0, T ], j ∈ {1, 2}, (3c)

ν 1 (0 -)∂ z U 1 (0 -, t) = ν 2 (0 + )∂ z U 2 (0 + , t), t ∈ [0, T ], (3d)
U 2 (0 + , t) = U 1 (0 -, t), t ∈ [0, T ], (3e)
where U j,0 (z) is the initial condition, and U j,g (t) is the geostrophic velocity component imposed at the external boundaries of the spatial domain.

Outline of the paper

When discretizing (3a) in numerical models, the turbulent viscosity is systematically treated implicitly in time due to the small vertical grid spacing (relative to the horizontal one) typically used in most configurations [START_REF] Lemarié | Stability constraints for oceanic numerical models: implications for the formulation of time and space discretizations[END_REF]. This implicit treatment yields significant challenges in solving the coupled problem at a numerical level. A possibility is to consider iterative SWR methods to solve (3), as motivated in Sec. 1.1.

The convergence properties of Schwarz algorithms with continuously variable coefficients in each subdomain have already been studied using the technique of separation of variables to decouple the space and time variables (or equivalently the two spatial variables for steady-state problems) directly in the physical space. However, previously obtained results are either restricted to cases where the coefficients vary in the direction parallel to the interface [START_REF] Gander | Optimized Schwarz methods for model problems with continuously variable coefficients[END_REF], or difficult to interpret because the resulting convergence factors are given by a numerical series of abstract eigenmodes and eigenvalues [START_REF] Lemarié | Toward an optimized global-in-time Schwarz algorithm for diffusion equations with discontinuous and spatially variable coefficients, part 2: the variable coefficients case[END_REF]. The seminal work of [28] showed, for the one-dimensional steady-state model problem -∂ z (ν(z)∂ z u) = F , that the absorbing boundary conditions depend on the full diffusion profile and not only on its local interfacial values. In the same spirit, the aim of the present paper is to further study the effect of continuously variable coefficients on the convergence properties of SWR algorithms. We study three specific expressions of the viscosity ν j (z), which have already been considered in the Ekman problem literature:

(i) the laminar case of [START_REF] Ekman | On the influence of the Earth's rotation in ocean-currents[END_REF] with constant viscosity (ν j (z) = ν j = cst);

(ii) as in [START_REF] Lewis | Time-dependent, coupled, Ekman boundary layer solutions incorporating Stokes drift[END_REF][START_REF] Madsen | A realistic model of the wind-induced Ekman boundary layer[END_REF], a more realistic case is studied by considering turbulent viscosities linearly increasing with depth (resp. with height);

(iii) as in [START_REF] Grisogono | A generalized Ekman layer profile with gradually varying eddy diffusivities[END_REF][START_REF] O'brien | A note on the vertical structure of the eddy exchange coefficient in the planetary boundary layer[END_REF], nearly parabolic shape functions for ν j (z) are introduced, consistently with what is used in climate simulations.

These studies addressed exclusively physical issues, and none of them tackled the numerical resolution of the associated coupled problem.

The paper is organized as follows. The well-posedness of the coupled Ekman problem (3) is studied in Sec. 2, leading to a SWR algorithm being formulated. Its convergence rates for various interface conditions and viscosity profiles are determined in Sec. 3. Their asymptotic behaviour w.r.t. the parameters of the problem are investigated in Sec. 4. Finally, numerical results illustrating the relevance of our theoretical analysis are shown in Sec. 5.

2 Formulation of the global-in-time Schwarz method and well-posedness

Well-posedness of the coupled problem

The coupled model (3) corresponds to a global parabolic problem with a discontinuous diffusion coefficient at the interface:

∂ t U + 0 -f f 0 U -∂ z (ν(z)∂ z U) = F, (z, t) ∈]h 1 , h 2 [×]0, T ], (4a)
U(z, 0) = U 0 (z), z ∈ [h 1 , h 2 ], (4b) U(h j , t) = U g (t), t ∈ [0, T ], j ∈ 1, 2. (4c)
In [START_REF] Lions | Problèmes aux limites non homogènes et applications[END_REF] (Chapter 4 and example 4.

7.1.), it is shown that if ν(z) > 0 is bounded in [h 1 , h 2 ], U 0 ∈ L 2 (Ω), U g ∈ H 3/4 ([0, T ]) and F ∈ L 2 (0, T ; L 2 (Ω)), then (4) has a unique weak solution in L 2 (0, T ; H 1 (Ω)) ∩ C 0 (0, T ; L 2 (Ω)
). This proof is based on an energy estimate.

The corresponding natural transmission conditions for (4) at z = 0 read:

ν 1 (0 -)∂ z U 1 (0 -, t) = ν 2 (0 + )∂ z U 2 (0 + , t), t ∈ [0, T ], (5a)
U 2 (0 + , t) = U 1 (0 -, t), t ∈ [0, T ], (5b)
where ν j = ν| Ωj . These are the interface conditions (3d) and (3e) of the original problem. ( 4) is thus equivalent to (3), hence it is well-posed in L 2 (0, T ; H 1 (Ω)) ∩ C 0 (0, T ; L 2 (Ω)).

Change of variables

In each medium (ocean and atmosphere), the evolution of the two velocity components, u j and v j , are coupled by the Coriolis term. In a matrix form, the systems (4a) for j ∈ {1, 2} can be rewritten as:

∂ t U j + AU j -∂ z (ν j (z) ∂ z U j ) = F j , with A = 0 -f f 0 . ( 6
)
In order to decouple the equations on ϕ 1 and ϕ 2 , we diagonalize A = PDP -1 , where:

D = if 0 0 -if , P = 1 √ 2 1 1 -i i and P -1 = 1 √ 2 1 i 1 -i .
Inserting this diagonalization and multiplying on the left by P -1 , (6) becomes:

∂ t Φ j + DΦ j -∂ z (ν j (z) ∂ z Φ j ) = F Φ j , with Φ j = P -1 U j = ϕ j (z, t) ϕ j (z, t) , (7)
which is equivalent to the following set of two independent scalar equations:

∂ t ϕ j + if ϕ j -∂ z (ν j (z) ∂ z ϕ j ) = F ϕ j , (8a)
∂ t ϕ j -if ϕ j -∂ z ν j (z) ∂ z ϕ j = F ϕ j = F ϕ j . (8b)
Rewritten with these new variables, the coupled Ekman layer problem (3) now reads:

∂ t Φ j + DΦ j -∂ z (ν j (z) ∂ z Φ j ) = F Φ j , (z, t) ∈ Ω j ×]0, T ] , j ∈ {1, 2}, (9a)
Φ j (z, 0) = Φ 0,j (z), z ∈ Ω j , j ∈ {1, 2}, (9b)
Φ j (h j , t) = Φ g,j (t), t ∈ [0, T ] , j ∈ {1, 2}, (9c)
Φ 1 (0 -, t) = Φ 2 (0 + , t), t ∈ [0, T ], (9d)
ν 1 (0 -) ∂ z Φ 1 (0 -, t) = ν 2 (0 + ) ∂ z Φ 2 (0 + , t), t ∈ [0, T ]. (9e)
Note that since (3) is well-posed, so is (9).

Schwarz Waveform Relaxation algorithms for the coupled Ekman problem

We consider that both submodels are independent from each other, thus prohibiting solving [START_REF] Ekman | On the influence of the Earth's rotation in ocean-currents[END_REF] with a single-model monolithic method. Let us introduce a SWR algorithm for solving the first component of [START_REF] Ekman | On the influence of the Earth's rotation in ocean-currents[END_REF], related to ϕ. The time interval [0, T ] is partitioned into N windows (T n ) 1≤n≤N , the iterative algorithm being applied successively on each window. For a given window T n , two sequences of subsolutions ϕ k j k∈N * (j ∈ {1, 2}) are defined by:

∂ t ϕ k j + ıf ϕ k j -∂ z ν j ∂ z ϕ k j = F ϕ j , (z, t) ∈ Ω j × T n , j ∈ {1, 2}, (10a)
ϕ k j (z, 0 Tn ) = ϕ 0,Tn j (z) , z ∈ Ω j , (10b)
ϕ k j (h j , t) = ϕ g,j (t), t ∈ T n , (10c)
B 1,1 ϕ k 1 0 -, t = B 1,2 ϕ k 2 0 + , t , t ∈ T n , (10d)
B 2,2 ϕ k 2 0 + , t = B 2,1 ϕ k-1 1 0 -, t , t ∈ T n . (10e)
In (10b), 0 Tn is the initial time in T n , and ϕ 0,Tn j are either given by (9b) (if n = 1), or by the final state of the solution on T n-1 (if n ≥ 2). In the interface conditions (10d) and (10e), B j,l are four interface operators to be chosen. The simplest choice would be B

1,1 = B 1,2 = Id, B 2,1 = ν 1 (0 -)∂ z and B 2,2 = ν 2 (0 + )∂ z ,
which, provided that the SWR algorithm converges, ensures that the coupling constraints (9d)-(9e) will be satisfied by the solution (ϕ 1 , ϕ 2 ). However, as discussed in Sec. 3, other choices are possible to accelerate the convergence of (10) toward the solution of ( 9). Here we study the "multiplicative" Schwarz method: first, the problem on ϕ k 2 is solved using ϕ k-1

1

; then, the problem on ϕ k 1 can be solved. Note that the well-posedness of the coupled problem, as well as SWR algorithm, applied to a similar problem with constant viscosity in each subdomain (but discontinuous at the interface) has been studied in [START_REF] Gander | A Schwarz waveform relaxation method for advectiondiffusion-reaction problems with discontinuous coefficients and non-matching grids[END_REF].

Well-posedness of the SWR algorithm

The coupling algorithm (10) is made of successive resolutions restricted to one subdomain (either (h 1 , 0) or (0, h 2 )). To ensure the well-posedness of the algorithm, it is first necessary to ensure the well-posedness of the subproblems, i.e. the existence and uniqueness of both subsolutions, and a sufficient regularity for these solutions to be applied as a boundary condition for the other subproblem. This leads to study the generic problem:

∂ t ϕ(z, t) + if ϕ(z, t) -∂ z (ν(z)∂ z ϕ(z, t)) = F ϕ (z, t), (z, t) ∈ (0, h)×]0, T ], (11a)
ϕ(z, 0) = ϕ 0 (z), z ∈ (0, h), (11b)
ϕ(h, t) = ϕ g (t), t ∈ [0, T ], (11c) Bϕ(0, t) = G(t), t ∈ [0, T ], (11d)
where B is a boundary operator. In the same way as in Sec. 2.1, the Coriolis effect has no impact on the well-posedness. Thus, according to [27], if the following conditions are satisfied:

(i) ν is bounded and strictly positive on (0, h);

(ii) F ϕ ∈ L 2 (0, T, L 2 (0, h));

(iii) ϕ g (t = 0) = ϕ 0 (z = h) and ϕ 0 (t = 0) = G(z = 0, t = 0);

(iv) ϕ 0 ∈ H 1 (0, h); (v) ϕ g ∈ H 3/4 (0, T ); (vi) B is a Dirichlet (B = Id) or a Neumann (B = ∂ z ) operator; (vii) G ∈ H 3/4 (0, T ) if B is a Dirichlet operator or G ∈ H 1/4 (0, T ) if B is a Neumann operator;
then there exists a unique solution of (11) in L 2 (0, T, H 2 (0, h))∩H 1 (0, T, L 2 (0, h)).

Both subproblems of the (10) coupled problem are thus well-posed. Moreover, the conditions (i)-(vii) listed above also guarantee that the trace Bϕ(0, t) satisfies both regularity conditions required for G in (vii). Therefore, the Schwarz algorithm on coupled problem (10) is well-posed. If B is a Robin operator (B = ∂ z + rId), the same result can be obtained by introducing further constraints (see Sec. 3.2.2). While less constraining regularity conditions guaranteeing the well-posedness of ( 11) exist (depending on B), from now on we will assume the conditions (i)-(vii) to be satisfied, since they offer enough flexibility given our context.

Analytical determination of convergence rates

Setting the (B j,l ) operators directly from (9d) and (9e) would mean using the physical coupling operators. This has been shown to be at best suboptimal (slow SWR convergence), at worst leading to non-converging algorithms, even with both subproblems being well-posed (e.g., in the case of non-overlapping domains [START_REF] Gander | Optimal convergence for overlapping and non-overlapping Schwarz waveform relaxation[END_REF]). In the following, we derive convergence rates for two different sets of boundary operators.

3.1 General form of convergence rates

Ordinary differential equation for Schwarz iterate errors

A usual approach to prove the convergence of SWR algorithms is to analytically determine e k ϕ,j = ϕ k jϕ j , the errors between SWR iterates and the exact solution ϕ j . In a given time interval [0, T ], these errors satisfy the Ekman problem (10) with null r.h.s., null initial conditions, and null boundary conditions at z = h j . Results introduced in Sec. 2.4 guarantee that e k ϕ,j has L 2 space regularity in [0, T ]. A continuous time Fourier transform can be performed, assuming that T → ∞ and extending all functions to zero for t < 0. The problem satisfied by the errors e k ϕ,j (z, ω) in Fourier space is:

ı(f + ω) e k ϕ,j -∂ z ν j (z)∂ z e k ϕ,j = 0, (z, ω) ∈ Ω j × R, (12a) e k ϕ,j (h j , ω) = 0, ω ∈ R, (12b) B 1,1 e k ϕ,1 0 -, ω = B 1,2 e k ϕ,2 0 + , ω , ω ∈ R, (12c) B 2,2 e k ϕ,2 0 + , ω = B 2,1 e k-1 ϕ,1 0 -, ω , ω ∈ R. (12d)
For each ω ∈ R, (12) is a set of two linear second-order ordinary differential equations on z for e k ϕ,1 , e k ϕ,2 . The differential operators in (12a) being independent of k, in each medium (j = 1, 2), there exists two basis functions f ν,j (z, ω) and g ν,j (z, ω) such that:

e k ϕ,j (z, ω) = A k ν,j (ω) f ν,j (z, ω) + B k ν,j (ω) g ν,j (z, ω), (z, ω) ∈ Ω j × R, (13)
where A k ν,j , B k ν,j is determined thanks to the boundary conditions (12b) and transmission conditions (12c)-(12d). In Sec. 3.3, explicit forms of the basis functions are given for the three physically meaningful forms of ν j introduced at the end of Sec. 1.3.

Convergence rates in the Fourier and the physical spaces

Convergence rates are the usual metric for assessing the asymptotic behavior of a SWR algorithm. They can be defined either in the Fourier space or in the physical space. Here we proceed as usually done with SWR algorithms: first, we define the convergence rate in the Fourier space for ϕ; second, we investigate its relation to the convergence rate in the physical space for U, which is our final interest.

In the Fourier space, the convergence rate of the SWR algorithm (10) can be defined as:

ρ k ϕ,j (ω) = e k ϕ,j (z = 0 ∓ , ω) e k-1 ϕ,j (z = 0 ∓ , ω) , ω ∈ R, k ∈ N * , j ∈ {1, 2}, (14)
where ∓ depends on j ∈ {1, 2}. ρ k ϕ,j (ω) represents the damping (or the amplification) of the error at the k th SWR iteration in the Fourier space, at the frequency ω. A sufficient condition for the convergence of the SWR algorithm is:

sup ω∈R ρ k ϕ,j (ω) < 1, ∀k ∈ N * , ∀j ∈ {1, 2}. (15)
The smaller ρ k ϕ,j , the faster convergence is reached. Condition [START_REF] Gander | Parareal Schwarz waveform relaxation methods[END_REF] is formulated in the Fourier space, yet our main interest is convergence in the physical space. Since t → e k ϕ,j (0 ∓ , t) and ω → e k ϕ,j (0 ∓ , ω) are L 2 functions, Parseval's theorem yields:

inf ω∈R ρ k ϕ,j (ω) ≤ e k ϕ,j (0 ∓ , t) 2 e k-1 ϕ,j (0 ∓ , t) 2 = e k ϕ,j (0 ∓ , ω) 2 e k-1 ϕ,j (0 ∓ , ω) 2 ≤ sup ω∈R ρ k ϕ,j (ω) . (16)
(15) is therefore a sufficient condition for (10) to converge in the L 2 physical space as well. In practice, we are interested in the numerical solutions of the time-dependent problem (3), using the physical variables (U 1 , U 2 ). The numerically observed error of this problem is:

E k,m j (z) = U k,m j (z) -U j (z, t m j ), 1 ≤ m ≤ M j , j ∈ {1, 2}, (17)
where m is a time step index, t m j is the physical time corresponding to m, and M j is the total number of time steps for the numerical approximation of U j in a given time window T n . The associated convergence rate of the SWR algorithm in the physical space is:

R k Uj = Mj m=1 E k,m j (0 ∓ ) 2 Mj m=1 E k-1,m j (0 ∓ ) 2 , k ∈ N * , j ∈ {1, 2}. (18)
Using the change of variables [START_REF] Dubois | The optimized Schwarz method with a coarse grid correction[END_REF], it is shown in App. 2 that:

R k Uj ≤ max    sup |ω|≤π/ min j∈{1,2} {δtj } ρ k ϕ,j (ω), sup |ω|≤π/ min j∈{1,2} {δtj } ρ k ϕ,j (ω)    + , (19)
where k ∈ N * , j ∈ {1, 2} and -→ Mj →+∞ 0.

(19) bounds R k Uj , the SWR algorithm convergence rate in the physical space, with its Fourier-space translation ρ k ϕ,j , which is easier to assess. Inequality [START_REF] Grisogono | A generalized Ekman layer profile with gradually varying eddy diffusivities[END_REF] will be used for the numerical results presented in Sec. 5. In the remainder, the interval of relevant values for ω from [START_REF] Grisogono | A generalized Ekman layer profile with gradually varying eddy diffusivities[END_REF], which depends on the time discretization, is:

I ω :=   - π min j∈{1,2} {δt j } ; π min j∈{1,2} {δt j }   = [-ω max ; ω max ]. ( 20
)
Remark 1 The upper bound [START_REF] Grisogono | A generalized Ekman layer profile with gradually varying eddy diffusivities[END_REF] does not account for the partitioning in time windows. In [START_REF] Gander | Dirichlet-Neumann and Neumann-Neumann waveform relaxation for the wave equation[END_REF], a link is established between the length of the time windows and the L ∞ convergence. However this result is valid for a particular case of interface conditions with a relaxation parameter, which is not directly translatable to our case. Indeed, the convergence factors arising from (10) do not benefit from the same regularity properties as those from [START_REF] Gander | Dirichlet-Neumann and Neumann-Neumann waveform relaxation for the wave equation[END_REF].

Convergence rates for different interface conditions

We now provide exact expressions of convergence rates [START_REF] Gander | A Schwarz waveform relaxation method for advectiondiffusion-reaction problems with discontinuous coefficients and non-matching grids[END_REF] for two different choices of interface operators (B j,l ).

Dirichlet-Neumann interface conditions

Implementing a SWR algorithm with the natural Dirichlet-Neumann (DN) transmission conditions (9d) and (9e) yields choosing B 1,1 = B 1,2 = Id, and B 2,j = ν j ∂ z .

The determination of the (A k ν,j , B k ν,j ) coefficients can be done by injecting the general form ( 13) for e k ϕ,j into boundary conditions (12c) and (12d). The following convergence rate, independent of both k ∈ N and j ∈ {1, 2}, can then be obtained:

ρ k,DN ϕ,j (ω) = ρ DN ϕ (ω) = λ |S 1 (ω) S 2 (ω)| , ω ∈ R (21)
where:

λ = ν 1 (0 -) ν 2 (0 + ) , (22a)
S 1 (ω) = f 1 (h 1 , ω) ∂ z g 1 (0, ω) -g 1 (h 1 , ω) ∂ z f 1 (0, ω) f 1 (h 1 , ω) g 1 (0, ω) -g 1 (h 1 , ω) f 1 (0, ω) , (22b)
S 2 (ω) = f 2 (h 2 , ω) g 2 (0, ω) -g 2 (h 2 , ω) f 2 (0, ω) f 2 (h 2 , ω) ∂ z g 2 (0, ω) -g 2 (h 2 , ω) ∂ z f 2 (0, ω) . (22c)
In the next subsections, we will give the exact expressions of S 1 and S 2 for different forms of ν j . If ν 1 and ν 2 are of similar nature (i.e. both constant, linear or parabolic), then S 1 = ψ(h 1 , ν 1 ) and S 2 = 1/ψ(h 2 , ν 2 ), where ψ depends on the particular form of ν j .

Robin-Robin interface conditions

We now consider the Robin-Robin (RR) transmission conditions:

B 1,1 = Id + qλ∂ z , B 1,2 = Id + q∂ z , B 2,2 = p Id + ∂ z , B 2,1 = p Id + λ∂ z , (23)
where (p, q) ∈ C 2 are free parameters and λ is defined in (22a). In order to be compatible with [START_REF] Gander | Optimized Schwarz methods for model problems with continuously variable coefficients[END_REF], the Robin-Robin B j,l operators must be in the form of ( 23), since any other choice of RR operators would not provide the same coupled solution as the natural DN transmission conditions (9d)-(9e).

Proposition 1 If:

(p, q) ∈ P = (p, q) ∈ C 2 Re(p) ≤ 0 and Re(q) ≥ 0 , (24)
then SWR algorithm (10), using Robin-Robin transmission conditions as defined in [START_REF] Lemarié | Analysis of ocean-atmosphere coupling algorithms: Consistency and stability[END_REF], is well-posed on each T n window.

The proof (not detailed here) is based on a priori energy estimates, can be directly derived from [START_REF] Bennequin | Optimized Schwarz Waveform Relaxation Methods for Convection Reaction Diffusion Problems[END_REF], and is available in [START_REF] Thery | Etude numérique des algorithmes de couplage océan-atmosphère avec prise en compte des paramétrisations physiques des couches limites[END_REF]. From now on, we will assume that (p, q) ∈ P, so that the SWR algorithm [START_REF] Gander | Optimized Schwarz methods for model problems with continuously variable coefficients[END_REF] with RR conditions is well-posed.

Similar to the DN case, the A k ν,j , B k ν,j decomposition coefficients are determined using the boundary conditions (12b), (12c) and (12d). The resulting convergence rate for Robin-Robin conditions, written in terms of the S j functions defined in [START_REF] Lemarié | Toward an optimized global-in-time Schwarz algorithm for diffusion equations with discontinuous and spatially variable coefficients, part 2: the variable coefficients case[END_REF], reads:

ρ RR ϕ (p, q, ω) = (S 2 + q)(p + λS 1 ) (1 + pS 2 )(1 + λqS 1 ) = ρ DN ϕ (ω) (1 + qS -1 2 )(1 + λ -1 pS -1 1 ) (1 + pS 2 )(1 + λqS 1 )
.

(25) As expected, setting p = q = 0 in (25) corresponds to the DN case, i.e. ρ DN ϕ (ω) = ρ RR ϕ (0, 0, ω).

Analytical expressions of the solution for various forms of ν

We now give analytical expressions of the solutions to [START_REF] Gander | Optimized Schwarz methods for model problems with continuously variable coefficients[END_REF] for choices of ν(z) (constant, affine, and parabolic) commonly used in Earth system applications. We aim at investigating and comparing ρ DN (ω) and ρ RR (p, q, ω), for such forms of ν j (z) assumed strictly positive and bounded in Ω j .

Remark 2 As explained in Sec. 2.2, convergence rates should be studied for both ϕ and its complex conjugate ϕ. However, since both subequations of ( 8) are conjugation of one another, Fourier transform properties yield:

ρ DN ϕ (ω) = ρ DN ϕ (-ω), (26a)
ρ RR ϕ (p, q, ω) = ρ RR ϕ (p, q, -ω). (26b)
Since I ω , the interval of interest for ω, is symmetric w.r.t. 0 (see ( 20)), we restrict our study to ϕ, with results on ϕ arising from (26).

Constant viscosity

In the case ν j (z) = ν c j , the basis functions (f j , g j ) are:

f cst j (z) = exp χ f + ω ν c j z and g cst j (z) = exp -χ f + ω ν c j z , (27)
where χ = e ı π 4 . Injecting ( 27) into ( 22) yields:

S cst 1 = -χ √ Fo 1 |h 1 | coth -χ Fo 1 , S cst 2 = -χ -1 |h 2 | √ Fo 2 tanh χ Fo 2 , (28)
where the nondimensional number Fo j = (f +ω)h 2 j ν c j has been introduced.

Affine viscosity profile

We now consider

ν j (z) = ν c j + b j z with the constraints b 1 < 0, b 2 > 0 leading to ν j (z) > 0 in Ω j .
The basis functions (f j , g j ) are:

f aff j (z) = I 0 2χ f + ω b 2 j ν j (z) and g aff j (z) = K 0 2χ f + ω b 2 j ν j (z) , (29)
where I 0 and K 0 are the modified 0 th -order Bessel functions. Injecting (29) into [START_REF] Lemarié | Toward an optimized global-in-time Schwarz algorithm for diffusion equations with discontinuous and spatially variable coefficients, part 2: the variable coefficients case[END_REF] yields:

S aff 1 = χ √ Fo 1 |h 1 | I 0 (2χ Fo 1 µ 1 (µ 1 + 1)) K 1 (2χ √ Fo 1 µ 1 ) + K 0 (2χ Fo 1 µ 1 (µ 1 + 1)) I 1 (2χ √ Fo 1 µ 1 ) I 0 (2χ Fo 1 µ 1 (µ 1 + 1)) K 0 (2χ √ Fo 1 µ 1 ) -I 0 (2χ √ Fo 1 µ 1 ) K 0 (2χ Fo 1 µ 1 (µ 1 + 1)) , (30a)
S aff 2 = -χ -1 |h 2 | √ Fo 2 I 0 (2χ Fo 2 µ 2 (µ 2 + 1)) K 0 (2χ √ Fo 2 µ 2 ) -I 0 (2χ √ Fo 2 µ 2 ) K 0 (2χ Fo 2 µ 2 (µ 2 + 1)) I 0 (2χ Fo 2 µ 2 (µ 2 + 1)) K 1 (2χ √ Fo 2 µ 2 ) + K 0 (2χ Fo 2 µ 2 (µ 2 + 1)) I 1 (2χ √ Fo 2 µ 2 ) , (30b)
where µ j is the following positive nondimensional number:

µ j = ν c j ∂ z ν j (0)h j = ν c j b j h j . (31)

Quadratic viscosity profile

Let us finally consider ν j (z) = ν c j + b j z + a j z 2 . To be consistent with the oceanatmosphere context, we impose a j < 0, ν j (h j ) = ν j (0) and ν j (z) > 0 in Ω j , which implies ν j (z) = ν c ja j h j z + a j z 2 . The basis functions (f j , g j ) are:

f par j (z) = P 0 ξj (η j (z)) = 2 F 1 ξ j + 1, -ξ j , 1, 1 2 - 1 2 η j (z) , (32a)
g par j (z) = P 0 ξj (-η j (z)) = 2 F 1 ξ j + 1, -ξ j , 1, 1 2 + 1 2 η j (z) , (32b)
where P 0 ξj are Legendre's functions and 2 F 1 (a, b, c, •) is the hypergeometric function [START_REF] Olver | Hypergeometric function[END_REF]. (32) also introduces ξ j ∈ C and η j (z) ∈ C for j ∈ {1, 2}, which are given in App. 1. Since we assume a j < 0 and ν j (z) > 0 in Ω j , ξ j and η j can be simplified to:

ξ j (ω) = - 1 2 1 ± 1 + 4ı (f + ω) a j = - 1 2 1 + 1 -4ıFo µ j , (33a)
η j (z) = ∂ z ν j (z) (a j h j ) 2 -4a j ν c j , (33b)
where µ j is defined like in [START_REF] Martin | Schwarz waveform relaxation algorithms for the linear viscous equatorial shallow water equations[END_REF]:

µ j = ν c j ∂ z ν(0)h j = ν c j |a j |h 2 j
. Computing the derivatives of the f par j and g par j functions defined in (32) yields:

S par 1 = (ξ 1 + 1)ξ 1 |h 1 | √ 1 + 4µ 1 P 0 ξ1 (η 1 )G ξ1 (η 1 ) + P 0 ξ1 (-η 1 )G ξ1 (-η 1 ) P 0 ξ1 (-η 1 ) 2 -P 0 ξ1 (η 1 ) 2 , (34a)
S par 2 = |h 2 | √ 1 + 4µ 2 (ξ 2 + 1)ξ 2 P 0 ξ2 (-η 2 ) 2 -P 0 ξ2 (η 2 ) 2 P 0 ξ2 (η 2 )G ξ2 (η 2 ) + P 0 ξ2 (-η 2 )G ξ2 (-η 2 ) , (34b)
where, for j ∈ {1, 2}, η j = η j (0) = (-1) j 1 + 4µ j = -η j (H j ), and

G ξj (±η j ) = 2 F 1 ξ j + 2, 1 -ξ j , 2, 1 2 - 1 2 (±η j ) . (35)
At this point we have prepared all the necessary ingredients to proceed to the actual convergence analysis for these three particular forms of ν j .

Convergence analysis

In this section, we study the behaviour of the convergence rates introduced in Sec. 3. Sec. 4.1 and 4.2 offer a review of known results on those convergence rates in specific cases, with the consistency of our results being checked in those limit cases. Sec. 4.3 focuses on the DN case and describes the derivation of upper bounds for the convergence factor in the linear and parabolic cases. Finally, Sec. 4.4 deals with the acceleration of SWR algorithm convergence using RR interface conditions.

Overview of the results obtained in simplified cases

Existing literature has already covered different results corresponding to particular cases of the more general problem (3). Here we described three such cases and express the corresponding S j quantities of interest. As emphasized by Sec. 5 of [28], if one considers the equation ∂ z (ν(z)∂ z ϕ) = 0, it is possible to guarantee convergence of the SWR algorithm in two iterations. This result can easily be extended to the case where ν(z) is discontinuous at the interface.

Using our notations, we get:

S sta 1,f=0 = 1 ν c 1 0 h1 (ν 1 (z)) -1 dz -1 , S sta 2,f=0 = -ν c 2 h2 0 (ν 2 (z)) -1 dz , (36)
with S sta j,f=0 = S j (ω = 0, f = 0, ν j (z)). The DN convergence factor thus reads ρ sta DN,f=0 = λ S sta 1,f=0 S sta 1,f=0 and the boundary conditions leading to convergence in two iterations (i.e., transparent boundary conditions) are conditions [START_REF] Lemarié | Analysis of ocean-atmosphere coupling algorithms: Consistency and stability[END_REF] with p = -λ S sta 1,f=0 and q = -S sta 2,f=0 . For our three particular viscosity profiles of interest, we easily find that: S sta,par 1,f=0 ) as functions of µ 2 (resp. µ 1 ) are displayed in Fig. 1a. Since the viscosity profiles tend to a constant value as µ j → +∞, the ratios tend to 1 for large values of µ j . However, for values of µ j smaller than 1, a significant impact of the variations in the viscosity profile on the values of S sta j can be observed. This would affect the convergence rate in the DN case as well as the optimal coefficients in the RR one.

S sta,cst 1,f=0 = 1 |h 1 | , S sta,cst 2,f=0 = -h 2 , S sta,aff 1,f=0 = 1 |h 1 | 1 µ 1 ln(1 + 1/µ 1 ) , S sta,aff 2,f=0 = -h 2 µ 2 ln(1 + 1/µ 2 ), S sta,par 1,f=0 = 1 |h 1 | √ 1 + 4µ 1 4µ 1 arccoth √ 1 + 4µ 1 , S sta,par 2,f=0 = -h 2 4µ 2 arccoth √ 1 + 4µ 2 √ 1 + 4µ 2 . ( 37
10 -5 10 -3 10 -1 10 1 10 3 µ 10 -3 10 -2 10 -1 10 0 µln 1 + 1 µ 4µ √ 1+4µ arcoth √ 1 + 4µ (

Stationary case with rotation and constant viscosity

(ω = 0, ν j = ν c j )
In the stationary (∂ t = 0) and constant viscosity case (ν j (z) = ν c j ), the values of S sta,cst j can easily be derived:

S sta,cst 1 = 1 |h 1 | χ Fo 0 1 tanh χ Fo 0 1 , S sta,cst 2 = -h 2 tanh χ Fo 0 2 χ Fo 0 2 , (38)
where Fo 0 j = f h 2 j /ν c j is the non-dimensional number introduced in (28), evaluated here for ω = 0. As shown in Fig. 1b, in the small Fo 0 j limit, the rotation's impact is negligible. The same remark withholds when Fo 0 1 ≈ Fo 0 2 with DN interface conditions, since then S sta,cst 1 S sta,cst 2 ≈ S sta,cst 1,f =0 S sta,cst 2,f =0 . In the general case, the rotation can become important in the convergence analysis as soon as Fo j > 1. This numerical conjecture implies that the relative importance of the Coriolis effect f , the fluid's viscosity ν c j , and the domain height h j , determines the convergence (or lack thereof) of SWR algorithm in this stationary case. In other words, in this case, SWR convergence is tributary to the particular setting of the physical problem: for example, convergence may be observed close to the equator (small f and Fo), but not at higher latitudes (large f and Fo).

(f = 0, ν j (z) = ν c j , |h j | → ∞)
This case is thoroughly studied in [10] for ν c 1 = ν c 2 , and in [START_REF] Lemarié | Toward an optimized global-in-time Schwarz algorithm for diffusion equations with discontinuous and spatially variable coefficients, part 1: The constant coefficients case[END_REF] for

ν c 1 = ν c 2 . With DN conditions, ρ ∞,cst DN = ν c 1 /ν c 2 = √ λ.
The transparent conditions cannot be determined locally in time (i.e., in the physical space), but it is possible to optimize the (p, q) parameters in order to minimize the convergence rate over a given range of temporal frequencies ω ∈ [ω min , ω max ]. Such an optimization problem is analytically solved in [START_REF] Lemarié | Toward an optimized global-in-time Schwarz algorithm for diffusion equations with discontinuous and spatially variable coefficients, part 1: The constant coefficients case[END_REF]. In Sec. 5, we use those optimized transmission conditions obtained for constant viscosity in the case of non-constant viscosity profiles, and investigate their efficiency.

Consistency between the different viscosity profiles

We now promptly check that the convergence rates derived in Sec. 3 encompass the limit cases discussed in Sec. 4.1. In particular, we consider (ω + f ) → 0, µ j → ∞, and |h j | → ∞. Including the rotation (i.e. f = 0) makes the Fo = 0 case important, since it corresponds to ω = -f , which is to be encountered in numerical climate applications. In the non-rotating case (f = 0), Fo = 0 is not encountered in practice, since it only occurs when ω = 0, and the stationary part of the SWR algorithm error stays at its initial state, which is zero.

Stationary limit without rotation (Fo j → 0)

In the constant viscosity case, the S cst j functions defined in (28) are such that lim Foj →0

S cst j = S sta,cst j,f =0 . For linear viscosity profiles, using Taylor expansion of Bessel functions yields:

lim Fo1→0 S aff 1 = 1 |h 1 | µ 1 ln (1 + 1/µ 1 ) = S sta,aff 1,f =0 , (39a)
lim Fo2→0 S aff 2 = -h 2 µ 2 ln (1 + 1/µ 2 ) = S sta,aff 2,f =0 , (39b)
which also holds in the parabolic case.

Similarly, assuming |h j | → ∞, the known results obtained under the infinite domain assumption can be obtained, both with and without rotation. It is worth mentioning that to maintain the consistency between the stationary and nonstationary cases for the convergence rates, we need to consider subdomains of finite sizes. Indeed, considering subdomains of infinite size implies that h 2 j ν c j /(f + ω), which can not hold as f + ω → 0.

Small viscosity gradient limit (µ j → ∞)

Since the parameters µ j are defined as µ j = νj (0)

hj ∂zνj (0) , the limit case ∂ z ν j (0) → 0 (i.e. µ j → ∞) should lead to results similar to those obtained in the constant viscosity case. Indeed, studying S par j and S aff j leads to:

lim µj →∞ S par j = lim µj →∞ S aff j = S cst j . ( 40
)
This result implies that for both DN and RR transmission conditions, in the µ j → ∞ limit, the convergence rates become equivalent for constant, affine and parabolic ν j .

Behaviour of the convergence rates with DN transmission conditions

At this point we have checked that the convergence rates have the proper stationary (i.e. small Fo j ) limit and small viscosity gradients (i.e. µ j → ∞) limit. The objective is now to characterize the parameter values for which the effect of the continuous viscosity variations on the convergence rate can be significant in the DN case. [START_REF] Lemarié | Toward an optimized global-in-time Schwarz algorithm for diffusion equations with discontinuous and spatially variable coefficients, part 2: the variable coefficients case[END_REF] suggested that the convergence is insensitive to variations of ν j (z) at large frequencies. Indeed, in the ω → ∞ limit, regardless of the viscosity profile, the S 1 S 2 product, with S j defined in either (28), [START_REF] Marti | A Schwarz iterative method to evaluate ocean-atmosphere coupling schemes. Implementation and diagnostics in IPSL-CM6-SW-VLR[END_REF], or [START_REF] Olver | Hypergeometric function[END_REF], converges towards the same value

ν c 2 /ν c 1 = λ -1/2 .
The high frequency asymptote of the DN convergence rate is thus the same as the one obtained with constant coefficients and under the infinite domain assumption.

On the contrary, the stationary case discussed in Sec. 4.1 suggests that variations of the viscosity can have a large impact on the convergence rate at lower frequencies. Since the DN convergence rate is λ|S 1 S 2 |, the variations in the viscosity profiles on the Dirichlet side (i.e. in (0, h 1 )) and on the Neumann side (i.e. in (0, h 2 )) will have opposite effects in terms of convergence speed. To further investigate this aspect, we introduce the following two quantities:

R aff,cst j = S aff j S cst j , R par,cst j = S par j S cst j , (41)
which are functions of the dimensionless parameters Fo j and µ j . In the stationary case without rotation, we had R aff,cst ary case, Fig. 2 suggests the following conjectures:

1 = 1/[µ 1 ln (1 + 1/µ 1 )] and R par,cst 1 = √ 1 + 4µ 1 4µ 1 arccoth √ 1 + 4µ 1 -1 . R aff,cst
1 ≤ R aff,cst 1 ≤ lim Fo1→0 R aff,cst 1 = S sta,aff 1,f =0 S sta,cst 1,f =0 = 1 µ 1 ln (1 + 1/µ 1 ) , ( 42a
)
1 ≤ R aff,cst 1 ≤ lim Fo1→0 R aff,cst 1 = S sta,par 1,f =0 S sta,cst 1,f =0 = √ 1 + 4µ 1 4µ 1 arccoth √ 1 + 4µ 1 . (42b)
In the stationary case without Coriolis effect, from [START_REF] Qaddouri | Optimized Schwarz methods with an overset grid for the shallow-water equations: preliminary results[END_REF], and assuming that min z ν var = ν var (0) = ν cst , it can be shown that S sta,var 1,f =0 > S sta,cst 1,f =0 . In the nonstationary case and with Coriolis effect, (42) can only be numerically checked and is assumed true from now on. Upper and lower bounds of R aff,cst 2 (resp. R par,cst ). Moreover:

ρ ν DN = ρ cst DN R ν,cst 1 R ν,cst 2 , (43)
which, combined to (42), leads to:

S sta,aff 2,f =0 S sta,cst 2,f =0 ρ cst DN ≤ ρ aff DN ≤ S sta,aff 1,f =0 S sta,cst 1,f =0 ρ cst DN (44a) S sta,par 2,f =0 S sta,cst 2,f =0 ρ cst DN ≤ ρ par DN ≤ S sta,par 1,f =0 S sta,cst 1,f =0 ρ cst DN (44b)
with S sta,aff 1,f =0 / S sta,cst 1,f =0 and S sta,par 1,f =0 / S sta,par 1,f =0 given by (42). ( 44) contains upper bounds for the convergence rate with DN interface conditions in both the linear and parabolic viscosity cases. (44) and Fig. 1a suggest that different behaviours should be expected between ρ cst DN and ρ aff DN , or ρ cst DN and ρ par DN , for example:

[µ 1 > 1 and µ 2 > 1] =⇒ ρ aff/par DN ≈ ρ cst DN ( 45
)
The numerical conjecture expressed in (45) suggest that, as in previously covered simplified cases (e.g. Sec. 4.1.2), the particular physical setting of the problem (translated into the signs of µ j -1), over which the user has no control, bears an impact on the convergence properties of DN-derived SWR algorithms. For example, with parabolic ν j , convergence is less likely to occur as h 2 increases and h 1 decreases, since such a configuration would favour small µ 1 and large µ 2 . Hence, in addition to the physics of the problem (translating into the values of ν j and ∂ z ν j ), its geometric setting (h j ) also impacts the potential convergence.

Upper and lower bounds for ρ cst DN We now describe the derivation of an upper bound for ρ cst DN to further exploit inequalities (44). First, injecting (28) into ( 21) and reformulating yields:

ρ cst DN = λ Q( √ 2Fo 2 ) Q( √ 2Fo 1 ) , with Q(x) = cosh x -cos x cosh x + cos x . (46)
Introducing x ( ,0) ≈ 2.36502 and x ( ,1) ≈ 5.49780 the two smallest positive roots of the transcendental equation tanh x+tan x = 0, and denoting β = Q(x ( ,0) )/Q(x ( ,1) ) ≈ 1.32165, it can then be proved (see App. 3) that:

       λ/β ≤ ρ cst DN ≤ λ h 2 h 1 if √ λ > h 2 h 1 , λ h 2 h 1 ≤ ρ cst DN ≤ √ λβ if √ λ ≤ h 2 h 1 . (47)
Note that these upper (resp. lower) bounds are reached maxima (resp. minima). Moreover, when studying convergence, ω is a common parameter in Fo 2 and Fo 1 and |λh 2 /h 1 | = λFo 2 /Fo 1 with Fo 2 /Fo 1 independent of ω.

Upper and lower bounds for ρ aff/par DN

Replacing ρ cst DN by its upper bound in inequalities (44), the following upper bound for the convergence rate in the DN case with affine viscosities can be obtained :

ρ aff DN ≤ U aff DN :=      λ h 2 h 1 (µ 1 ln (1 + 1/µ 1 )) -1 if Fo 2 > Fo 1 , λβ (µ 1 ln (1 + 1/µ 1 )) -1 if Fo 2 ≤ Fo 1 . (48)
Similarly, with parabolic viscosities:

ρ par DN ≤ U par DN :=          λ h 2 h 1 √ 1 + 4µ 1 4µ 1 arccoth √ 1 + 4µ 1 if Fo 2 > Fo 1 , λβ(1 + 4µ 1 ) 4µ 1 arccoth √ 1 + 4µ 1 if Fo 2 ≤ Fo 1 . (49)
So far, we have considered that ω may vary in R. However, for a discretized problem, ω ∈ I ω , defined in [START_REF] Keyes | Multiphysics simulations: Challenges and opportunities[END_REF]. Limiting ω to bounded values leads to reformulating (42) as:

R aff,cst 1 (ω max ) ≤ R aff,cst 1 ≤ R aff,cst 1 (-f ), (50a)
R par,cst 1 (ω max ) ≤ R par,cst 1 ≤ R par,cst 1 (-f ), (50b)
since Fo 1 (-f ) = 0 and Fo 2 (ω max ) = max ω∈Iω (|Fo 2 (ω)|), and since R ν,cst j decreases with Fo j (see Fig. 2). Therefore, |R ν,cst

2 (ω)| ≤ R ν,cst 2
(ω max ) can be inferred from (50). The following more restrictive upper bounds can then be obtained by injecting this into (43):

ρ aff DN ≤ U aff DN R aff,cst 2 (ω max ), ρ par DN ≤ U par DN R par,cst 1 (ω max ), (51)
where U aff DN (resp. U par DN ) is the upper bound of ρ aff DN (resp. U par DN ) defined in (48) (resp. (49)), and R aff,cst 2 (ω max ), R par,cst 2 (ω max ) -→ ωmax→∞ 1 -. To our knowledge, similar upper bounds had already been derived only in the case of constant coefficients and subdomains of infinite size. In that case, they are all equal to ρ ∞,c DN = √ λ.

Convergence analysis for RR transmission conditions

The exact transparent conditions (leading to ρ = 0 for all ω, hence convergence in two iterations) can easily be found in the Fourier space. For our model problem (10), considering the interface operators

B 1,1 = Id + λ Λ 2 , B 1,2 = Id + Λ 2 , B 2,1 = Λ 1 + ∂ z and B 2,2 = Λ 1 + λ∂ z , convergence in two iterations is obtained for Λ 1 (ϕ 2 ) = F -1 (-λS 1 ϕ 2 ) and Λ 2 (ϕ 1 ) = F -1 (-S 2 ∂ z ϕ 1 )
where F -1 denotes the inverse Fourier transform. In the stationary case (ω = 0), the Λ j operators can be found in the physical space (Λ 1 = -λS sta 1 , Λ 2 = -S sta 2 ∂ z ), but in the general case ω = 0, these conditions are nonlocal in time. Nevertheless, it is possible to improve the convergence speed by using RR interface conditions as given in [START_REF] Lemarié | Analysis of ocean-atmosphere coupling algorithms: Consistency and stability[END_REF], which amounts to assimilating Λ 1 to pId and Λ 2 to q∂ z , with (p, q) ∈ P so that the problem is well-posed. A classical way to find suitable values for (p, q) is to solve the optimization problem: Find (p * , q * ) ∈ P such that max ω∈Iω ρ RR (ω, p * , q * ) = min (p,q)∈P max ω∈Iω ρ RR (ω, p, q) . (52) (52) has been analytically solved for simple cases with constant viscosity and infinite domain [e.g. 10, 21]. For more complicated settings, the min-max problem can only be solved numerically. In our case, with spatially-variable viscosity coefficients and bounded domains, the optimization can only be done numerically (see Sec. [START_REF] Brassington | Progress and challenges in short-to medium-range coupled prediction[END_REF]

.3 for specific examples).

We now aim at finding values of (p, q) ∈ P which would at least guarantee the convergence of the algorithm, without necessarily being the optimal choice. Let us first recall that the convergence rate with RR interface conditions is:

ρ RR (p, q, ω) = (S 2 (ω) + q)(p + λS 1 (ω)) (1 + pS 2 (ω))(1 + λqS 1 (ω)) . (53)
Convergence (i.e. ρ RR < 1) then amounts to:

|pq + λS 1 S 2 + λqS 1 + pS 2 | 2 < |1 + pqλS 1 S 2 + λqS 1 + pS 2 | 2 . ( 54
)
Ensuring SWR convergence thus boils down to finding sufficient conditions on (p, q) ∈ P for satisfying (54). In order to simplify this endeavor, we will from now on assume that:

(p, q) ∈ P ∩ (R * ) 2 = R * -× R * + , (55)
which is a limitation guided by pragmatic motivations. By splitting (54) into real and imaginary parts and relying on (55), it is possible to show that (54) is equivalent to: viscosity case (using S cst j defined in ( 28)), we have only been able to numerically check it in the affine and parabolic cases (see Fig. 3). Hence, from now on, we assume the left condition of (57) to be true, so that its right condition is satisfied as well. By injecting this into (56), we obtain:

ρ 2 DN -1 (pq + 1) + 2p Re (λS 1 ) |S 2 | 2 -Re (S 2 ) ϑ1 +2q Re (S 2 ) |λS 1 | 2 -Re (λS 1 )
(p, q) ∈ R * -× R * + and pq = -1 ⇒ {∀ω ∈ R, ρ RR (p, q, ω) < 1} . ( 58
)
The left condition of (58) is not optimal in terms of convergence speed, partly because of the limitations induced by ( 55), but it does ensure the convergence of the SWR algorithm with RR interface conditions, and its simple formulation can easily be implemented. Moreover, (58) illustrates the benefits of investigating RR-based SWR algorithms. Indeed, for DN-based ones, ρ DN is directly derived from the problem's setting, hence convergence cannot be guaranteed; for RR-based ones, the pq + 1 term in (56) permits modulating ρ RR and thus potentially converging for problem settings under which DN-based algorithms diverge.

We now aim at extending the analytically certified region of convergence in the (p, q) ∈ R * -× R * + space beyond the particular line pq = -1. For doing so, we still assume that the left condition of (57) is satisfied. Thanks to the bounds on ρ DN determined in Sec. 4.3, in some cases, it is possible to predict the constant sign of ρ 2 DN -1, and then to choose (p, q) ∈ R * -× R * + such that pq + 1 has the opposite sign of ρ 2 DN -1. Yet, in other cases, the sign of ρ 2 DN -1 can change with respect to ω. However, this issue can be circumvented at theoretical (constant ν j ) or practical (affine or parabolic ν j ) levels.

Proposition 2 If the viscosity profiles are constant on each subdomain (ν = ν c j ), and if (p, q) ∈ R * -× R * + are such that:

-p ≤ 2 Γ 0 λ |h 1 | and q ≥ (Γ 0 ) 2 2 |h 2 |, (59)
where Γ 0 ≈ 1.45529 is independent of the problem setting, then the SWR algorithm (10) with the RR transmission conditions defined in [START_REF] Lemarié | Analysis of ocean-atmosphere coupling algorithms: Consistency and stability[END_REF] converges for all ω ∈ R.

Proposition 3 If the viscosity profiles are constant on each subdomain (ν = ν c j ), and if (p, q) ∈ R * -× R * + are such that:

-p ≥ (Γ 0 ) 2 2 λ S cst 1 (ω max ) and q ≤ 2 Γ 0 S cst 2 (ω max ) , (60)
then the SWR algorithm (10) with the RR transmission conditions defined in [START_REF] Lemarié | Analysis of ocean-atmosphere coupling algorithms: Consistency and stability[END_REF] converges for all ω ∈ I ω .

Conjecture 1 Proposition 2 holds in the case of affine or parabolic viscosities. This means for all three choices of viscosity profiles, the choice of RR parameters providing the exact transparent conditions in the stationary case with the constant viscosity profiles ν cst j = ν j (0) (i.e. p 0 = -λ/|h 1 | and q 0 = |h 2 |) leads to SWR convergence (ρ RR (ω, p 0 , q 0 ) < 1).

Conjecture 2 In the cases of affine or parabolic ν j , if (p, q) ∈ R * -× R * + are such that either of these conditions are met:

-p ≤ |S ν 1 (ω max )| |S cst 1 (ω max )| (Γ 0 ) 2 2 λ S cst,sta 1,f =0 and q ≥ |S ν 2 (ω max )| |S cst 2 (ω max )| 2 Γ 0 S cst,sta 2,f =0 , (61a)
-p ≥ S ν,sta

1,f =0 S cst,sta 1,f =0 2 Γ 0 λ S cst 1 (ω max ) and q ≤ S ν,sta 2,f =0 S cst,sta 2,f =0 (Γ 0 ) 2 2 S cst 2 (ω max ) , (61b)
then the SWR algorithm (10) with the RR transmission conditions defined in [START_REF] Lemarié | Analysis of ocean-atmosphere coupling algorithms: Consistency and stability[END_REF] converges for all ω ∈ I ω . 59) and ( 61) can be substituted with these less constraining sets of conditions:

Conjecture 3 Moreover, if ω → |Re [S 1 (ω)] | is increasing and ω → |Re [S 2 (ω)] | is decreasing, then (
-p ≤ 2 Γ 0 λ S sta,ν 1,f =0 and q ≥ (Γ 0 ) 2 2 S sta,ν 2,f =0 , (62a)
or

-p ≥ (Γ 0 ) 2 2 λ |S ν 1 (ω max )| and q ≤ 2 Γ 0 |S ν 2 (ω max )| . (62b)
In App. 4, Props. 2 and 3 are proved and Conjs. 1, 2 and 3 are discussed. While Prop.

2 and (62a) offer sufficient condition for theoretical SWR convergence (ω ∈ R), Prop. 3 and (62b) are limited to its discrete implementation (ω ∈ I ω ). Interestingly, for all three choices of viscosity profiles, the choice of RR parameters providing the exact transparent conditions in the stationary case (i.e. p 0 = -λ S sta,ν 1,f =0 and q 0 = S sta,ν 2,f =0 ) leads to SWR convergence (ρ RR (ω, p 0 , q 0 ) < 1).

Numerical experiments

In this section, we illustrate the different properties highlighted in the previous sections with numerical experiments. In particular, we numerically investigate the influence of continuously variable coefficients on the convergence in the case of DN and RR interface conditions in view of our theoretical results. We directly simulate the equations for the errors, i.e. we set

F 1 = F 2 = 0, U 1,g = U 2,g = 0 in (3). z 0 z -1 z 1 z 2 z -2 Ω 1 Ω 2

Discretization

Because viscosity coefficients ν j (z) can have large variations in space, we use a nonuniform grid with increased resolution in the vicinity of the interface. On subdomain Ω j of length |h j |, the location of the cell center z k for grid cell k is

z k = h c σ k + (|h j | -h c ) sinh(σ k θ s ) sinh θ s , σ k = k N , k = 0, . . . , N ;
where N + 1 is the number of grid points for the discretization on Ω j , and (h c , θ s ) two parameters controlling the grid stretching. For our experiments, we choose h c = 10 -3 |h j |, θ s = 10 -1 N and N 10, which leads to a resolution of ∆z ≈ 10 -3 |h j | near the interface and ∆z ≈ 10 -1 |h j | near the external boundaries. As shown in Fig. 4, the interface is located in the middle of a grid cell at z 0 = 0. Regarding the discretization in time, in order to circumvent the stability issues of the Euler forward scheme to integrate the Coriolis term in (3), a Forward-Backward approach is used. For a velocity vector U = (u, v), this scheme reads:

u k = u n k + ∆t f v n k , v k = v n k -∆t f u k .
The order of integration of u and v is inverted from one time-step to the other to minimize splitting errors. The diffusion is then discretized using a standard implicit finite-difference scheme. This leads to a tridiagonal system for interior grid points, which reads (for u, with the same applying for v):

a k u n+1 k-1 + b k u n+1 k + c k u n+1 k+1 = u k , k = 1, . . . , N -1, (63)
where

a k = - ν k-1/2 ∆t ∆z k ∆z k-1/2 , c k = - ν k+1/2 ∆t ∆z k ∆z k+1/2
, and b k = 1a kc k , with

∆z k = z k+1/2 -z k-1/2 and ∆z k+1/2 = z k+1 -z k . Same applies to the v velocity component.
Since the interface is located at a cell center (Fig. 4), the evaluation of the gradient (∂ z u) 0 at the interface (required for Neumann and Robin conditions) is obtained via extrapolation with a second-order formula. Simulations have been performed using a multiplicative Schwarz algorithm initialized with a random guess at the interface, so that the observed error widely samples the temporal frequencies that can be represented by the computational grid.

Numerical experiments for Dirichlet-Neumann interface conditions

We compare the theoretical convergence rate ρ ϕ (ω) on I ω , with the observed one ρ k U (z = 0) defined in [START_REF] Gander | Dirichlet-Neumann waveform relaxation methods for parabolic and hyperbolic problems in multiple subdomains[END_REF], for the three viscosity profiles, and for fixed parameters h j , ν c j and ∂ z ν j (0). This comparison is illustrated by Fig. 5 and6. The left panels represent the theoretical convergence rates ρ cst DN (ω), ρ aff DN (ω) and ρ par DN (ω) for different parameter values, and the convergence rates obtained from numerical simulations. Since the latter vary from one Schwarz iteration to the next one, they are not restricted to a single line but are materialized by a shaded area delimited by the minimum and maximum observed values over I ω . Since we consider f = 0, the theoretical convergence rates are not symmetrical with respect to ω = 0, hence the figures feature two curves for each rates: one for ρ ϕ (ω > 0) and one for ρ ϕ (ω < 0), with |ω| as the x-axis. Several remarks can be drawn. First, the numerically observed convergence rates satisfy the theoretical bound given by [START_REF] Grisogono | A generalized Ekman layer profile with gradually varying eddy diffusivities[END_REF]. Indeed, the maximum of the curves is systematically above the grey shaded areas, which are bounded by min

k ρ k U2 (z = 0) and max k ρ k U2 (z = 0) (
where k is the iteration number). As expected, the convergence (or divergence) of the SWR algorithm depends on both the physical configuration and the nature of the viscosity profile (constant, affine or parabolic). For instance, as shown in the right panels of Fig. 5 and6, there are situations where the algorithm converges only for the constant viscosity case and diverges in the affine and parabolic cases (while the exact opposite behaviour can be found for other parameter values). These examples illustrate that the spatial variations of the diffusion coefficient is determining for the SWR algorithm convergence, and that results from one particular choice cannot be generalized to other ones.

Figure 6 confirms the consistency between the numerically observed convergence rate, computed in the physical space, and the theoretical convergence rate, computed in the Fourier space. The observed convergence rates (right panel) are smaller during the first iterations. This is due to the rapid damping of error components corresponding to the frequencies with fast convergence. As iterations progress, the errors at slower-converging frequencies have more influence on the L 2 norm, and thus on the observed convergence rate. This results in observed convergence rates which increasing at each iteration. Consistently, when the variations of the theoretical convergence rate in the Fourier space are small, the observed convergence rate is more stable throughout the iterations.

Figure 6 also shows that the peak corresponding to ω + f = 0 in the theoretical convergence rates indeed impacts the convergence at a numerical level. Such peak appears in the theoretical derivation only when considering finite size domains . Without these hypotheses, the curve for the convergence rate would be flat, corresponding to ρ(ω > 0). This example illustrates that the observed convergence rate is larger than the one given by the flat portion of the theoretical convergence rate, because its maximum value is influenced by the asymptotic case ω + f → 0. Adding the rotation (f = 0) shifts this peak by -f in the ω space. In the case where |f | < π/δt, the bounds for ρ obs are the same as with the Coriolis effect. But if |f | > π/δt, or in the stationary case, this shift could change the bounds and thus the convergence of the algorithm.

h 1 = -50 m, ν c 1 = 0.8 m 2 s -1 , ∂ z ν 1 (0) = -0.006 m s -1 , h 2 = 100.0 m, ν c 2 = 0.09 m 2 s -1 , ∂ z ν 2 (0) = 0.4 m s -1 , f = 5 × 10 -5 s -1 , ∆t = 2000 s.

Numerical experiments for RR interface

Here we focus on the calculation of Robin coefficients for optimizing the SWR algorithm convergence, as explained in Sec. 4.4. While previous studies have already determined the optimized coefficients in various cases with infinite domain, constants viscosity profiles and without Coriolis effect [e.g. 10, 11, 21], the assumptions made during their determination may lead to such coefficients being irrelevant in the cases investigated below. In Fig. 7, the behaviour of the SWR algorithm for a parabolic profile is shown, with four different choices of Robin coefficients leading to four different convergence rates:

(a) ρ RR,par 0 , corresponding to optimized coefficients analytically determined in [START_REF] Lemarié | Toward an optimized global-in-time Schwarz algorithm for diffusion equations with discontinuous and spatially variable coefficients, part 1: The constant coefficients case[END_REF] for infinite domain, constant viscosity and no Coriolis term (the corresponding convergence rate will be referred to as ρ RR,par , same as (b) with solving (52) on ρ par RR . As in the DN case, the numerically observed convergence rates agree well with the theoretical bound given by [START_REF] Grisogono | A generalized Ekman layer profile with gradually varying eddy diffusivities[END_REF]. Figure 7a displays the theoretical convergence rates ρ RR,par s (s ∈ {0 . . . 3}). This figure clearly illustrates that significant gain in performance can be obtained through optimizing Robin coefficients as consistently as possible with the practical problem of interest. As confirmed by the experiments shown in Fig. 7b, optimizing the Robin coefficients using a convergence rate which accounts for the parabolic viscosity profile, finite size subdomains and the Coriolis effect provides a very efficient algorithm. On the contrary, other alternatives, where simplifications had been made, are at best suboptimal and may even lead to divergence. Once again, neglecting the variations of the viscosity in the convergence analysis can conduct to critically erroneous choices of Robin coefficients.

Besides the determination of optimized Robin coefficients, Sec. 4.4 dealt with different ways of choosing acceptable coefficients without numerically solving the costly (52) min-max problem. Figure 8 maps the maximum (w.r.t. ω) of the convergence rate as a function of p < 0 and q > 0, for the three types of considered viscosity profiles, with specific symbols localizing the four different couples of (p, q) values discussed in previous paragraph, and used in Fig. 7. (p, q) values corresponding to optimized parameter values in the nonrotating and stationary case (i.e., p 0 = -λS sta,par 1,f =0 and q 0 = -S sta,par 2,f =0 ), are also displayed. The left panel of Fig. 8 underlines that using Robin coefficients obtained from optimizing the constant or affine viscosity case does not ensure good convergence of the model with a parabolic viscosity profile. These maps also display the areas of sufficient conditions for RR-based SWR algorithm convergence, established in Sec. 4.4. As explained in Sec. 4.4, the (p 0 , q 0 ) couple from the nonrotating and stationary case is located close to the white zone, which ensures convergence for all three viscosity profiles without requiring solving a costly optimization problem.

Conclusion

In this paper, we have analyzed the convergence of global-in-time Schwarz waveform relaxation algorithms on the coupled Ekman layer problem with continuously variable viscosity coefficients. We have obtained new theoretical results for problems featuring several delicacies which are very often neglected in existing convergence analyses: viscosity profiles are allowed to vary in space and the Coriolis effect as- (for constant ν j ) and Conjecture 3 (for affine and parabolic ν j ); light grey zones to observed convergence without having those conditions satisfied; dark grey zones to divergence. The first four symbols correspond to the values of (p, q) used in Fig. 7. The fifth symbol corresponds to (p 0 , q 0 ) obtained in the stationary case without rotation (Sec. 4.4). The parameter values are the same as in Fig. 7. sociated with the Earth's rotation is accounted for. First, we have emphasized that including the Coriolis effect amounts to investigating a reaction-diffusion equation in the complex space. We then have considered two types of SWR interface conditions (Dirichlet-Neumann (DN) and Robin-Robin (RR)), and three types of viscosity profiles (constant, affine and parabolic). Fourier space convergence rates have been given for all three cases, with indications on how such results are linked to empirical convergence rates observed from numerically simulated problems. In our physically realistic setting, we have defined analytical bounds of SWR convergence rate for constant viscosity profiles for the DN interface conditions, while for the RR case we have conjectured sufficient conditions ensuring convergence. Similar results extended to the affine and parabolic viscosity cases have also been numerically assessed. More importantly, we have emphasized that spatial variations in viscosity can be determining on the convergence of SWR algorithm, and that neglecting them can mislead preliminary convergence analyses. This work is part of a community effort to mathematically and numerically evaluate the coupling methods presently used in realistic climate models. Several open perspectives are left. First, some of this work's conclusions rely on numerical conjectures that have not been proved, but only observed. Second, generalizing our study to cubic viscosity profiles would be interesting, as this would include the profile proposed by [START_REF] O'brien | A note on the vertical structure of the eddy exchange coefficient in the planetary boundary layer[END_REF] which is used as a baseline in broadly-used turbulence schemes. In the same way as Bessel (resp. hypergeometric) functions are used for describing SWR algorithm solutions for affine (resp. parabolic) viscosity profiles, Heun's function is a good candidate for studying the cubic case. However, preliminary numerical studies have suffered from numerical instabilities arising from evaluating such functions. It should also be borne in mind that other recently proposed waveform relaxation algorithms derived from Dirichlet-Neumann conditions [e.g., [START_REF] Gander | Dirichlet-neumann and neumann-neumann waveform relaxation algorithms for parabolic problems[END_REF][START_REF] Gander | Dirichlet-Neumann waveform relaxation methods for parabolic and hyperbolic problems in multiple subdomains[END_REF] have not been investigated here, and may lead to convergence properties able to compete with the RR-based SWR algorithms presented here.

On the longer term, a significant challenge would be to investigate SWR algorithms in the context of viscosity profiles resulting from a turbulent kinetic energy (TKE) closure. Such methods, which are increasingly used in climate models, rely on solving an additional prognostic equation and result in space and time varying viscosity profiles depending nonlinearly on model variables. Investigating idealized testcases are a necessary and instructive endeavor for accurately seizing the main numerical challenges arising from model coupling.

(z) = 2a √ b 2 -4ac z + b 2a and ξ = - 1 2 1 ± 1 + 4ı (f + ω) a , (64)
where the choice of the sign in (64) does not bear any consequence. In particular, one can show that using relevant quadratic viscosity functions, (e.g. satisfying ν(z) > 0 for all z ∈ Ω and (a, b, c) ∈ R * -× R 2 ), η ∈ R and satisfies -1 < η(z) < 1 for z ∈ Ω. The functions given by ( 32) are then numerically satisfactory solutions to the Legendre ODE [meaning that computing a reasonable amount of sum terms leads to accurate values, 8], and thus to (12a). In our numerical results, Legendre functions have been computed using the hypergeometric function 2 F 1 [START_REF] Olver | Hypergeometric function[END_REF] through:

P 0 ξ (η) = 2 F 1 (ξ + 1, -ξ, 1; η).

On the observed discretized convergence factor

Let us consider the numerical solution of the coupled problem (3) on a given time window T n using a finite difference discretization in time. The observed error can be written as (with U j = (u j , v j )):

E k,m j (z) = U k,m j (z) -U j (z, t m j ) 1 ≤ m ≤ M j , j ∈ {1, 2}, k ∈ N * , ( 65
)
where m is an index for the time step, t m j is the physical time corresponding to time step m, and M j is the total number of time steps for the numerical approximation of U j in T n . A convergence rate characterizing the behaviour of E k,m j = (E k,m u,j , E k,m v,j ) as a function of k is:

R k Uj = E k,m j (0 ∓ ) m 2 E k-1,m j (0 ∓ ) m 2 = Mj m=1 E k,m u,j (0 ∓ ) 2 + E k,m v,j (0 ∓ ) 2 1/2 Mj m=1 E k-1,m u,j (0 ∓ ) 2 + E k-1,m v,j (0 ∓ ) 2 1/2 , k ∈ N * , j ∈ {1, 2}. (66)
Also the discrete errors on the complex variables (ϕ, ϕ) are:

e k,m j (z) = e k,m ϕ,j (z) e k,m ϕ,j (z) = ϕ k,m j (z) -ϕ(z, t j m ) ϕ k,m j (z) -ϕ(z, t j m ) , (67)
where P is given in Sec. 2.2. Since P is a unitary matrix (i.e. P -1 = P T ), {δtj } . This is due to the Nyquist-Shannon sampling theorem [START_REF] Dubois | The optimized Schwarz method with a coarse grid correction[END_REF]. e k,m ϕ,j (0 ∓ ) corresponds to the Fourier transform of the continuous error êk ϕ on I ω . The convergence factor denoted ρ k,∞ ϕ,j can be linked to the errors in the Fourier space by the Parseval's theorem: The derivative of Q(x) cancels for x such that tanh(x ) + tan(x ) = 0. This transcendental equation has no root on [0, π/2] (and 0 ≤ Q(x) ≤ 1 on that interval) and then one root x ( ,k) per interval [π/2 + kπ, 3π/2 + kπ] (k ≥ 0). For even values of k, Q(x ( ,k) ) is a local maximum and Q(x) ≥ 1, while for odd values of k, Q(x ( ,k) ) is a local minimum and Q(x) ≤ 1. Moreover it can be shown that Q(x ( ,0) ) > Q(x ( ,2) ) > Q(x ( ,4) ) > . . . > 1 and Q(x ( ,1) ) < Q(x ( ,3) ) < Q(x ( ,5) ) < . . . < 1. This implies that the maximum of Q(x) is Q(x ( ,0) ) and thus that Q(x) ≤ Q(x ( ,0) ), x ≥ 0.

If we now introduce α = Fo 2 /Fo 1 and consider x = √ 2Fo 1 , our problem is to find the upper bound of Q(αx) Q(x) with α ≤ 1 (because we consider the case Fo 2 < Fo 1 here). It is straightforward to prove that Q(αx) ≤ Q(x ( ,0) ), for x ≥ 0, α ≤ 1 and Q(αx) behaves like Q(x) on the subintervals [ π 2α + kπ, 3π 2α + kπ]. We can then proceed sub-interval by sub-interval: (i) on [0, π/2] Q(x) increases more rapidly than Q(αx). We thus have Q(αx) Q(x) ≤ 1 as soon as α ≤ 1; (ii) on [π/2, 3π/2] Q(x) is larger than 1, meaning that 1 Q(x) ≤ 1. We thus have Q(αx) Q(x) ≤ Q(x ( ,0) ) whatever the value of α; (iii) on [3π/2, 5π/2] Q(x) has a local minimum Q(x ( ,1) ) which is smaller than 1. We thus have 1 Q(x) ≤ 1 Q(x ( ,1) ) and the worst case scenario is when this coincides with values of α such that the maximum of Q(αx) occurs for x ∈ [3π/2, 5π/2]. We thus obtain: 1) ) .

Q(αx) Q(x) ≤ Q(x ( ,0) ) Q(x ( ,
Because of the ordering of the values Q(x ( ,k) ) given earlier (i.e. Q(x ( ,k) ) gets closer and closer to 1 as k increases), it is not needed to go beyond x = 5π/2 since the upper bound could not be larger for x > 5π/2.

3.2 Proof for the upper bound (47) for Fo 2 ≥ Fo 1

For this proof, we first study the behaviour of Q(x) x 2 . The sign of the derivative of Q(x)

x 2

with respect to x is the same as cos 2 (x)cosh 2 (x) + x(cosh(x) sin(x) + sinh(x) cos(x)), which is smaller or equal to zero for x ≥ 0, Q(x)

x 2 is thus a decreasing function of x for x ≥ 0. In the case α ≥ 1 (i.e. Fo 2 ≥ Fo 1 ) we have αx ≤ x, hence Q(αx)

α 2 x 2 ≤ Q(x) x 2
which leads to Q(αx) Q(x) ≤ α 2 where x = √ 2Fo 1 and α = Fo 2 /Fo 1 . We can rewrite this last inequality in terms of Fo j to obtain

ρ cst DN ≤ λ Fo 2 Fo 1 . (73)
Note that this also proves that |S

4. 1 . 1

 11 Stationary case without rotation (ω = 0, f = 0)

Fig. 1 :

 1 Fig. 1: Behaviours of various S j in the stationary case, (a) without and (b) with rotation .

4. 1 . 3

 13 Nonstationary case on infinite domains without rotation and with constant viscosity

2 and R par,cst 2 can 2 R par,cst 1 Fig. 2 : 1 ( 1 (

 2221211 Fig. 2: R aff,cst 1

2 ) 1 (

 21 are obtained by taking the inverse of R aff,cst resp. R par,cst 1

2 Fig. 3 :S aff 1 - 1

 2311 Fig. 3: Re S aff 1 -1 (left) and Re (S par 1 )-1 (right) with respect to Fo and µ. Note that both axes and the colorbar are in logarithmic scale. The hatched area corresponds to values of (Fo, µ) for which the computations of Bessel functions in S aff 1 suffer from numerical instability.

Fig. 4 :

 4 Fig. 4: Grid arrangement for the discretization of model problem (3) on subdomains Ω 1 and Ω 2 . Vertical lines indicate interfaces between grid cells while ellipses indicate cell centers. The vertical dashed line represents the interface between Ω 1 and Ω 2 .

2 ||eFig. 5 :

 25 Fig.5: Left panel: convergence rates with respect to |ω| ≤ π/∆t for three different viscosity profiles: constant (solid lines), parabolic (dotted lines) and affine (dashed lines). The two curves correspond to the theoretical convergence rate ρ(ω > 0) and ρ(ω < 0). The grey or hatched areas correspond to the values reached by the observed convergence rate in L 2 norm. Right panel: L 2 norm of the error on (u, v)(z = 0) as a function of the iteration number. Here the parameter values are:h 1 = -50 m, ν c 1 = 0.8 m 2 s -1 , ∂ z ν 1 (0) = -0.006 m s -1 , h 2 = 100.0 m, ν c 2 = 0.09 m 2 s -1 , ∂ z ν 2 (0) = 0.4 m s -1 , f = 5 × 10 -5 s -1 , ∆t = 2000 s.

1 , 2 ,Fig. 6 :

 126 Fig. 6: Left panel: same as left panel of Fig. 5 but for different parameter values. Right panel: observed convergence rates in L 2 norm on (u, v)(z = 0) as a function of the iterations. Both panels use the same y-axis. Here the parameter values are:h 1 = -50 m, ν c 1 = 0.012 m 2 s -1 , ∂ z ν 1 (0) = -0.04 m s -1 , h 2 = 200 m, ν c 2 = 0.06 m 2 s -1 , ∂ z ν 2 (0) = 0.01 m s -1 , f = 5 × 10 -5 s -1 and ∆t = 2000 s. (d) ρ RR,par 3

Fig. 7 :

 7 Fig. 7: (a): Theoretical convergence rate with respect to |ω| ≤ π/∆t for a given parabolic viscosity profile. The four couples of curves correspond to different choices of Robin coefficients (see Sec. 5.3). (b): Evolution of the L 2 -norm of the error on (u, v)(z = 0) with respect to the iterations for numerical simulations in the same setting as reported in (a). (c)-(f): details of the four cases from (a), one by one. Grey zones correspond to the amplitude of the values reached by the observed convergence rate over the iterations and the temporal frequencies. Parameter values are: h 1 = -50 m, ν c 1 = 0.06 m 2 s -1 , ∂ z ν 1 (0) = -0.001 m s -1 , h 2 = 200 m, ν c 2 = 0.012 m 2 s -1 , ∂ z ν 2 (0) = 0.04 m s -1 , f = 5 × 10 -4 s -1 , ∆t = 1000 s.

Fig. 8 :

 8 Fig. 8: Contour lines: max ω log 10 ρ RR (ω) respect to (-p, q) ∈ R 2 + for different viscosity profiles: (a) parabolic; (b) affine; (c) constant. The bold line indicates max ω log 10 ρ RR (ω) = 0; thin full (resp. dashed) lines are drawn with a ±0.2 increment. White zones correspond to areas of convergence guaranteed by Props. 2-3(for constant ν j ) and Conjecture 3 (for affine and parabolic ν j ); light grey zones to observed convergence without having those conditions satisfied; dark grey zones to divergence. The first four symbols correspond to the values of (p, q) used in Fig.7. The fifth symbol corresponds to (p 0 , q 0 ) obtained in the stationary case without rotation (Sec. 4.4). The parameter values are the same as in Fig.7.

  http://www.theses.fr/s184608 39. Vallis GK (2006) Atmospheric and Oceanic Fluid Dynamics: Fundamentals and Large-scale Circulation. Cambridge University Press 1 Quadratic viscosity basis functions Eq. (12a) using the quadratic viscosity ν(z) = c+bz +az 2 for z ∈ Ω, can be reduced to the Legendre ODE 1η 2 d 2 y dη 2 -2η dy dη + ξ (ξ + 1) y = 0 with η

2 1/ 2

 22 k ∈ N * , j ∈ {1, 2}. for k ∈ N * , j ∈ {1, 2}. When M j → +∞ (i.e. simulating (10) on a infinite time window), the temporal grid can only generate modes which frequencies lie in I ω := -

.3cosh x -cos x cosh x + cos x . 3 . 1

 31 as M j → ∞ (thanks to §2.1), (71) can be rewritten as ρ Finally the same relationships being true for e k,m ϕ,j , we end up with: Proof for the upper bound (47)In this appendix we provide the proof for the upper bound (47) on ρ cst DN for constant viscosity and Dirichlet-Neumann interface conditions given in Sec. 4.3. We first reformulate ρ cst DN as ρ cst DN = √ λ Proof for the upper bound (47) for Fo 2 < Fo 1

cst 2 ( 2 = 1 = |h 1 | - 1 . 4 Lemma 1 3 Fig. 9 : 1 Re S aff 1 ;

 22111413911 Fig. 9: Left: log 10 Im S aff 1

Acknowledgements

The authors acknowledge the support of the French national research agency (ANR) through contract ANR-16-CE01-0007.

then the following combination of conditions on (p, q) ∈ R * -× R * + guarantees the convergence of the RR-based SWR algorithm, for all ω ∈ R:

Proof (of Lemma 1) By reorganising (56), ρ RR (p, q, ω) < 1 is equivalent to :

Since pq < 0 and per hypothesis (57), we have ϑ 3 < 0. Let us now exploit the conditions (74) to obtain (76). Using ρ 2 DN = |λS 1 ||S 2 |, and recalling that Re(S 2 ) < 0 and q > 0, then

In the same way, because q > 0 and Re(S 1 ) > 0, -p < 2λ| Re(S 1 )| is a sufficient condition for ϑ 5 > 0.

Proof (of Prop. 2) As previously mentioned, we rely on Lemma 1. In the case of constant viscosities, then, using (28) and extensive reformulation (not detailed, analoguous to (46)) yields:

(77) Analogously to what is done in App. 3, it can then be proved that:

where y is the smallest positive root of tanh(y)tan(y) = 0. Hence, (74a) is satisfied. Moreover, still in App. 3, it is proved that ∀ω ∈ R,

Hence, (74b) and (74c) can be obtained by defining:

Injecting ( 79) into (75), and defining Γ 0 := 1 + (γ cst ) 2 , lead to (59).

Proof (of Prop.

3) The proof is analogous to that of Prop. 2, hence extensive detail is not provided. (78) still holds in that case with the same γ cst . Moreover, since here ω ∈ I ω R, we know that there exists

Analogously to the proof of Lemma 1, and using (80), it can then be shown that (76), and thus SWR convergence, are satisfied as soon as:

| is a suitable choice in (80). In the same way, since

Reinjecting this into (81) and relying on Γ 0 := 1 + (γ cst ) 2 lead to (60), which, as previously mentioned, is a sufficient condition for SWR convergence with RR transmission conditions for ω ∈ I ω .

Discussion on Conjectures 1, 2 and 3

As soon as the viscosity profiles are no longer constant, the validity of ( 57) and (74) can only be numerically assessed. From Fig. 3 we can conjecture that Re(S aff 1 ) ≥ 1 (resp. Re(S par 1 ) ≥ 1), which suggests that (57) and (74b) are satisfied on S 1 . Since S 2 (µ, Fo) = -|h 2 |/(|h 1 |S 1 (µ, Fo)), this would imply that (57) and (74c) are also satisfied for S 2 . Figure 9 shows numerical values of log 10 (γ), with γ = |Im(S ν 1 )|/|Re(S ν 1 )| in the cases where ν is affine and parabolic. It supports the idea that this ratio does not significantly depart from its constant viscosity value, which has been analytically found before. Then using result from App. 3 and (50) we find a lower bound to | Re(S ν 1 (ω))|:

Analogously, an upper bound to | Re(S ν 2 (ω))| can be found. Both bounds justify the first part of Conj. 1.

If ω ∈ I ω , then the bounds on | Re(S ν 1 (ω))| can be adjusted: