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

HAL
open science

# Semi-discrete analysis of a simplified air-sea coupling problem with nonlinear coupling conditions

Simon Clement, Florian Lemarié, Eric Blayo

► To cite this version:

Simon Clement, Florian Lemarié, Eric Blayo. Semi-discrete analysis of a simplified air-sea coupling problem with nonlinear coupling conditions. 27th International Domain Decomposition Conference, Jul 2022, Prague, Czech Republic. hal-04019469

HAL Id: hal-04019469

https://hal.science/hal-04019469

Submitted on 8 Mar 2023

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

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Semi-discrete analysis of a simplified air-sea coupling problem with nonlinear coupling conditions

Simon Clement, Florian Lemarié, and Eric Blayo

## 1 Introduction

This paper addresses the mathematical properties of a simplified nonlinear coupling problem representing the air-sea exchanges: it is shown in [4] that there is room for improvement in the coupling methods of state-of-the-art Earth System Models. Key ingredients in the ocean-atmosphere coupling are the computation, exchange and diffusion of turbulent fluxes at the interface. Those ingredients are represented by a nonlinear coupling condition which depends on the space discretization. We focus here on showing the existence of strong unsteady solutions of the semi-discrete in space problem with a nonlinear interface condition. Moreover we study the convergence of Schwarz Waveform Relaxation (SWR) applied to this coupled problem. In [2], the existence of unsteady solutions of nonlinear turbulent models for oceanic surface mixing layers is proven with the help of the inverse function theorem. After introducing the coupled problem in §2, its well-posedness is discussed by applying the inverse function theorem in §3. A convergence analysis of SWR is then pursued in §4 and complemented by numerical experiments detailed in §5.

<footnote>Simon Clement
Univ. Grenoble Alpes, CNRS, Inria, Grenoble INP, LJK, Grenoble, France, e-mail: simon.clement2@univ-grenoble-alpes.fr</footnote>

<footnote>Florian Lemarié
Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, LJK, Grenoble, France, e-mail: florian.lemarie@inria.fr</footnote>

<footnote>Eric Blayo
Univ. Grenoble Alpes, CNRS, Inria, Grenoble INP, LJK, Grenoble, France, e-mail: eric.blayo@univ-grenoble-alpes.fr</footnote>
## 2 Simplified air-sea coupled problem

We examine the solutions $U_a, U_o$ of coupled 1D linear reaction-diffusion equations which is a proxy for coupled ocean-atmosphere problems [3, 5]:

$$ (\partial_t + i f) U_j = v_j \partial_{zz} U_j + i f u_G^j, \quad (j = a, o), \quad (1) $$

where initial, boundary and coupling conditions will be specified together with the discretization. The reaction-diffusion equation (1) has a constant viscosity $v_j$ representing the turbulent vertical mixing and a complex reaction term if accounting for the Coriolis effect: $U_j$ are complex and include the horizontal direction of winds and currents (e.g. [5]). The Coriolis parameter $f$ is a constant determined by the latitude of the 1D vertical column considered. There is also a nudging term $u_G^j$ in (1) which pulls the solution towards the geostrophic equilibrium, i.e. a balance between the Coriolis force and the pressure gradient. In this study $u_G^a, u_G^o$ are constant.

The coupling between ocean and atmosphere actually excludes from the computational domains a surface layer $\Omega_{sl} = [z_{-\frac{1}{2}}, z_{\frac{1}{2}}]$ which is located between the first grid points of the two domains and contains the interface $z_0$. The size of this surface layer being linked with the discretization in space, we investigate the properties of the semi-discrete in space coupled problem. The solution $U_j$ ($j = o, a$) at the grid point $z_{m+\frac{1}{2}}$ ($m$ is the space index) is denoted $U_{j,m+\frac{1}{2}}$ and we denote $\phi_j$ the derivatives in space which are approximated with finite differences:

$$ \phi_{j,m} = \frac{1}{h_j} (U_{j,m+\frac{1}{2}} - U_{j,m-\frac{1}{2}}), \quad m \neq 0, \quad (2) $$

where $h_j$ is the space step assumed constant in each subdomain. The function $\phi_{j,0}$ will be determined by the boundary conditions involved in the coupling. The semi-discrete in space coupled problem is the following (for $\frac{H_o}{h_o} < m + \frac{1}{2} < \frac{H_a}{h_a}$, where $H_j$ are the size of the spatial domains):

$$ (\partial_t + i f) U_{j,m+\frac{1}{2}} = v_j \frac{\phi_{j,m+1} - \phi_{j,m}}{h_j} + i f u_G^j, \quad t \in ]0, T] \quad (3a) $$

$$ U_{j,m+\frac{1}{2}}|_{t=0} = U_0, \quad U_j|_{z=H_j} = U_j^\infty, \quad t \in ]0, T] \quad (3b) $$

$$ v_a \phi_{a,0} = C_D |U_{a,\frac{1}{2}} - U_{o,-\frac{1}{2}}|(U_{a,\frac{1}{2}} - U_{o,-\frac{1}{2}}), \quad t \in ]0, T] \quad (3c) $$

$$ \rho_o v_o \phi_{o,0} = \rho_a v_a \phi_{a,0}, \quad t \in ]0, T] \quad (3d) $$

The initial condition $U_0$ is chosen as the steady state (derived in Section 3.1) and the geostrophic winds and currents are prescribed as boundary conditions: $U_j^\infty = u_G^j$. The spatial extent of the domains will be considered sufficiently large ($H_j \to \infty$) but all the results can be easily extended to finite domains. The coupling conditions are composed of a quadratic friction law (3c) and of a flux continuity (3d). Those conditions are representative of the ones that can be found in more realistic models:
they are simpler, because $C_D$ and $\nu_j$ are assumed to be constant instead of depending themselves on $U_{a, \frac{1}{2}} - U_{o, -\frac{1}{2}}$. The densities $\rho_j$ are such that $\frac{\rho_a}{\rho_o} \approx 10^{-3}$.

### 3 Well-posedness

In this section the focus is on the well-posedness of (3). First the steady state is given, and the existence and uniqueness of solutions are proven in a neighborhood of this steady state. We assume in this section that $f \neq 0$: it can be proven that there is otherwise no bounded steady state of (3).

#### 3.1 Steady state

The derivation of the steady state of (3a) is somewhat similar to [6] and we obtain $\phi_{o,-m}^e = A^e(\lambda_o^e + 1)^m$ and $\phi_{a,m}^e = B^e(\lambda_a^e + 1)^m$ with $\lambda_j^e = \frac{1}{2}(\chi_j - \sqrt{\chi_j}\sqrt{\chi_j + 4})$, $\chi_j = \frac{ifh_j^2}{\nu_j}$. The continuity of the flux (3d) gives $\frac{A^e}{B^e}$ and finding a steady state solution amounts to find $\tilde{x} \in \mathbb{C}$ such that

$$ \tilde{x} - (u_G^a - u_G^o) = d|\tilde{x}|\tilde{x}, \quad (4) $$

where $d = C_D \left( \frac{\lambda_a^e}{ifh_a} + \frac{\rho_a}{\rho_o} \frac{\lambda_a^e}{ifh_o} \right)$ and $\tilde{x} = A^e \frac{\rho_o \nu_o}{\rho_a C_D} d + (u_G^a - u_G^o)$. The steady state corresponds to a solution $\tilde{x}$ whose modulus $|\tilde{x}|$ is a real and non-negative root of a polynomial:

$$ |\tilde{x}| = \frac{d_R}{2|d|^2} \pm \frac{\sqrt{\zeta + \gamma}}{2} \pm \frac{1}{2} \sqrt{2\zeta - \gamma \pm \frac{2d_R^3 - 2d_R|d|^2}{|d|^6\sqrt{\zeta + \gamma}}}, \quad d_R = \Re(d), \quad (5a) $$

$$ \gamma = \frac{\left(\frac{s}{2} + \frac{1}{2}\sqrt{s^2 - 4\beta^3}\right)^{\frac{1}{3}}}{3|d|^2} + \frac{\beta}{3|d|^2 \left(\frac{s}{2} + \frac{1}{2}\sqrt{s^2 - 4\beta^3}\right)^{\frac{1}{3}}}, \quad \zeta = -\frac{2}{3|d|^2} + \frac{d_R^2}{|d|^4}, \quad (5b) $$

$$ s = 2 + (u_G^a - u_G^o)^2 (72|d|^2 - 108d_R^2), \quad \beta = 1 - 12|d|^2 (u_G^a - u_G^o)^2, \quad (5c) $$

where the first $\pm$ and the third one in (5a) are necessarily the same. For our parameters there is only one combination of $\pm$ in (5a) for which $|\tilde{x}|$ is real and non-negative. There is hence only one steady solution of (3).

Finally we recover $\tilde{x}$ from (4): $\tilde{x} = \frac{u_G^a - u_G^o}{1 - d|\tilde{x}|}$, then $A^e = \frac{C_D \rho_a}{\rho_o \nu_o d} (\tilde{x} - (u_G^a - u_G^o))$ and $B^e = \frac{\rho_o \nu_o}{\rho_a \nu_a} A^e$. The steady state $U_j^e$ is given by
$$ 
\begin{align} 
U_{o,-m-1/2}^e &= u_G^o - \frac{v_o \lambda_o^e}{i f h_o} (1 + \lambda_o^e)^m A^e, \\ 
U_{a,m+1/2}^e &= u_G^a + \frac{v_a \lambda_a^e}{i f h_a} (1 + \lambda_a^e)^m B^e. 
\end{align} 
\tag{6} $$

Fig. 1 shows that this analysis exactly fits the numerical solution.

<figure><img src="image_2.png" /><figcaption>Fig. 1 Stationary solution profile in the ocean (bottom) and the atmosphere (top); A numerical steady state computed with 10 Schwarz iterations is the continuous red lines and the theoretical steady state obtained is displayed with dashed blue lines. Notice that the surface layer is not explicitly computed.</figcaption></figure>

## 3.2 Existence of solutions of the nonlinear semi-discrete in space problem

The method used by [2] to prove the existence and unicity of a solution in the neighborhood of a steady state can be used to deal with several types of nonlinearities. In particular, we can prove the existence and uniqueness of a solution to the problem (3) close to the steady state (i.e. with initial condition, boundary conditions and nudging terms close to $U_0, U_j^\infty, u_G^j$), thanks to the following steps:

1. The existence of a steady state $U^e$ is discussed in §3.1.

2. The well-posedness of the coupled problem with the linearized transmission conditions is proven in Appendix 6.A of [3].

3. The use of the inverse function theorem can be done in four steps:
a. concatenate the state vectors in a single vector $\mathbf{U} = \{U_a, U_o, \phi_a|_{z=0}\} \in \mathcal{U}$ where $\mathcal{U} = (L^2([0,T]))^{M_a+M_o+1}$ and $M_a, M_o$ are the number of grid levels in the subdomains. The functions $\phi_a|_{z=0}$ and $\phi_o$ are not in $\mathbf{U}$ because they can be expressed as linear combinations of elements of $\mathbf{U}$.

b. Define a mapping $\mathbf{\Phi}: \mathcal{U} \to \mathbf{Y}$ such that

$$ \mathbf{\Phi}(\mathbf{U}) = \{(\partial_t + \mathrm{i}f)U_a - \nu_a \partial_z \phi_a - g_a, (\partial_t + \mathrm{i}f)U_o - \nu_o \partial_z \phi_o - g_o, \\ U_a(H_a, t) - U_a^\infty, U_o(H_o, t) - U_o^\infty, U_a|_{t=0} - U_a^e, U_o|_{t=0} - U_o^e, \\ \nu_a \phi_a|_{z=0} - C_D \left| U_{a, \frac{1}{2}} - U_{o, -\frac{1}{2}} \right| \left( U_{a, \frac{1}{2}} - U_{o, -\frac{1}{2}} \right)\}, \quad (7) $$

where $\partial_z \phi_j$ is to be understood in the finite difference sense. Let us draw some important remarks about $\mathbf{\Phi}$:

* The equation $\phi_{o,0} = \frac{\rho_a \nu_a}{\rho_o \nu_o} \phi_{a,0}$ at interface is implicit: in (7), $\partial_z \phi_o$ at the first grid level is $\frac{\rho_a \nu_a}{\rho_o \nu_o} \phi_{a,0} - \frac{1}{h_o} (U_{o, -\frac{1}{2}} - U_{o, -\frac{3}{2}})$.

* The function $\mathbf{\Phi}$ is such that $\mathbf{\Phi}(\mathbf{U}^e) = 0$ where $\mathbf{U}^e$ is the steady state;

* The codomain $\mathbf{Y}$ is

$$ \mathbf{Y} = (L^2([0, T]))^{M_a-1} \times (L^2([0, T]))^{M_o-1} \times (L^2([0, T]))^3 \times \mathbb{R}^{M_a+M_o}. \quad (8) $$

* Finding $\mathbf{\Phi}^{-1}(y)$ is equivalent to solving the nonlinear semi-discrete problem (3) if the component of $y$ corresponding to the interface condition is zero (the other components correspond to other forcing terms, boundary conditions and initial condition). The idea of the proof is that if $\mathbf{\Phi}$ is invertible around $\mathbf{U}^e$ then the nonlinear semi-discrete problem (3) is invertible. Moreover, the inverse function theorem also tells us that $\mathbf{\Phi}^{-1}$ is continuous: this means that around the equilibrium state, the problem (3) is well-posed: it has a unique solution that depends continuously on the initial data.

c. Prove that $\mathbf{\Phi}$ is $C^1$ in a neighborhood of $\mathbf{U}^e$. The function $\mathbf{\Phi}$ is linear except for the transmission condition. Besides, the nonlinearity in this transmission condition is the function $x \mapsto |x|x$, which is continuously differentiable in a ball that does not contain zero. It is then straightforward to show that $\mathbf{\Phi}$ is $C^1$ and that its differential $D\mathbf{\Phi}(\mathbf{U}^e)$ is given by the linearized problem (a rigorous proof that can be directly adapted here is given in [2]).

d. Prove that $D\mathbf{\Phi}(\mathbf{U}^e)$ is an isomorphism: this is where the well-posedness of the coupled problem with linearized transmission conditions intervenes. The differential $D\mathbf{\Phi}(\mathbf{U}^e)$ corresponds indeed to the linearized problem with additional input data.

The next section uses the idea of considering the linearized problem around the steady state in the context of examining the convergence of a SWR algorithm.
4 Convergence analysis

In this section we conduct a convergence analysis of the SWR algorithm applied to
the coupled problem (initial and boundary conditions are similar to (3) and omitted):

$$
(\partial_t + \mathrm{if}) U_{j,m+\frac{1}{2}}^k = v_j \frac{\phi_{j,m+1}^k - \phi_{j,m}^k}{h_j} + \mathrm{if} u_G^j, \quad (9a)
$$

$$
\nu_a \phi_{a,0}^k = C_D \left| U_{a, \frac{1}{2}}^{k-1} - U_{o, -\frac{1}{2}}^{k-1} \right| \left( U_{a, \frac{1}{2}}^{k-1+\theta} - U_{o, -\frac{1}{2}}^{k-1} \right), \quad (9b)
$$

$$
\rho_o \nu_o \phi_{o,0}^k = \rho_a \nu_a \phi_{a,0}^k, \quad (9c)
$$

where $k$ is an iteration index and $U_{a,1/2}^{k-1+\theta} = (1-\theta)U_{a,1/2}^{k-1} + \theta U_{a,1/2}^k$. $\theta > 0$ is a
relaxation parameter: the goal of the convergence analysis pursued in this section is
to find an adequate value of $\theta$ leading to a fast convergence. The time windows is also
assumed to be of infinite size: this hypothesis is necessary to carry a Fourier analysis
of the convergence; finite time windows were handled at the semi-discrete in time
level [1] but the method is not straightforward to extend to discretized equations in
space.

The analysis of the SWR algorithm with the nonlinear transmission condition
cannot be directly pursued through a Fourier transform. We thus consider the lin-
earization of the problem around a steady state $U_j^e$, $\phi_j^e$ defined in (6). Assuming
that $U_j^k$ is in a neighborhood of $U_j^e$, the modulus in (9b) is non-zero and is thus
differentiable (close to zero, one could also smooth the modulus). Differences with
the steady state at the interface are noted $\delta\phi_j^k = \phi_{j,0}^k - \phi_{j,0}^e$, $\delta U_a^k = U_{a,\frac{1}{2}}^k - U_{a,\frac{1}{2}}^e$ and
$\delta U_o^k = U_{o,-\frac{1}{2}}^k - U_{o,-\frac{1}{2}}^e$ (note that we omit the $\frac{1}{2}$ in $\delta U_j^k$).

The linearized transmission operator involves the complex conjugate $\overline{\delta U_j^{k-1}}$:

$$
\nu_a \delta \phi_a^k = \alpha^e \left( \left( \frac{3}{2} - \theta \right) \delta U_a^{k-1} + \theta \delta U_a^k - \frac{3}{2} \delta U_o^{k-1} + \frac{\mathcal{R}^e}{2} \delta \overline{U_a^{k-1} - \delta U_o^{k-1}} \right) \quad (10)
$$

with $\alpha^e = C_D |U_{a,1/2}^e - U_{o,-1/2}^e|$ and $\mathcal{R}^e = \frac{U_{a,1/2}^e - U_{o,-1/2}^e}{U_{a,1/2}^e - U_{o,-1/2}^e}$. The relation (10) is used in
the convergence analysis instead of (9b).

We now follow [6] to derive a convergence factor of the SWR method applied
to the linearized transmission condition. It yields notably that the Fourier transform
of $U_{a,m+1/2}^k$ is $\tilde{U}_{a,m+1/2}^k = B_k(\lambda_a + 1)^m$ with $\lambda_j = \frac{1}{2}(\tilde{\chi}_j - \sqrt{\tilde{\chi}_j}\sqrt{\tilde{\chi}_j + 4})$ and $\tilde{\chi}_j = \frac{(f+\omega)ih_j^2}{\nu_j}$. We find that the evolution of $B_k$ is:

$$
\begin{align*}
B_{k+1}(\omega) &= a_1(\omega)B_k(\omega) + a_2(\omega)\overline{B_k(-\omega)}, \quad \text{where} \quad a_1 = \alpha^e \frac{\frac{3}{2}\mu(\omega) - \theta}{\frac{\nu_a}{h_a}(\lambda_a - \tilde{\chi}_a) - \alpha^e\theta}, \\
a_2 &= \alpha^e \frac{\frac{\mathcal{R}^e}{2}\overline{\mu(-\omega)}}{\frac{\nu_a}{h_a}(\lambda_a - \tilde{\chi}_a) - \alpha^e\theta} \quad \text{and} \quad \mu(\omega) = 1 - \frac{\lambda_a - \tilde{\chi}_a}{\lambda_o}\epsilon\frac{\nu_a h_o}{\nu_o h_a}.
\end{align*}
$$
Note that the variable $\bar{i\omega} = -i\omega$ appears when using the Fourier transform on $\delta U_a^{k-1} - \delta U_o^{k-1}$. As a consequence, the convergence factor $\xi^q$ in the linearized quadratic friction case differs from one iteration to another: it is a function of $\frac{B_{k-1}(-\omega)}{B_{k-1}(\omega)}$. We need to examine the evolution of both $B_{k+1}(\omega)$, $\overline{B_{k+1}(-\omega)}$:

$$ \left( \frac{B(\omega)}{B(-\omega)} \right)_{k+1} = \mathbf{M} \left( \frac{B(\omega)}{B(-\omega)} \right)_k, \quad \mathbf{M} = \left( \frac{a_1(\omega)}{a_2(-\omega)} \frac{a_2(\omega)}{a_1(-\omega)} \right) \quad (11) $$

The singular values (shown in Fig. 2) of $\mathbf{M}$ can be studied instead of the convergence factor. One can see on Fig. 2 that its two singular values $\tilde{\xi}_1$ and $\tilde{\xi}_2$ are different for small frequencies, especially around the frequencies $f$ and $-f$.

One can hence expect that for frequencies close to $f$ and $-f$ the “convergence factor” $\sqrt{|B_k(\omega)|^2+|B_{k-1}(\omega)|^2}$ will be different from one iteration to another and will be between $\tilde{\xi}_1$ and $\tilde{\xi}_2$.

We optimize only the maximum over the frequencies of the largest singular value $\tilde{\xi}_1$ (see Fig. 3) and find that the optimal value of $\theta$ is slightly smaller than 1.5. It is

<figure><img src="image_3.png" /><figcaption>**Fig. 3** Singular values $\tilde{\xi}_1$, $\tilde{\xi}_2$ of $\mathbf{M}$ maximized over the set of discrete frequencies $\{-\frac{\pi}{\Delta t}, \dots, \frac{\pi}{T}, 0, \frac{\pi}{T}, \dots, \frac{\pi}{\Delta t}\}$ as a function of $\theta$. The vertical dashed line highlights the minimum of $\max \tilde{\xi}_1$. The windows length $T$ is one day and $\Delta t = 60$ s.</figcaption></figure>

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

also seen in Fig. 3 that for $\theta < \frac{1}{2}$, $\tilde{\xi}_1 > 1$, which means that the SWR algorithm does not converge. On the contrary, for $\theta$ larger than $\frac{1}{2}$ both singular values are smaller than one for all the frequencies, and the convergence factor is then bounded by $\tilde{\xi}_1$.
# 5 Numerical experiments

Parameters of the numerical experiments are $C_D = 1.2 \times 10^{-3}$, $h_a = 20$ m, $h_o = 2$ m, $H_o = H_a = 2000$ m. The Coriolis parameter is $f = 10^{-4} \text{ s}^{-1}$ and the diffusivities are $v_a = 1 \text{ m}^2 \text{ s}^{-1}$, $v_o = 3 \times 10^{-3} \text{ m}^2 \text{ s}^{-1}$. The boundary conditions and nudging terms $U_j^\infty = u_G^j$ are set to constant values of $10 \text{ m s}^{-1}$ in the atmosphere and $0.1 \text{ m s}^{-1}$ in the ocean, while the initial condition is the steady state $U_0(z) = U_j^e(z)$. SWR is initialized at the interface with a white noise around the interface value of the steady state. Fig. 1 and Fig. 2 show that with those parameters, the theoretical results are coherent with the numerical experiments. Moreover Fig. 4 shows the evolution of the

<figure><img src="image_5.png" /><figcaption>**Fig. 4** Evolution of the $L^2$ norm of the errors for $\theta = \frac{1}{2}$ (green) and $\theta = \frac{3}{2}$ (purple). The singular value $\tilde{\xi}_1$ gives an upper bound of the error represented by dashed lines.</figcaption></figure>

error for two choices of $\theta$. As it is expected the choice $\theta = \frac{1}{2}$ leads to a convergence rate of approximately 1 whereas a relatively fast convergence is obtained with $\theta = \frac{3}{2}$.

# References

[1] Arnoult A, Japhet C, Omnes P (2022) Discrete-time analysis of Schwarz wave-form relaxation convergence. $\langle \text{hal-03746438} \rangle$

[2] Chacon-Rebollo T, Gomez-Marmol M, Rubino S (2014) On the existence and asymptotic stability of solutions for unsteady mixing-layer models. Discrete Contin Dyn Syst 34(2):421–436

[3] Clement S (2022) Numerical analysis for a combined space-time discretization of air-sea exchanges and their parameterizations. Phd thesis, UGA, $\langle \text{tel-03822632} \rangle$

[4] Marti O, Nguyen S, Braconnot P, Valcke S, Lemarié F, Blayo E (2021) A Schwarz iterative method to evaluate ocean–atmosphere coupling schemes: implementation and diagnostics in IPSL-CM6-SW-VLR. Geosci Model Dev 14(5):2959–2975

[5] Thery S, Pelletier C, Lemarié F, Blayo E (2021) Analysis of Schwarz waveform relaxation for the coupled Ekman boundary layer problem with continuously variable coefficients. Numer Algorithms 89:1145–1181

[6] Wu SL, Al-Khaleel MD (2014) Semi-discrete Schwarz waveform relaxation algorithms for reaction diffusion equations. BIT Numer Math 54(3):831–866