# Analysis and optimization of Schwarz algorithms for ocean-sea ice-atmosphere coupling

Pierre Lozano

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

## To cite this version:

Pierre Lozano. Analysis and optimization of Schwarz algorithms for ocean-sea ice-atmosphere coupling. Numerical Analysis [math.NA]. 2022. hal-05344227

HAL Id: hal-05344227

https://inria.hal.science/hal-05344227v1

Submitted on 3 Nov 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Master Research Project Report

## Analysis and optimization of Schwarz algorithms for ocean-sea ice-atmosphere coupling

Pierre Lozano

March 1, 2022 – July 29, 2022

**Laboratoire Jean Kuntzmann**
AIRSEA Team
700 avenue Centrale
38401 Saint Martin d'Hères

**Scientific advisors**
Eric Blayo (UGA)
Florian Lemarié (Inria)
Olivier Marti (LSCE)
# Contents

<table><tr><td>1</td><td>Introduction</td><td>2</td></tr><tr><td>2</td><td>Schwarz Algorithm</td><td>5</td></tr><tr><td>2.1</td><td>Coupling of Physical Models</td><td>5</td></tr><tr><td>2.2</td><td>Algorithms</td><td>6</td></tr><tr><td>2.3</td><td>Convergence Analysis</td><td>6</td></tr><tr><td>2.4</td><td>Coupled Stationary 1D Diffusion Equation</td><td>7</td></tr><tr><td>2.5</td><td>Coupled Nonstationary 1D Diffusion Equation</td><td>9</td></tr><tr><td>3</td><td>Atmosphere-Ice-Ocean Model</td><td>10</td></tr><tr><td>3.1</td><td>Motivations</td><td>11</td></tr><tr><td>3.2</td><td>The Model</td><td>11</td></tr><tr><td>3.2.1</td><td>A Balance Between Heat Fluxes</td><td>14</td></tr><tr><td>3.2.2</td><td>The System of Equations</td><td>15</td></tr><tr><td>3.3</td><td>Solving the Ice Thickness Evolution Equation</td><td>16</td></tr><tr><td>4</td><td>Applying the Schwarz Iterative Method</td><td>17</td></tr><tr><td>4.1</td><td>Further Assumptions</td><td>18</td></tr><tr><td>4.2</td><td>Iterative Algorithm</td><td>19</td></tr><tr><td>4.3</td><td>Errors and Convergence Analysis</td><td>19</td></tr><tr><td>5</td><td>Numerical Simulations</td><td>20</td></tr><tr><td>5.1</td><td>Physical Settings</td><td>20</td></tr><tr><td>5.2</td><td>Results</td><td>22</td></tr><tr><td>5.3</td><td>More General Case</td><td>22</td></tr><tr><td>6</td><td>Concerning Oscillations</td><td>23</td></tr><tr><td>7</td><td>Conclusion and Future Work</td><td>25</td></tr><tr><td>8</td><td>Schwarz Optimization</td><td>26</td></tr><tr><td>A</td><td>Sturm-Liouville Theory</td><td>28</td></tr><tr><td>A.1</td><td>Sturm-Liouville Operator</td><td>28</td></tr><tr><td>A.2</td><td>Main Results</td><td>29</td></tr><tr><td>A.3</td><td>The Case of the 1D Nonstationary Diffusion Equation</td><td>30</td></tr><tr><td>B</td><td>Chini's Equation</td><td>32</td></tr><tr><td>B.1</td><td>Existence of a Solution</td><td>33</td></tr><tr><td>B.2</td><td>Lambert W-function</td><td>33</td></tr><tr><td>B.2.1</td><td>Simple Case</td><td>34</td></tr><tr><td>B.2.2</td><td>Autonomous Case</td><td>34</td></tr></table>
# 1 Introduction

State-of-the-art Earth system models consist in the interaction of different components such as ocean, atmosphere, soils or ice caps. These components interact with each other at their boundaries, thus creating a strong coupling between them. An overview of the components taken into account in a climate model is presented in figure 2. In this thesis, we will focus on the interactions between atmosphere, ocean and sea ice.

<figure><img src="image_3.png" /><figcaption>Figure 2: A climate model schematic.</figcaption></figure>

To numerically achieve this coupling, usual algorithms generally lack synchronization. These temporal inconsistencies come from the fact that when one component evolves, either it is isolated from the others or the others are fixed. There are two main algorithms: a parallel one and a sequential one. Consider as an example the coupling between the ocean and the atmosphere. Figure 3 shows how both algorithms work during one coupling time window. A coupling time window is different from a time step: it corresponds to several time steps for each models.

In the parallel version (see 3a), ocean and atmosphere run in parallel. As interface flux for the coupling window $[t, t + \Delta t]$, the ocean sees the average flux on the previous coupling window $[t - \Delta t, t]$, which is computed by the atmosphere using the ocean values on $[t - 2\Delta t, t - \Delta t]$. It works similarly for the boundary condition of the atmosphere. Concerning the sequential atmosphere-first version (see 3b), ocean and atmosphere run one after the other starting with the atmosphere. The interface flux used as boundary condition by the atmosphere on $[t, t + \Delta t]$ relies on the time average ocean values on $[t - \Delta t, t]$. Then, the ocean sees as boundary condition for the coupling window $[t, t + \Delta t]$ the time average flux on $[t, t + \Delta t]$ computed using the atmospheric values.

The question that arises from these mathematical inconsistencies is the possible impact on the obtained result. To address this synchronisation issue, iterative algorithms such as the Schwarz iterative method can be used. The idea is to keep running the ocean and the atmosphere on the same coupling window and update the interface fluxes
<figure><img src="image_4.png" /><figcaption>**Figure 3:** The two usual algorithms used to achieve coupling of models. **A** and **O** are the state vectors of the atmosphere and the ocean respectively. They encompass several physical quantities such as temperature, velocity, humidity or salinity. **f** corresponds to the exchanged fluxes at the ocean-atmosphere interface Ω. **F**_A and **F**_O are partial derivative operator modeling the physics of each component.</figcaption></figure>

used as boundary condition until convergence. As before, we can choose to run both
components in parallel or in sequential. An illustration is presented in figure 4 where
ocean and atmosphere are run in parallel.

<figure><img src="image_5.png" /><figcaption>**Figure 4:** Schwarz algorithm for the ocean-atmosphere coupling. The subscript *k* corresponds to the *k*-th Schwarz iteration while *\* denotes the converged solution.</figcaption></figure>

The usual algorithms actually correspond to a single iteration of the Schwarz method.
More details are given on those methods in section 2.

It has recently<fnref n="1" /> been shown using Schwarz algorithms in large scale climate models,

<footnote n="1">¹During the ANR project COCOA (comprehension coupling approach for the ocean and the atmosphere) ended in September 2021.</footnote>
such as CNRM-CM<fnref n="2" /> or IPSL-CM<fnref n="3" />, that this dephasing has a harmful impact on the climate models [1], notably on the sea surface temperature as shown on figure 5a. With iterations obviously comes the question of the computational cost, especially for these types of simulations which are very demanding in terms of computer resources. What has been observed is that for the ocean-atmosphere coupling, only few iterations are sufficient to significantly reduce the error made with a single iteration. However, more iterations are needed for the algorithm to converge in regions with sea ice as we can see on figure 5b. We also notice that with a longer coupling period Δt, the number of iterations is even greater. Besides, oscillations of the ice temperature are often observed as the iterations proceed.

<figure><img src="image_6.png" /><figcaption>(a) Impact of the lagged algorithms on the sea surface temperature. The monthly average for an ocean basin over a set of 50 simulations can differ by more than 1°C</figcaption></figure>

<caption>(b) Number of Schwarz iterations needed to reach convergence</caption>

**Figure 5:** Impact of the iterative algorithm on the sea surface temperature and the number of iterations.

The goal of the project is thus to understand and explain in a scientific manner these two phenomena:

* slow convergence in regions with sea ice

* oscillations of the ice temperature during convergence

To do so, we propose to analyse the Schwarz algorithm applied on a simple one dimensional vertical model of the ocean-ice-atmosphere coupling problem. The three media model is presented in section 3. Analysis of the iterative algorithm is derived in section 4. Numerical simulations will then be run on this model in order to validate the theoretically computed rate of convergence. The results are exposed in section 5. Finally, early conclusion on the achieved work and perspective are provided in section 7.

<footnote n="2">²CNRM-CM: Centre National de Recherches Météorologiques - Climate Model</footnote>

<footnote n="3">³IPSL-CM: Institut Pierre Simon Laplace - Climate Model</footnote>
## 2 Schwarz Algorithm

Before going into more details on the thermodynamics of the three components at stake, we present in this section an iterative algorithm used to solve domain decomposition problems.

### 2.1 Coupling of Physical Models

The Schwarz algorithms are iterative methods to solve a problem by splitting it in subproblems with specific boundary conditions. Consider for example a domain $\Omega \subset \mathbb{R}^n$. We denote its border by $\partial\Omega$. Assume we want to solve the problem (2.1) where $\mathcal{L}$ and $\mathcal{B}$ are partial derivative operators and $f$ and $g$ are given functions.

$$ \begin{cases} \mathcal{L}u = f & \text{on } \Omega \\ \mathcal{B}u = g & \text{on } \partial\Omega \end{cases} \qquad (2.1) $$

For simplicity we will only consider a two-subdomains decomposition. Namely we have $\Omega = \Omega_1 \cup \Omega_2$. This decomposition can be overlapping or not. In the case of physical coupling, the subdomains do not overlap since the physics and thus the partial derivative operators are different on the two subdomains. We introduce the domains interface $\Gamma = \partial\Omega_1 \cap \partial\Omega_2$. We also define the external boundaries by $\partial\Omega_i^{\text{ext}} = \partial\Omega_i \setminus \Gamma$. In the most general case, this gives us two subproblems (2.2a) and (2.2b).

$$ \begin{cases} \mathcal{L}_1 u_1 = f_1 & \text{on } \Omega_1 \\ \mathcal{B}_1 u_1 = g_1 & \text{on } \partial\Omega_1^{\text{ext}} \\ \mathcal{I}_{1,1} u_1 = \mathcal{I}_{1,2} u_2 & \text{on } \Gamma \end{cases} \qquad (2.2a) $$

$$ \begin{cases} \mathcal{L}_2 u_2 = f_2 & \text{on } \Omega_2 \\ \mathcal{B}_2 u_2 = g_2 & \text{on } \partial\Omega_2^{\text{ext}} \\ \mathcal{I}_{2,2} u_2 = \mathcal{I}_{2,1} u_1 & \text{on } \Gamma \end{cases} \qquad (2.2b) $$

<figure><img src="image_7.png" /><figcaption>Figure 6: Non overlapping domain decomposition for coupling of physical models.</figcaption></figure>

$\mathcal{I}_{1,1}$, $\mathcal{I}_{1,2}$, $\mathcal{I}_{2,1}$, and $\mathcal{I}_{2,2}$ are partial derivative operators corresponding to interface boundary conditions. They can be of Dirichlet, Neumann or Robin type for instance. A basic example illustrating the domain decomposition for models coupling is presented on figure 6.
## 2.2 Algorithms

There are two versions of the Schwarz algorithm: the Schwarz sequential method and the Schwarz parallel method.

**Algorithm 1** Schwarz sequential method

Require: $u_2^0$ given on $\Gamma$, $N$
$n = 0$

while $n < N$ do

solve

$$ \begin{cases} \mathcal{L}_1 u_1^n = f_1 & \text{on } \Omega_1 \\ \mathcal{B}_1 u_1^n = g_1 & \text{on } \partial \Omega_1^{\text{ext}} \\ C_{1,1} u_1^n = C_{1,2} u_2^{n-1} & \text{on } \Gamma \end{cases} \quad (2.3a) $$

then solve

$$ \begin{cases} \mathcal{L}_2 u_2^n = f_2 & \text{on } \Omega_2 \\ \mathcal{B}_2 u_2^n = g_2 & \text{on } \partial \Omega_2^{\text{ext}} \\ C_{2,2} u_2^n = C_{2,1} u_1^n & \text{on } \Gamma \end{cases} \quad (2.3b) $$

end while

**Algorithm 2** Schwarz parallel method

Require: $u_1^0$ and $u_2^0$ given on $\Gamma$, $N$
$n = 0$

while $n < N$ do

solve

$$ \begin{cases} \mathcal{L}_1 u_1^n = f_1 & \text{on } \Omega_1 \\ \mathcal{B}_1 u_1^n = g_1 & \text{on } \partial \Omega_1^{\text{ext}} \\ C_{1,1} u_1^n = C_{1,2} u_2^{n-1} & \text{on } \Gamma \end{cases} \quad (2.4a) $$

and solve

$$ \begin{cases} \mathcal{L}_2 u_2^n = f_2 & \text{on } \Omega_2 \\ \mathcal{B}_2 u_2^n = g_2 & \text{on } \partial \Omega_2^{\text{ext}} \\ C_{2,2} u_2^n = C_{2,1} u_1^{n-1} & \text{on } \Gamma \end{cases} \quad (2.4b) $$

end while

$C_{1,1}$, $C_{1,2}$, $C_{2,1}$, and $C_{2,2}$ are partial derivative operators to be chosen. We see that the Schwarz parallel method is naturally parallel, hence the name. However, it usually requires a larger number of iterations (typically twice more) to converge compared to the Schwarz sequential method.

## 2.3 Convergence Analysis

In order to investigate the convergence speed of the Schwarz algorithm, it is useful to take a look at the errors evolution. Let us set $e_i^n = u_i^n - u_i$ and assume the partial derivative operators $\mathcal{L}_i$ and $\mathcal{B}_i$ are linear. For the sequential case, by subtracting (2.2a) to (2.3a) and (2.2b) to (2.3b), the errors satisfy

$$ \begin{cases} \mathcal{L}_1 e_1^n = 0 & \text{on } \Omega_1 \\ \mathcal{B}_1 e_1^n = 0 & \text{on } \partial \Omega_1^{\text{ext}} \\ C_{1,1} e_1^n = C_{1,2} e_2^{n-1} - (C_{1,1} - \mathcal{I}_{1,1}) u_1 + (C_{1,2} - \mathcal{I}_{1,2}) u_2 & \text{on } \Gamma \end{cases} $$

$$ \begin{cases} \mathcal{L}_2 e_2^n = 0 & \text{on } \Omega_2 \\ \mathcal{B}_2 e_2^n = 0 & \text{on } \partial \Omega_2^{\text{ext}} \\ C_{2,2} e_2^n = C_{2,1} e_1^n - (C_{2,2} - \mathcal{I}_{2,2}) u_2 + (C_{2,1} - \mathcal{I}_{2,1}) u_1 & \text{on } \Gamma \end{cases} $$

This points out the fact that we would get benefits from choosing $C = \mathcal{I}$.
## 2.4 Coupled Stationary 1D Diffusion Equation

For a better understanding of the Schwarz methods, let us provide a simple example of a coupling problem. Consider $\Omega = \Omega_o \cup \Omega_a$, with $\Omega_o = [H_o, 0]$ and $\Omega_a = [0, H_a]$. We want to solve the stationary coupled diffusion problem (2.6) with Neumann-Dirichlet interface conditions ensuring the continuity of temperature and the right exchange of heat fluxes.

$$ egin{cases} -K_o \frac{\partial^2 u_o}{\partial z^2} = f_o & \text{on } [H_o, 0] \\ u_o(H_o) = g_o(H_o) \\ k_o \frac{\partial u_o}{\partial z}(0) = k_a \frac{\partial u_a}{\partial z}(0) \end{cases} \quad \text{and} \quad 
\begin{cases} -K_a \frac{\partial^2 u_a}{\partial z^2} = f_a & \text{on } [0, H_a] \\ u_a(H_a) = g_a(H_a) \\ u_a(0) = u_o(0) \end{cases} \quad (2.6) $$

In this case, $\mathcal{L}_o = -K_o \frac{\partial^2}{\partial z^2}$, $\mathcal{L}_a = -K_a \frac{\partial^2}{\partial z^2}$, $\mathcal{B}_o = \mathcal{B}_a = \text{Id}$, $\mathcal{C}_{o,o} = k_o \frac{\partial}{\partial z}$, $\mathcal{C}_{o,a} = k_a \frac{\partial}{\partial z}$ and $\mathcal{C}_{a,a} = \mathcal{C}_{a,o} = \text{Id}$. If we apply the Schwarz sequential method, the errors $e_o^n$ and $e_a^n$ then satisfy

$$ egin{cases} -K_o \frac{\partial^2 e_o^n}{\partial z^2} = 0 & \text{on } \Omega_o \\ e_o^n(H_o) = 0 \\ k_o \frac{\partial e_o^n}{\partial z}(0) = k_a \frac{\partial e_a^{n-1}}{\partial z}(0) \end{cases} \quad \text{and} \quad 
\begin{cases} -K_a \frac{\partial^2 e_a^n}{\partial z^2} = 0 & \text{on } \Omega_a \\ e_a^n(H_a) = 0 \\ e_a^n(0) = e_o^n(0) \end{cases} $$

This means that the errors are of the following form

$$ e_o^n(z) = a_n(z - H_o) \quad \text{and} \quad e_a^n(z) = c_n(z - H_a) $$

Moreover, from the interface conditions we get

$$ a_n = \frac{k_a}{k_o} c_{n-1} \quad \text{and} \quad c_n = \frac{H_o}{H_a} a_n $$

We can deduce the rates of convergence of the Schwarz algorithm $\rho_o^n$ and $\rho_a^n$ in the two subdomains $\Omega_o$ and $\Omega_a$.

$$ \rho_a^n(z) = \rho_o^n(z) = \frac{|a_n|}{|a_{n-1}|} = \frac{k_a}{k_o} \frac{|H_o|}{H_a} = \rho^{\text{stat}} $$

We notice that the rate of convergence is assymetric. If we had imposed Dirichlet-Neumann interface conditions, we would have obtained the multiplicative inverse for the rate of convergence. A choice is thus necessary to ensure the convergence of the algorithm.

Let us check that we observe this result numerically. For physical values, we will consider the setting presented in table 1.

We further impose $\frac{\partial e_o^n}{\partial z}(0) = 0.1$ for the initialization. The observed errors and theoretical rate of convergence can be seen in figure 7, where we perform 15 iterations of the Schwarz sequential method.
<table><caption>**Table 1:** Physical settings in the ocean and atmosphere.</caption><thead><tr><th>Physical constant</th><th>Symbol</th><th>Unit</th><th>Ocean</th><th>Atmosphere</th></tr></thead><tbody><tr><td>Height</td><td><i>H</i></td><td>m</td><td>-50</td><td>200</td></tr><tr><td>Thermal conductivity</td><td><i>k</i></td><td>J.s<sup>-1</sup>.m<sup>-1</sup>.K<sup>-1</sup></td><td>0.570</td><td>0.02364</td></tr><tr><td>Thermal diffusivity</td><td><i>K</i></td><td>m<sup>2</sup>.s<sup>-1</sup></td><td>1.39 &times; 10<sup>-7</sup></td><td>1.82 &times; 10<sup>-5</sup></td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Figure 7:** Observed errors (blue and orange lines) and theoretical convergence rate (dashed green lines) of the Schwarz sequential method for the coupled 1D stationary diffusion problem. x-axis: iteration number. y-axis: error in $L^2$ norm.</figcaption></figure>
## 2.5 Coupled Nonstationary 1D Diffusion Equation

For the nonstationary case, the partial derivative operators $\mathcal{L}_o$ and $\mathcal{L}_a$ simply become
$\mathcal{L}_o = \frac{\partial}{\partial t} - K_o \frac{\partial^2}{\partial z^2}$ and $\mathcal{L}_a = \frac{\partial}{\partial t} - K_a \frac{\partial^2}{\partial z^2}$. The evolution of the errors using the Schwarz sequential method is now given by

$$ 
\left\{ 
\begin{array}{ll} 
\displaystyle \frac{\partial e_o^n}{\partial t} - K_o \frac{\partial^2 e_o^n}{\partial z^2} = 0 & \text{on } [H_o, 0] \times ]0, +\infty[ \\ \\ 
\displaystyle e_o^n(z, 0) = 0 & \text{on } [H_o, 0] \\ \\ 
\displaystyle e_o^n(H_o, t) = 0 & \text{on } ]0, +\infty[ \\ \\ 
\displaystyle k_o \frac{\partial e_o^n}{\partial z}(0, t) = k_a \frac{\partial e_a^{n-1}}{\partial z}(0, t) & \text{on } ]0, +\infty[ 
\end{array} 
\right. 
\quad \text{and} 
\quad 
\left\{ 
\begin{array}{ll} 
\displaystyle \frac{\partial e_a^n}{\partial t} - K_a \frac{\partial^2 e_a^n}{\partial z^2} = 0 & \text{on } ]0, H_a] \times ]0, +\infty[ \\ \\ 
\displaystyle e_a^n(z, 0) = 0 & \text{on } ]0, H_a] \\ \\ 
\displaystyle e_a^n(H_a, t) = 0 & \text{on } ]0, +\infty[ \\ \\ 
\displaystyle e_a^n(0, t) = e_o^n(0, t) & \text{on } ]0, +\infty[ 
\end{array} 
\right. 
\quad (2.7) 
$$

Let us define the Fourier transform with respect to $t$ by

$$ \hat{w}(z, \omega) = \frac{1}{2\pi} \int_{\mathbb{R}} w(z, t) e^{-i\omega t} dt $$

In the spectral domain, (2.7) then becomes

$$ 
\left\{ 
\begin{array}{ll} 
\displaystyle \left( i\omega - K_o \frac{\partial^2}{\partial z^2} \right) \hat{e}_o^n = 0 & \text{on } \Omega_o \\ \\ 
\displaystyle \hat{e}_o^n = 0 & \text{on } \partial\Omega_o^{\text{ext}} \quad \text{and} \\ \\ 
\displaystyle k_o \frac{\partial \hat{e}_o^n}{\partial z} = k_a \frac{\partial \widehat{e_a^{n-1}}}{\partial z} & \text{on } \Gamma 
\end{array} 
\right. 
\quad \text{and} 
\quad 
\left\{ 
\begin{array}{ll} 
\displaystyle \left( i\omega - K_a \frac{\partial^2}{\partial z^2} \right) \hat{e}_a^n = 0 & \text{on } \Omega_a \\ \\ 
\displaystyle \hat{e}_a^n = 0 & \text{on } \partial\Omega_a^{\text{ext}} \\ \\ 
\displaystyle \hat{e}_a^n = \hat{e}_o^n & \text{on } \Gamma 
\end{array} 
\right. 
$$

This means that the errors are of the following form

$$ \hat{e}_o^n(z, \omega) = A_n(\omega) \sinh \left( \sqrt{\frac{i\omega}{K_o}} (z - H_o) \right) \quad \text{and} \quad \hat{e}_a^n(z, \omega) = C_n(\omega) \sinh \left( \sqrt{\frac{i\omega}{K_a}} (z - H_a) \right) $$

This time, we get from the interface conditions that

$$ k_o A_n(\omega) \sqrt{\frac{i\omega}{K_o}} \cosh \left( \sqrt{\frac{i\omega}{K_o}} H_o \right) = k_a C_{n-1}(\omega) \sqrt{\frac{i\omega}{K_a}} \cosh \left( \sqrt{\frac{i\omega}{K_a}} H_a \right) $$

and

$$ A_n(\omega) \sinh \left( \sqrt{\frac{i\omega}{K_o}} H_o \right) = C_n(\omega) \sinh \left( \sqrt{\frac{i\omega}{K_a}} H_a \right) $$
So by taking the ratio of the two equations, we deduce that

$$ \frac{A_n(\omega)}{A_{n-1}(\omega)} = \frac{C_n(\omega)}{C_{n-1}(\omega)} = \frac{k_a \sqrt{\frac{i\omega}{K_a}} \cosh\left(\sqrt{\frac{i\omega}{K_a}} H_a\right)}{k_o \sqrt{\frac{i\omega}{K_o}} \cosh\left(\sqrt{\frac{i\omega}{K_o}} H_o\right)} \frac{\sinh\left(\sqrt{\frac{i\omega}{K_o}} H_o\right)}{\sinh\left(\sqrt{\frac{i\omega}{K_o}} H_a\right)} \\ = \frac{k_a}{k_o} \sqrt{\frac{K_o}{K_a}} \frac{\tanh\left(\sqrt{\frac{i\omega}{K_o}} H_o\right)}{\tanh\left(\sqrt{\frac{i\omega}{K_o}} H_a\right)} $$

This means that the rates of convergence in $\Omega_o$ and $\Omega_a$ are given by

$$ \rho_a^n(z, \omega) = \rho_o^n(z, \omega) = \frac{|C_n(\omega)|}{|C_{n-1}(\omega)|} = \frac{k_a}{k_o} \sqrt{\frac{K_o}{K_a}} \left| \frac{\tanh\left(\sqrt{\frac{i\omega}{K_o}} H_o\right)}{\tanh\left(\sqrt{\frac{i\omega}{K_o}} H_a\right)} \right| = \rho(\omega) \quad (2.8) $$

**Asymptotic cases** However, we know that for a small enough $x$

$$ \tanh(x) = x + o(x) $$

So for low frequencies, we deduce that

$$ \lim_{\omega \to 0} \rho(\omega) = \frac{k_a}{k_o} \sqrt{\frac{K_o}{K_a}} \frac{\sqrt{\frac{i\omega}{K_o}} |H_o|}{\sqrt{\frac{i\omega}{K_o}} H_a} = \frac{k_a |H_o|}{k_o H_a} = \rho^{\text{stat}} $$

It is logical to retrieve the rate of convergence of the stationary case since low frequencies mean slow changes in time. Moreover for high frequencies, we have

$$ \lim_{\omega \to +\infty} \rho(\omega) = \frac{k_a}{k_o} \sqrt{\frac{K_o}{K_a}} = \rho^{\text{nonstat}} $$

This can be observed in figure 8 where we plotted $\rho(\omega)$ for the setting of constant values we defined earlier. The numerically obtained errors are also shown in the same figure, with 15 iterations. We observe that the convergence rate is given by the high frequency asymptote of $\rho(\omega)$.

## 3 Atmosphere-Ice-Ocean Model

Now that we introduced the main algorithm for solving coupled problems, we propose to model interactions between atmosphere, sea ice and ocean in this section. The study of this three media coupling is interesting since adding a layer of sea ice means reducing the direct interactions between the atmosphere and the ocean. Indeed, the ocean has a huge ability to absorb heat and it is thus naturally a climate regulator. Besides, the presence of ice increases the amount of reflected solar energy back in the atmosphere. This is why it is important to succeed in modeling sea ice.
<figure><img src="image_9.png" /><figcaption>(a) Theoretical rate of convergence (solid line) with low and high frequencies asymptotes (dashed green and red lines)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(b) Observed errors (solid blue and orange lines) between minimum (brown dashdots) and maximum (pink dashdots) theoretical rate of convergence</figcaption></figure>

**Figure 8:** Theoretical convergence rate and observed errors of the Schwarz sequential method for the coupled 1D nonstationary diffusion problem.

## 3.1 Motivations

For the moment, the algorithms used in climate models to solve the coupling problems actually correspond to a single iteration of the Schwarz algorithm. What has been shown is that performing only two iterations would reduce the errors made, especially the ones linked to the diurnal cycle as it can be observed on figure 9 from [1].

The authors took as reference solution the fully converged solution of the Schwarz algorithm. However performing two iterations means double the computational cost. This is clearly not acceptable. Besides, it has been observed that more than two iterations are needed for the algorithm to converge in sea-ice regions as figure 10a shows. Moreover for some icy points in the globe, strong oscillations appear in the convergence of the ice surface temperature. This phenomenon is shown in figure 10b where the amplitude of oscillations is measured.

To investigate these two convergence issues, we will derive a simplified unidimensional vertical model *hopefully* retaining the slow convergence phenomena. Analysis of the Schwarz algorithm will be conducted on this simpler model.

## 3.2 The Model

In this section, we establish an air-ice-ocean model by adapting the work of MAYKUT & UNTERSTEINER in [2] and SEMTNER in [3]. The sea ice is seen as a flat slab above the ocean surface. The model principle is based on the physical balance between the different heat fluxes involved at both upper and lower ice interface. We take different fluxes into account to establish our model, as illustrated in figure 11. From this energy balance, the rate of accretion or ablation of the ice over time is deduced. We suppose
<figure><img src="image_11.png" /><figcaption>**Figure 9:** Histograms of errors as a function of the Roman local time and error classes for two iterations of the Schwarz algorithm. Those results were obtained with IPSL-CM.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10:** Convergence issues in sea-ice regions.</figcaption></figure>

Credits: MARTI Olivier (Laboratoire des Sciences du Climat et de l'Environnement, Gif-sur-Yvette)
that the whole portion of ocean is covered by sea ice, there is no fraction of ice. So ocean and atmosphere do not share any boundary.

**Net surface heat flux** At the surface of the ice, there is the net surface heat flux $F_{net}^{ice}$ which is the sum of several fluxes.

$$F_{net}^{ice} = \underbrace{SW_{net}}_{\text{solar}} + \underbrace{F_{lat} + F_{sen} + LW_{net}}_{\text{non-solar}} \quad (3.1)$$

where $SW_{net}$ is the net surface shortwave radiative flux

$$SW_{net} = (1 - \alpha_{ice})SW^{in} \quad (3.2)$$

with $\alpha_{ice} \simeq 0.8$ the ice surface albedo which will reflect a portion of the sun incoming shortwave radiation $SW^{in}$.

**Latent and sensible heat flux** The latent heat flux $F_{lat}$ (associated to vaporization of water and fusion of ice) is not taken

into account in this model. $F_{sen}$ is the sensible heat flux (linked to the difference of temperature between the ice surface and the atmosphere)

<figure><img src="image_13.png" /><figcaption>**Figure 11:** Schematic of the atmosphere-ice-ocean model.</figcaption></figure>

$$\begin{aligned} F_{sen} &= \rho_a c_a C_{sen} W_a (T_a - T_{su}) \\ &= \alpha_a (T_a - T_{su}) \end{aligned} \quad (3.3)$$

where $\rho_a = 1.292$ is the air density (in kg.m$^{-3}$), $c_a = 1006$ is the air specific heat capacity (in J.kg$^{-1}$.K$^{-1}$), $C_{sen} = 1.4 \times 10^{-3}$ is the heat transfer coefficient for sensible heat, $W_a$ is the scalar wind speed at the bottom of the atmosphere (in m.s$^{-1}$) and $T_a$ is the temperature in the atmosphere just above the ice surface.

**Longwave radiative flux** $LW_{net}$ is the net surface longwave radiative flux: the sum of incoming longwave radiation and black-body radiation (Stefan-Boltzmann law).

$$LW_{net} = \epsilon(LW^{in} - \sigma_{SB}T_{su,K}^{4}) \quad (3.4)$$

where $\epsilon \simeq 0.98$ is the emissivity of the surface, $\sigma_{SB} = 5.67 \times 10^{-8}$ is the Stefan-Boltzmann constant (in J.s$^{-1}$.m$^{-2}$.K$^{-4}$), $LW^{in}$ is the incoming longwave radiation (that we will neglect) and $T_{su,K}$ is the surface temperature in Kelvin. Note that for temperatures between -40°C and 0°C, typically the ones we are working with, the outgoing longwave radiation can be linearized [4]. The approximation can be seen on figure 12.

$$\sigma_{SB}T_{su,K}^{4} \simeq A + BT_{su} \quad (3.5)$$
where $A = 309.7913$ (in J.s$^{-1}$.m$^{-2}$) and $B = 3.6939$ (in J.s$^{-1}$.m$^{-2}$.°C$^{-1}$). In the end, if we denote the ice thickness by $h(t)$, the net surface flux reads

$$F_{net}^{ice}(t) = SW_{net} + \alpha_a(T_a(z = h(t), t) - T_{su}(t)) - \epsilon(A + BT_{su}(t)) \quad (3.6)$$

**Ice inner conductive flux** We assume that the vertical gradient of temperature is constant in the ice and that the ice bottom temperature is at freezing point $T_b = -1.8°C$. Then the inner conductive heat flux in the ice $F_c$ is given by

$$\begin{aligned} F_c(z,t) &= -k_i \frac{\partial T}{\partial z}(z,t) \\ &= -\frac{k_i}{h(t)}(T_{su}(t) - T_b) = F_c(t) \end{aligned} \quad (3.7)$$

with $k_i = 2.0344$ the thermal conductivity of ice (in J.s$^{-1}$.m$^{-1}$.K$^{-1}$).

### 3.2.1 A Balance Between Heat Fluxes

We have two heat flux interface conditions. At the top of the ice, there is a first energy balance between the net surface heat flux and the ice inner conductive flux.

<figure><img src="image_14.png" /><figcaption>**Figure 12:** Linear approximation of the outgoing longwave radiation.</figcaption></figure>

$$F_{net}^{ice} + F_c = 0 \quad (3.8)$$

Using (3.6) and (3.7), we deduce an expression for the ice surface temperature which has to be below the melting point of water $T_f$.

$$T_{su}(t) = \min \left\{ \frac{SW_{net}(t) - \epsilon A + \alpha_a T_a(h(t), t) + \frac{k_i}{h(t)} T_b}{\frac{k_i}{h(t)} + \epsilon B + \alpha_a}, T_f \right\} \quad (3.9)$$

In the case where $T_{su}(t) = T_f$, ice is melting from the top using the surplus of heat flux.

$$-\rho_i L_{fus} \frac{dh}{dt} = F_{net}^{ice}(T_{su} = T_f) + F_c(T_{su} = T_f) \neq 0 \quad (3.10)$$

where $\rho_i = 917$ is the ice density (in kg.m$^{-3}$), $L_{fus} = 3.0033 \times 10^5$ is the latent heat of fusion of sea ice (in J.kg$^{-1}$) and $h$ is the ice thickness. At the bottom of the ice, either ablation or accretion can occur depending on the ocean surface temperature below the ice.

$$-\rho_i L_{fus} \frac{dh}{dt} = F_{oce} - F_c \quad (3.11)$$
where $F_{oce}$ is the ocean heat flux.

$$ 
\begin{aligned}
F_{oce}(T_o) &= \rho_o c_o C_h u^*(T_o - T_b) \\
&= \alpha_o (T_o - T_b)
\end{aligned}
\qquad (3.12) $$

where $\rho_o = 1024.458$ is the ocean density (in kg.m$^{-3}$), $c_o = 4001.896$ is the ocean specific heat capacity (in J.kg$^{-1}$.K$^{-1}$), $C_h = 7.5 \times 10^{-3}$ is the heat transfer coefficient (or turbulent Stanton number), $u^*$ is the interface friction velocity (in m.s$^{-1}$) and $T_o$ is the temperature at the top of the ocean. The general ice thickness evolution is given by

$$ -\rho_i L_{fus} \frac{dh}{dt} = F_{net}^{ice} + F_{oce} \qquad (3.13) $$

### 3.2.2 The System of Equations

In this section we combine the diffusion equations in the atmosphere and the ocean with the previous equations on the ice surface temperature and the ice thickness. We do not consider any diffusion equation in the ice. Ice rather plays the part of a buffer interface between the atmosphere and the ocean. This has an impact on the ocean which now stands alone as it always sees ice above at a constant temperature.

As an interface condition between atmosphere and ice, we impose the net surface heat flux $F_{net}^{ice}$ to be equal to the atmosphere heat flux. This results in a Robin type boundary condition at $z = h(t)$.

$$ 
\left\{
\begin{array}{ll}
\displaystyle \frac{\partial T_a}{\partial t}(z,t) - K_a \frac{\partial^2 T_a}{\partial z^2}(z,t) = 0 & \text{on } (h(t), H_a) \times (0, T) \\[2ex]
T_a(H_a, t) = g_a(t) & \text{on } (0, T) \\[2ex]
\displaystyle k_a \frac{\partial T_a}{\partial z}(h(t), t) - \alpha_a T_a(h(t), t) = SW_{net}(t) - \epsilon A - (\alpha_a + \epsilon B) T_{su}(t) & \text{on } (0, T) \\[2ex]
T_a(z, 0) \text{ given} & \text{on } (h(t), H_a)
\end{array}
\right.
\qquad (3.14) $$

Concerning the interface between ocean and ice, we impose the ocean heat flux to be equal to $F_{oce}$, which also gives a boundary condition of Robin type at $z = 0$.

$$ 
\left\{
\begin{array}{ll}
\displaystyle \frac{\partial T_o}{\partial t}(z,t) - K_o \frac{\partial^2 T_o}{\partial z^2}(z,t) = 0 & \text{on } (H_o, 0) \times (0, T) \\[2ex]
T_o(H_o, t) = g_o(t) & \text{on } (0, T) \\[2ex]
\displaystyle k_o \frac{\partial T_o}{\partial z}(0, t) - \alpha_o T_o(0, t) = -\alpha_o T_b & \text{on } (0, T) \\[2ex]
T_o(z, 0) \text{ given} & \text{on } (H_o, 0)
\end{array}
\right.
\qquad (3.15) $$
For the ice thickness evolution, we rewrite (3.13).

$$
\begin{cases}
-\rho_i L_{fus} \frac{dh}{dt}(t) = SW_{net}(t) - \epsilon(A + BT_{su}(t)) + \alpha_a(T_a(h(t),t) - T_{su}(t)) + \alpha_o(T_o(0,t) - T_b) & \text{on } (0,T) \\
h(0) \text{ given}
\end{cases}
\tag{3.16}
$$

with

$$
T_{su}(t) = \min \left\{ \frac{SW_{net}(t) - \epsilon A + \alpha_a T_a(h(t),t) + \frac{k_i}{h(t)} T_b}{\frac{k_i}{h(t)} + \epsilon B + \alpha_a}, T_f \right\}
$$

We observe that the ocean is isolated from the two other domains. It evolves on its own.
On the other hand, the equations seem difficult to solve: they are non linear and rather
complex. In section 3.3, we will try to solve them or at least a simplified version of them.

3.3 Solving the Ice Thickness Evolution Equation

In this section we focus on (3.16) and try to give some theoretical results on its solutions.
Since the min in the computation of the ice surface temperature is inconvenient, we
assume that the ice surface temperature is always below the melting point of water
$T_f$. This is not a strong assumption if we consider a cold season where the atmosphere
temperature is low. Then (3.8) is verified: there is no melting of the ice at the top. The
ice thickness evolution equation (3.13) becomes

$$
\begin{align}
-\rho_i L_{fus} \frac{dh}{dt}(t) &= F_{oce} - F_c \nonumber \\
&= \alpha_o (T_o(0,t) - T_b) + \frac{k_i}{h(t)} (T_{su}(t) - T_b) \nonumber \\
&= \alpha_o (T_o(0,t) - T_b) + \frac{SW_{net}(t) - \epsilon A + \alpha_a T_a(h(t),t) - (\epsilon B + \alpha_a)T_b}{1 + \frac{\epsilon B + \alpha_a}{k_i} h(t)} \tag{3.17}
\end{align}
$$

If we set $\beta = \frac{\epsilon B + \alpha_a}{k_i}$ and $u(t) = 1 + \beta h(t)$, then (3.17) turns into

$$
\frac{du}{dt}(t) + \frac{a(t)}{u(t)} = b(t) \tag{3.18}
$$

with

$$
a(t) = \frac{\beta}{\rho_i L_{fus}} (SW_{net}(t) - \epsilon A + \alpha_a T_a(h(t),t) - k_i \beta T_b)
$$

$$
b(t) = \frac{-\beta}{\rho_i L_{fus}} \alpha_o (T_o(0,t) - T_b)
$$

This type of non-linear differential equation is called Chini’s Equation. Unfortunately,
we only know solutions of these equations in specific cases: when *a* and *b* are constant
<figure><img src="image_15.png" /><figcaption>Figure 13: Lambert W-function. It verifies $z = W(z)e^{W(z)}$ for all $z \ge -e^{-1}$</figcaption></figure>

for instance. The solutions are linked to the Lambert W-function which is plotted on figure 13. More details are provided in appendix B. As we do not have a general solution for this type of equation and that we cannot make strong assumptions on the functions a and b, we have to make another assumption so as to analyse the Schwarz algorithm. On the other hand, the system of equation in the atmosphere becomes

$$ 
\left\{
\begin{array}{ll}
\displaystyle \frac{\partial T_a}{\partial t}(z,t) - K_a \frac{\partial^2 T_a}{\partial z^2}(z,t) = 0 & \text{on } (h(t), H_a) \times (0, T) \\[1.5ex]
T_a(H_a, t) = g_a(t) & \text{on } (0, T) \\[1.5ex]
\displaystyle \frac{\partial T_a}{\partial z}(h(t), t) - \frac{\alpha_a}{k_a} T_a(h(t), t) = & \displaystyle -\frac{\alpha_a}{k_a} \frac{u(t) - 1}{u(t)} T_a(h(t), t) \\
& + \frac{SW_{net}(t) - \epsilon A - k_i \beta T_b}{k_a u(t)} & \text{on } (0, T) \\[1.5ex]
T_a(z, 0) \text{ given} & \text{on } (h(t), H_a)
\end{array}
\right.
\qquad (3.19)
$$

## 4 Applying the Schwarz Iterative Method

Let us now propose a Schwarz algorithm to solve the coupled problem presented in subsection 3.2.2, under the assumption of section 3.3, which is having an ice surface temperature below the melting point of water.

The iterative algorithm is as follows: assume we are given a first guess atmosphere temperature at the bottom $T_a^0(h(t), t)$. We iteratively solve the ice thickness evolution
equation giving the rate of ablation or accretion of ice

$$ \frac{du^n}{dt}(t) + \frac{a^{n-1}(t)}{u^n(t)} = b(t) \qquad (4.1) $$

with

$$ a^{n-1}(t) = \frac{\beta}{\rho_i L_{fus}} (SW_{net}(t) - \epsilon A + \alpha_a T_a^{n-1}(h(t), t) - k_i \beta T_b) $$

Then we solve the diffusion equation in the atmosphere

$$ \left\{ \begin{array}{ll} \left(\frac{\partial}{\partial t} - K_a \frac{\partial^2}{\partial z^2}\right) T_a^n(z,t) = 0 & \text{on } (h(t), H_a) \times (0, T) \\ \\ T_a^n(H_a, t) = g_a(t) & \text{on } (0, T) \\ \\ \left(\frac{\partial}{\partial z} - \frac{\alpha_a}{k_a} Id\right) T_a^n(h(t), t) = & -\frac{\alpha_a}{k_a} \frac{u^n(t) - 1}{u^n(t)} T_a^{n-1}(h(t), t) \\ & + \frac{SW_{net}(t) - \epsilon A - k_i \beta T_b}{k_a u^n(t)} \quad \text{on } (0, T) \\ \\ T_a^n(z, 0) \quad \text{given} & \text{on } (h(t), H_a) \end{array} \right. \qquad (4.2) $$

These equations are non linear and we thus cannot express the errors in the two media. We have to linearize the operators.

## 4.1 Further Assumptions

We make several approximations for both numerical and analytical purposes.

* Since the ice thickness h is negligible when compared to the atmosphere height $H_a$, we will rather solve (3.14) on $(0, H_a)$.

* We further assume that the ice thickness is almost a constant on the Schwarz time window $(0, T)$. Namely

$$ h(t) = h_0 + \delta h(t) \quad \text{with} \quad h_0 = h(0) \quad \text{and} \quad \delta h(t) \ll h_0 $$

which gives in terms of u:

$$ u(t) = u_0 + \delta u(t) \quad \text{with} \quad u_0 = 1 + \beta h_0 \quad \text{and} \quad \delta u(t) \ll u_0 $$

This allows us to linearize (3.18) and apply the analysis methodology we presented in section 2 for which linearity of partial derivative operators was a limiting assumption.
## 4.2 Iterative Algorithm

Assume we are given a first guess atmosphere temperature at the bottom $T_a^0(0, t)$. We iteratively solve the ice thickness evolution equation giving the rate of ablation or accretion of ice

$$ \frac{du^n}{dt}(t) + \frac{a^{n-1}(t)}{u_0} = b(t) \qquad (4.3) $$

with

$$ a^{n-1}(t) = \frac{\beta}{\rho_i L_{fus}} (SW_{net}(t) - \epsilon A + \alpha_a T_a^{n-1}(0, t) - k_i \beta T_b) $$

And we solve the diffusion equation in the atmosphere

$$ \left\{ \begin{array}{ll} \left(\frac{\partial}{\partial t} - K_a \frac{\partial^2}{\partial z^2}\right) T_a^n(z, t) = 0 & \text{on } (0, H_a) \times (0, T) \\ \\ T_a^n(H_a, t) = g_a(t) & \text{on } (0, T) \\ \\ \left(\frac{\partial}{\partial z} - \frac{\alpha_a}{k_a} Id\right) T_a^n(0, t) = -\frac{\alpha_a}{k_a} \frac{u_0 - 1}{u_0} T_a^{n-1}(0, t) + \frac{SW_{net}(t) - \epsilon A - k_i \beta T_b}{k_a u_0} & \text{on } (0, T) \\ \\ T_a^n(z, 0) \quad \text{given} & \text{on } (0, H_a) \end{array} \right. \qquad (4.4) $$

## 4.3 Errors and Convergence Analysis

The following step in the analysis of the Schwarz algorithms is to focus on the errors. We consider here the errors $e_a^n = T_a^n - T_a$ and $e_u^n = u^n - u$. Errors verify the same type of equations as the quantity to which they are related. Namely, a simple differential equation for $e_u^n$.

$$ \frac{de_u^n}{dt}(t) + \frac{\alpha_a \beta}{u_0 \rho_i L_{fus}} e_a^{n-1}(0, t) = 0 \quad \text{on } (0, T) \qquad (4.5) $$

and a homogeneous diffusion equation for $e_a^n$

$$ \left\{ \begin{array}{ll} \left(\frac{\partial}{\partial t} - K_a \frac{\partial^2}{\partial z^2}\right) e_a^n(z, t) = 0 & \text{on } (0, H_a) \times (0, T) \\ \\ e_a^n(H_a, t) = 0 & \text{on } (0, T) \\ \\ \left(\frac{\partial}{\partial z} - \frac{\alpha_a}{k_a} Id\right) e_a^n(0, t) = -\frac{\alpha_a}{k_a} \frac{u_0 - 1}{u_0} e_a^{n-1}(0, t) & \text{on } (0, T) \\ \\ e_a^n(z, 0) = 0 & \text{on } (0, H_a) \end{array} \right. \qquad (4.6) $$

In order to transform the PDE in ODE, we will work in the spectral domain. We denote by $\hat{f}$ the Fourier transform of a function $f$ with respect to the time. From (4.6), we have

$$ \hat{e}_a^n(z, \omega) = A_n(\omega) \sinh\left(\sqrt{\frac{i\omega}{K_a}}(z - H_a)\right) $$
With the interface condition, we deduce the rate of convergence in the atmosphere. This is the same convergence rate for the ice thickness since the errors are linearly linked in the spectral domain.

$$ \rho_a^n(z, \omega) = \frac{|A_n(\omega)|}{|A_{n-1}(\omega)|} = \alpha_a \frac{u_0 - 1}{u_0} \frac{|\sinh(\sqrt{\frac{i\omega}{K_a}} H_a)|}{|\alpha_a \sinh(\sqrt{\frac{i\omega}{K_a}} H_a) + k_a \sqrt{\frac{i\omega}{K_a}} \cosh(\sqrt{\frac{i\omega}{K_a}} H_a)|} \quad (4.7) $$

In a slightly simpler form, the rate of convergence is given by

$$ \rho_a^n(z, \omega) = \alpha_a \frac{u_0 - 1}{u_0} \frac{|\tanh(\sqrt{\frac{i\omega}{K_a}} H_a)|}{|\alpha_a \tanh(\sqrt{\frac{i\omega}{K_a}} H_a) + k_a \sqrt{\frac{i\omega}{K_a}}|} = \rho^{\text{AIO}}(\omega, h_0) \quad (4.8) $$

**Asymptotic cases** So for low frequencies, we deduce that

$$ \lim_{\omega \to 0} \rho^{\text{AIO}}(\omega, h_0) = \frac{u_0 - 1}{u_0} \frac{\alpha_a H_a}{\alpha_a H_a + k_a} $$

We notice that when the ice thickness tends toward zero, the convergence rate vanishes. We can conclude that the ocean-atmosphere model without ice is not a limit case of this three media model. Moreover for high frequencies, we have

$$ \lim_{\omega \to +\infty} \rho^{\text{AIO}}(\omega, h_0) = 0 $$

We can compare the rate of convergence in absence of ice (2.8) with the rates of convergence (4.8) for different ice thickness $h_0$. This is what is shown on figure 14a. Moreover, the rate of convergence in presence of ice is growing with the ice thickness. It can be observed on figure 14b that when the ice thickness is larger than 0.2 meter, the Schwarz algorithm converges very slowly. This is in keeping with what has been observed in large scale simulation. Indeed, the slow convergence appears in dense ice regions.

# 5 Numerical Simulations

We will now provide some numerical results to validate the theoretical rate of convergence we derived earlier.

## 5.1 Physical Settings

We work with the setting presented in table 1 on a time window of 12 hours divided in 24 time steps. The incoming shortwave radiation during this time window is shown in figure 15a. As for the first guess of ice surface temperature needed to start the Schwarz algorithm, our choice is presented in figure 15b.

For the ocean and atmosphere temperature profile at initial time $t = 0$, we choose two constant profiles

$$ T_o(z, 0) = T_b \quad \text{and} \quad T_a^n(z, 0) = \frac{u_0 - 1}{u_0} T_a^{n-1}(0, 0) - \frac{SW_{net}(0) - \epsilon A - k_i \beta T_b}{\alpha_a u_0} $$
<figure><img src="image_16.png" /><figcaption>(a) Comparison of the convergence rates of the ocean-atmosphere coupling and ocean-ice-atmosphere coupling</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(b) Theoretical convergence rate of Schwarz algorithm in the atmosphere with respect to the ice thickness</figcaption></figure>

**Figure 14:** An overview of the theoretical rate of convergence of the ocean-ice-atmosphere coupling model.

<figure><img src="image_18.png" /><figcaption>(a) Incoming shortwave radiation $SW^{in}$</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(b) Ice surface temperature first guess $T_{su}^0$</figcaption></figure>

**Figure 15:** Incoming shortwave radiation and first guess for ice surface temperature as a function of time.
## 5.2 Results

Since we do not have the exact solution, we cannot observe the error numerically. However, we can have a look to the difference between solutions at two successive iterations. This is what we call *observed convergence rate* in figure 16. Note that this is consistent to compare those two quantities in logarithmic scale. Indeed if the error at iteration *n* is given by $e^n = \rho^n e^0$, we have

$$ \log(e^n) = n \log(\rho) + \log(e^0) $$

and

$$ \log(u^n - u^{n-1}) = \log(e^n - e^{n-1}) \\ = (n-1) \log(\rho) + \log((\rho-1)e^0) $$

Since we are only interested in the slope, it is justified to compare the two *errors*.

<figure><img src="image_20.png" /><figcaption>**Figure 16:** Theoretical (dashed lines) and observed (solid lines) convergence rate of Schwarz algorithm for 3 different ice thickness. The theoretical rate of convergence is $\rho^{AIO}(0, h_0)$.</figcaption></figure>

The numerically observed convergence rates are rather consistent with the computations we derived. Nevertheless, it can be observed that the theoretical rate of convergence we computed is *optimistic* since the solutions seem to be converging more slowly.

## 5.3 More General Case

In the case we make no assumptions on the ice surface temperature nor the behaviour of the ice thickness, we cannot derive the rate of convergence. But we can still run
simulations. The ice surface temperature first guess $T_{su}^0$ is now represented in figure 17a. As shown in figure 17b, the observed convergence rate is close to our previous computations. This is probably due to the fact that the ice surface temperature converges to the constant $T_b$. Thus our assumptions are valid in the general case too.

<figure><img src="image_21.png" /><figcaption>(a) Ice surface temperature first guess as a function of time</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>(b) Theoretical (dash lines) and observed (solid lines) convergence rate of Schwarz algorithm for 3 different initial ice thickness</figcaption></figure>

**Figure 17:** Setting and results obtained in the general case with no approximation.

# 6 Concerning Oscillations

The second issue we observed during the iterative process in large scale models was strong oscillations in the convergence. As shown in figure 18, there is a kind of 4 iterations oscillation pattern. To explain these oscillations, we try to provide an analytic formulation of the Schwarz algorithm error and see if there is an oscillating term in the expression. So we want to solve a nonstationary homogeneous diffusion problem with inhomogeneous boundary conditions (4.6). To do so, we first transform this problem into one with homogeneous boundary conditions using the Shifting Method. We split the error $e_a^n$ into two parts $e_a^n(z, t) = v^n(z, t) + w^n(z, t)$

<figure><img src="image_23.png" /><figcaption>**Figure 18:** Oscillations in the temperature at a specific point. From an IPSL-CM simulation.</figcaption></figure>
with the functions $v^n$ and $w^n$ verifying

$$ 
\begin{aligned}
w^n(z,t) &= - \left(1 + \frac{\alpha_a}{k_a} H_a\right)^{-1} \frac{\alpha_a}{k_a} \frac{u_0 - 1}{u_0} (z - H_a) e_a^{n-1}(0,t) \\
&= C(z - H_a) e_a^{n-1}(0,t)
\end{aligned} 
\qquad (6.1) $$

i.e. $w^n$ satisfies the boundary conditions

$$ 
\left\{ 
\begin{array}{ll} 
\left(\frac{\partial}{\partial z} - \frac{\alpha_a}{k_a} Id\right) w^n(0,t) = -\frac{\alpha_a}{k_a} \frac{u_0 - 1}{u_0} e_a^{n-1}(0,t) & \text{on } (0,T) \\ 
w^n(H_a,t) = 0 & \text{on } (0,T) 
\end{array} 
\right. 
$$

and

$$ 
\left\{ 
\begin{array}{ll} 
\left(\frac{\partial}{\partial t} - K_a \frac{\partial^2}{\partial z^2}\right) v^n(z,t) = -C(z - H_a) \frac{\partial e_a^{n-1}}{\partial t}(0,t) & \text{on } (0, H_a) \times (0, T) \\ 
v^n(H_a, t) = 0 & \text{on } (0, T) \\ 
\left(\frac{\partial}{\partial z} - \frac{\alpha_a}{k_a} Id\right) v^n(0,t) = 0 & \text{on } (0, T) \\ 
v^n(z, 0) = 0 & \text{on } (0, H_a) 
\end{array} 
\right. 
\qquad (6.2) $$

Then we can rely on the Sturm-Liouville Theory to express the solution of this new problem. More details on the theory are given in appendix A.

$$ v^n(z,t) = \sum_{k=0}^{+\infty} Z_k(z)T_k(t) $$

with

$$ Z_k(z) = \sin(\mu_k z) + \mu_k \frac{k_a}{\alpha_a} \cos(\mu_k z) \quad \text{and} \quad T_k(t) = \int_0^t f_k(s) e^{-K_a \mu_k^2 (t-s)} ds $$

where

$$ f_k(s) = -C \frac{\int_0^{H_a} (z - H_a) Z_k(z) dz}{\langle Z_k, Z_k \rangle} \frac{\partial e_a^{n-1}}{\partial t}(0, s) $$

and $\mu_k$ solves

$$ \mu + \frac{\alpha_a}{k_a} \tan(\mu H_a) = 0 $$

We cannot really draw conclusions on the oscillations from this expression. There is no apparent oscillating term. This is in keeping with the results we obtained, where no oscillations were observed. Nevertheless, some other leads remains to be explored and we could add some specificities to our simple model, such as a time dependent albedo, the fraction of ice, salinity, solar radiation penetrating the ice surface... Oscillations might come from one of these.
# 7 Conclusion and Future Work

The goal of this report was to answer two scientific problems. On the one hand, explain the slow convergence of the Schwarz iterative algorithm applied to an ocean-ice-atmosphere coupling problem. On the other hand, solve oscillations issue that often appears as the iterations proceed.

To answer the first question we developed a simplified one dimensional model. For the moment, it allows us to understand why the Schwarz algorithm needs more iterations to converge in ocean regions covered by ice. Indeed, we have shown theoretically that the thicker the ice slab the closer the convergence rate of the algorithm is to one. This theoretical result has been confirmed by numerical experiments.

However we observe none of the oscillating issue we mentioned earlier. From figure 19, the informations we have are that these oscillations are present in sea-ice regions where the ice fraction is rather constant. Furthermore, they appear for high ice temperature close to the melting point. However there is no link between the oscillations and the ice albedo. This is somewhat surprising because when the ice melts at the surface, a thin layer of water appears and a change in the surface albedo comes with it. A next step would thus be to include more physical phenomena in the model or rerun simulations on large scale models without certain physical parts and try to identify which one could be the cause of the oscillations.

<figure><img src="image_24.png" /><figcaption>**Figure 19:** Relations between ice temperature oscillations and ice fraction (left), ice temperature (middle), ice albedo (right).</figcaption></figure>

Another perspective would be to improve the convergence rate in presence of ice. A possible solution is to tune the interface conditions since they are directly linked to the convergence rate. We could also try to find more relevant initial guesses to start the Schwarz iterative process.
## 8 Schwarz Optimization

In order to optimize the convergence of the Schwarz algorithm in the case where we made assumptions on the behaviours of the ice surface temperature and the ice thickness, we can add degrees of freedom to the model at the interface. Of course they must not disturb the physics of the model. At the interface between ice and atmosphere we had

$$ \left( \frac{\partial}{\partial z} - \frac{\alpha_a}{k_a} Id \right) T_a^n(0,t) = - \frac{\alpha_a u_0 - 1}{k_a u_0} T_a^{n-1}(0,t) + \frac{SW_{net}(t) - \epsilon A - k_i \beta T_b}{k_a u_0} $$

which we turn into

$$ \left( \frac{\partial}{\partial z} - \frac{\alpha_a}{k_a} Id \right) T_a^n(0,t) + \lambda T_a^n(0,t) = - \frac{\alpha_a u_0 - 1}{k_a u_0} T_a^{n-1}(0,t) + \frac{SW_{net}(t) - \epsilon A - k_i \beta T_b}{k_a u_0} + \lambda T_a^{n-1}(0,t) $$

In terms of errors, this gives us

$$ \left( \frac{\partial}{\partial z} - \left( \frac{\alpha_a}{k_a} - \lambda \right) Id \right) e_a^n(0,t) = - \left( \frac{\alpha_a u_0 - 1}{k_a u_0} - \lambda \right) e_a^{n-1}(0,t) $$

We easily conclude that the optimal choice for the free parameter is $\lambda = \frac{\alpha_a u_0 - 1}{k_a u_0}$. The interface condition is now given by

$$ \left( \frac{\partial}{\partial z} - \frac{\alpha_a}{k_a u_0} Id \right) T_a^n(0,t) = \begin{cases} \frac{SW_{net}(t) - \epsilon A - k_i \beta T_b}{k_a u_0} & \\ - \frac{\alpha_a}{k_a} \left( \frac{u^n(t) - 1}{u^n(t)} - \frac{u_0 - 1}{u_0} \right) T_a^{n-1}(0,t) + \frac{SW_{net}(t) - \epsilon A - k_i \beta T_b}{k_a u_0} & \\ \frac{\alpha_a}{k_a} \frac{u_0 - 1}{u_0} T_a^{n-1}(0,t) + \frac{SW_{net}(t) - \epsilon A - k_i \beta T_{su}^{n-1}(t)}{k_a} & \end{cases} $$

depending on how many assumptions we make. We notice that we now also need a first guess in atmosphere temperature at the bottom for the most general case.
<figure><img src="image_25.png" /><figcaption>(a) No min in the computation of the ice surface temperature and almost constant ice thickness</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>(b) No min in the computation of the ice surface temperature</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>(c) No assumption</figcaption></figure>

**Figure 20:** Comparison of convergence rates with (right) or without (left) optimisation for different assumptions. Theoretical (dashed lines) and observed (solid lines) rates of convergence.
# A Sturm-Liouville Theory

The Sturm-Liouville theory is the study of real second-order linear ordinary differential equations of the form

$$ \frac{d}{dx} \left( p(x) \frac{du}{dx}(x) \right) + (\lambda r(x) - q(x)) u(x) = 0 \quad (A.1) $$

where $r$ and $q$ are continuous real functions and $p$ is a continuously differentiable real function. The unknowns are the function $u$ and the constant parameter $\lambda$. Along with boundary conditions, this becomes a Sturm-Liouville problem.

$$ \left\{ \begin{array}{ll} \frac{d}{dx} \left( p(x) \frac{du}{dx}(x) \right) + (\lambda r(x) - q(x)) u(x) = 0 & \text{on } [a, b] \\ \mathcal{B}_a u = 0 & \\ \mathcal{B}_b u = 0 & \end{array} \right. \quad (\text{Sturm-Liouville problem}) $$

A Sturm-Liouville problem is said to be regular if the functions $p$ and $r$ are positive in $[a, b]$. We will only consider separated boundary conditions of the form (A.2).

$$ \left\{ \begin{array}{l} \mathcal{B}_a u = \alpha u(a) + \alpha' u'(a) = 0 \\ \mathcal{B}_b u = \beta u(b) + \beta' u'(b) = 0 \end{array} \right. \quad (A.2) $$

with $\alpha, \beta, \alpha', \beta' \in \mathbb{R}$. We further assume that the boundary conditions are nontrivial i.e. $|\alpha| + |\alpha'| > 0$ and $|\beta| + |\beta'| > 0$.

## A.1 Sturm-Liouville Operator

To the Sturm-Liouville problem, we associate the Sturm-Liouville linear operator $\mathcal{S}$ defined in (A.3) on the Hilbert space $\mathcal{H} = L^2([a, b], r(x)dx)$.

$$ \forall u \in \mathcal{H}, \quad \mathcal{S}u = - \frac{(pu')' + qu}{r} \quad (A.3) $$

$\mathcal{H}$ is endowed with the weighted inner product $\langle \cdot, \cdot \rangle_r$ defined as

$$ \forall u, v \in \mathcal{H}, \quad \langle u, v \rangle_r = \int_a^b u(x)^* v(x) r(x) dx $$

So we actually want to solve the eigenvalue problem associated to the linear operator $\mathcal{S}$. This is to say, we want to find couples $(u, \lambda)$ such that $u$ satisfies the boundary conditions (A.2) and the couple verifies

$$ \forall x \in [a, b], \quad (\mathcal{S}u)(x) = \lambda u(x) $$

For such couples, $u$ is called the eigenfunction associated to the eigenvalue $\lambda$.

## Property 1 (Self-adjointness property)
The operator $\mathcal{S}$ is self-adjoint in $\mathcal{H}$ on $D = \{u \in \mathcal{H} / B_a u = B_b u = 0\}$.

Let $u, v \in D$. By definition of $\mathcal{S}$ and the inner product, we have

$$ \begin{aligned} \langle \mathcal{S}u, v \rangle_r - \langle u, \mathcal{S}v \rangle_r &= \int_a^b -(pu')'v + quv + (pv')'u - quv \, dx \\ &= \int_a^b (pv')'u - (pu')'v \, dx \\ &= p(v'u - vu') \Big|_a^b \end{aligned} $$

Since $B_a u = B_a v = 0$, we can show that $p(a)(v'(a)u(a) - v(a)u'(a)) = 0$ and similarly with the boundary condition in $x = b$. So for any $u, v \in D$, $\langle \mathcal{S}u, v \rangle_r = \langle u, \mathcal{S}v \rangle_r$.

## A.2 Main Results

### Theorem 1 (Eigenfunction theorem)

Any regular Sturm-Liouville problem has an infinite sequence of strictly increasing real simple eigenvalues $(\lambda_k)_{k \in \mathbb{N}}$ with $\lim_{k \to +\infty} \lambda_k = +\infty$. The eigenfunction $u_k$ corresponding to the eigenvalue $\lambda_k$ has exactly $k$ zeros in $(a, b)$. It is unique up to a constant factor.

We will prove some of the results stated in this theorem.

* We propose to show that the eigenvalues are real. Assume that $(u, \lambda)$ is a solution to the eigenvalue problem and that $u = u_1 + iu_2$ and $\lambda = \lambda_1 + i\lambda_2$. Let us show that $\lambda_2 = 0$. By definition we have $\mathcal{S}u = \lambda u$, so $\langle u, \mathcal{S}u \rangle_r = \langle u, \lambda u \rangle_r$. If we consider the imaginary part, we get

$$ \int_a^b u_2(pu_1')' - u_1(pu_2')' \, dx = \lambda_2 \int_a^b (u_1^2 + u_2^2)r \, dx $$

The integrand in the left hand side is actually the derivative of $u_2pu_1' - u_1pu_2' = -pW(u_1, u_2)$. So we have

$$ \lambda_2 \int_a^b (u_1^2 + u_2^2)r \, dx = p(a)W(u_1, u_2)(a) - p(b)W(u_1, u_2)(b) $$

However, if $W(u_1, u_2)(a) \neq 0$ then the boundary conditions (A.2) are trivial. With the same reasoning in $x = b$ and since $r > 0$, we deduce that $\lambda_2 = 0$.

* Now we show that each eigenvalue is simple. This means that the eigenfunctions associated to each eigenvalue are linearly dependent. Let $u_1$ and $u_2$ be two eigenfunctions associated to the eigenvalue $\lambda$. Suppose that they are linearly independent. So for any constant $C_1$ and $C_2$, $u = C_1u_1 + C_2u_2$ is a solution of the
same ODE. Furthermore, since the boundary conditions (A.2) are linear, $u$ satisfies them as well. By considering $u$ the solution such that $u(a) = 1$ and $u'(a) = 0$, we have $\alpha = 0$. If we now consider the solution $u$ such that $u(a) = 0$ and $u'(a) = 1$, we get $\alpha' = 0$. This contradicts the assumption of nontrivial boundary conditions.

**Theorem 2 (Orthonormal basis)**

The normalized eigenfunctions $(u_k)_{k \in \mathbb{N}}$ of a regular Sturm-Liouville problem form an orthonormal basis of $\mathcal{H}$.

To show this result, we will use the Self-adjointness property of the operator $\mathcal{S}$. Let $u_k$ and $u_\ell$ be two eigenfunctions associated to two distinct eigenvalues $\lambda_k$ and $\lambda_\ell$ respectively. Then we have that

$$ (\lambda_k - \lambda_\ell) \langle u_k, u_\ell \rangle_r = \langle \lambda_k u_k, u_\ell \rangle_r - \langle \lambda_\ell u_k, u_\ell \rangle_r \\ = \langle \mathcal{S} u_k, u_\ell \rangle_r - \langle u_k, \lambda_\ell u_\ell \rangle_r $$

But since $\mathcal{S}$ is self-adjoint, we deduce that

$$ (\lambda_k - \lambda_\ell) \langle u_k, u_\ell \rangle_r = \langle u_k, \mathcal{S} u_\ell \rangle_r - \langle u_k, \lambda_\ell u_\ell \rangle_r \\ = \langle u_k, \lambda_\ell u_\ell \rangle_r - \langle u_k, \lambda_\ell u_\ell \rangle_r \\ = 0 $$

So $\langle u_k, u_\ell \rangle_r = 0$ for $k \neq \ell$. To obtain an orthonormal basis, we just have to normalize the $u_k$'s to unit.

This orthonormal basis allows us to generalize the Fourier series and write any function $u \in D$ as a series.

$$ u(x) = \sum_{k=0}^{+\infty} \langle u, u_k \rangle_r u_k(x) $$

### A.3 The Case of the 1D Nonstationary Diffusion Equation

Consider the 1D nonstationary homogeneous diffusion equation with homogeneous boundary conditions.

$$ \left\{ \begin{array}{ll} \displaystyle \frac{\partial u}{\partial t} - \nu \frac{\partial^2 u}{\partial z^2} = 0 & \text{on } (0, L) \times (0, T) \\ \mathcal{B}_0 u(0, t) = 0 & \text{on } (0, T) \\ \mathcal{B}_L u(L, t) = 0 & \text{on } (0, T) \\ u(z, 0) = u_0(z) & \text{on } (0, L) \end{array} \right. \qquad (A.4) $$
It can be solved using a separation of variables technique. More precisely, we look for $u(z,t) = Z(z)T(t)$. From the first equation, we deduce that

$$ \forall z \in (0, L), \forall t \in (0, T), \frac{T'(t)}{T(t)} = \nu \frac{Z''(z)}{Z(z)} = c $$

We have managed to turn a PDE into two ODE. The non negative case $c \ge 0$ leads to trivial solutions, so we consider $c = -\lambda < 0$. This gives us for $Z$

$$ \left\{ \begin{array}{l} \nu Z''(z) + \lambda Z(z) = 0 \quad \text{on } (0, L) \\ B_0 Z(0) = 0 \\ B_L Z(L) = 0 \end{array} \right. \qquad (A.5) $$

This is a Sturm-Liouville problem with $p \equiv \nu$, $r \equiv 1$ and $q \equiv 0$. Let us denote $Z_k$ the eigenfunction associated to the eigenvalue $\lambda_k$. We have $T_k(t) = C_k e^{-\lambda_k t}$. So using the superposotion theorem, we deduce that the general solution of (A.4) is given by

$$ u(z,t) = \sum_{k=0}^{+\infty} \alpha_k Z_k(z) T_k(t) $$

where the $\alpha_k$'s are determined using the initial condition $u_0$ and the fact that $(Z_k)_{k \in \mathbb{N}}$ is a basis.

$$ \alpha_k \langle Z_k, Z_k \rangle_r = \langle u_0, Z_k \rangle_r $$

Of course, the eigenvalues and their eigenfunctions depend on the boundary conditions.

* Dirichlet-Dirichlet boundary conditions:

This is the case where $B_0 = B_L = Id$. The eigenfunctions and eigenvalues are given by

$$ Z_k = \sin \left( \frac{k\pi}{L} z \right) \quad \text{and} \quad \lambda_k = \frac{k^2\pi^2\nu}{L^2} $$

* Neumann-Neumann boundary conditions:

This is the case where $B_0 = B_L = \partial_z$. The eigenfunctions and eigenvalues are given by

$$ Z_k = \cos \left( \frac{k\pi}{L} z \right) \quad \text{and} \quad \lambda_k = \frac{k^2\pi^2\nu}{L^2} $$

* Dirichlet-Robin boundary conditions:

This is the case where $B_0 = Id$ and $B_L = \partial_z + r_L Id$. The eigenfunctions and eigenvalues are given by

$$ Z_k = \sin(\mu_k z) \quad \text{and} \quad \lambda_k = \nu \mu_k^2 $$

where $\mu_k$ solves

$$ \mu + r_L \tan(\mu L) = 0 $$
• Robin-Dirichlet boundary conditions:

This is the case where $B_0 = \partial_z + r_0 Id$ and $B_L = Id$. The eigenfunctions and eigenvalues are given by

$$Z_k = \sin(\mu_k z) - \frac{\mu_k}{r_0} \cos(\mu_k z) \quad \text{and} \quad \lambda_k = \nu \mu_k^2$$

where $\mu_k$ solves

$$\mu - r_0 \tan(\mu L) = 0$$

Assume now we want to solve the same equation (A.4) but with a source term $f$.

$$\left\{
\begin{array}{ll}
\displaystyle \frac{\partial u}{\partial t} - \nu \frac{\partial^2 u}{\partial z^2} = f(z,t) & \text{on } (0,L) \times (0,T) \\
\displaystyle B_0 u(0,t) = 0 & \text{on } (0,T) \\
\displaystyle B_L u(L,t) = 0 & \text{on } (0,T) \\
\displaystyle u(z,0) = u_0(z) & \text{on } (0,L)
\end{array}
\right.
\qquad (A.6)$$

We can project $f$ and $u_0$ on the orthonormal basis $(Z_k)_{k \in \mathbb{N}}$ deduced from (A.5).

$$f(z,t) = \sum_{k=0}^{+\infty} f_k(t) Z_k(z) \quad \text{and} \quad u_0(z) = \sum_{k=0}^{+\infty} u_{0,k} Z_k(z)$$

The coefficients $f_k$ and $u_{0,k}$ are given by

$$f_k(t) = \frac{\langle f(.,t), Z_k \rangle_r}{\langle Z_k, Z_k \rangle_r} \quad \text{and} \quad u_{0,k} = \frac{\langle u_0, Z_k \rangle_r}{\langle Z_k, Z_k \rangle_r}$$

Using the separation of variables technique, if we look for $u$ of the form $u(z,t) = \sum_{k=0}^{+\infty} Z_k(z)T_k(t)$, the boundary conditions are verified by construction of the $Z_k$'s and we find

$$\left\{
\begin{array}{ll}
T'_{k}(t) + \lambda_{k}T_{k}(t) = f_{k}(t) & \text{on } (0,T) \\
T_{k}(0) = u_{0,k}
\end{array}
\right.
\qquad (A.7)$$

The solutions of this first order ODE are

$$T_k(t) = u_{0,k}e^{-\lambda_k t} + \int_0^t f_k(s)e^{-\lambda_k(t-s)} ds$$

## B Chini's Equation

Chini's equations are first-order differential equations of the form

$$y'(t) + a_n(t)y^n(t) + a_1(t)y(t) = a_0(t)$$

Here we are interested in the case $n = -1$ and $a_1 \equiv 0$.
## B.1 Existence of a Solution

Let $\Omega = (0, T) \times (1, +\infty)$ and consider the function $f$ defined as follows

$$f: \begin{array}{rcl} \Omega & \to & \mathbb{R} \\ (t, y) & \mapsto & b(t) - \frac{a(t)}{y(t)} \end{array}$$

We want to solve a first-order nonlinear ordinary differential equation of the form

$$\left\{ \begin{array}{ll} y' = f(t, y) & (t, y) \in \Omega \\ y(t_0) = y_0 & (t_0, y_0) \in \Omega \end{array} \right. \qquad (B.1)$$

We have that $\frac{\partial f}{\partial y} = \frac{a}{y^2}$ is continuous on $\Omega$. So $f$ is locally Lipschitz continuous in $y$ and the Cauchy-Lipschitz theorem gives us the existence of a unique solution defined on an interval.

## B.2 Lambert W-function

The Lambert W-function is defined as the inverse function of the strictly increasing function $f$.

$$f: \begin{array}{rcl} [-1, +\infty) & \to & [-e^{-1}, +\infty) \\ w & \mapsto & we^w \end{array}$$

This is to say

$$\forall z \ge -e^{-1}, \quad f(W(z)) = W(z)e^{W(z)} = z \qquad (B.2)$$

This implies $W(z) = 0 \iff z = 0$. Using (B.2), we can numerically compute $W$ by applying a fix point algorithm. The function $W$ is plotted in figure 21. By the inverse

<figure><img src="image_28.png" /><figcaption>Figure 21: Lambert W-function.</figcaption></figure>
function rule, the derivative of W is given by

$$ \forall z > -e^{-1}, W'(z) = \frac{1}{(1 + W(z))e^{W(z)}} \qquad (B.3) $$

From (B.2), we also deduce another relation for the derivative.

$$ \forall z > -e^{-1}, zW'(z) = \frac{W(z)}{1 + W(z)} = 1 - \frac{1}{1 + W(z)} \qquad (B.4) $$

### B.2.1 Simple Case

Let us now consider the function $g: z \mapsto \frac{1 + W(e^{\lambda^2 A(z)})}{\lambda}$, with $\lambda \in \mathbb{R}$ and $A$ a continuously differentiable function. The derivative of $g$ is given by

$$ \forall z \in \mathbb{R}, g'(z) = \lambda a(z) W'(e^{\lambda^2 A(z)}) e^{\lambda^2 A(z)} \qquad (B.5) $$

where $a = A'$. From (B.4), we deduce that

$$ \forall z \in \mathbb{R}, g'(z) = \lambda a(z) \left( 1 - \frac{1}{1 + W(e^{\lambda^2 A(z)})} \right) \qquad (B.6) $$

Hence, $g$ verifies the following differential equation

$$ \forall z \in \mathbb{R}, g'(z) + \frac{a(z)}{g(z)} = \lambda a(z) \qquad (B.7) $$

### B.2.2 Autonomous Case

Now consider the function $g: z \mapsto \frac{a}{b} \left( 1 + W \left( e^{\frac{b^2}{a} z} \right) \right)$, with $a, b \in \mathbb{R}$. This time, the derivative of $g$ is given by

$$ \forall z \in \mathbb{R}, g'(z) = bW' \left( e^{\frac{b^2}{a} z} \right) e^{\frac{b^2}{a} z} \qquad (B.8) $$

Still from (B.4), we deduce that

$$ \forall z \in \mathbb{R}, g'(z) = b \left( 1 - \frac{1}{1 + W \left( e^{\frac{b^2}{a} z} \right)} \right) \qquad (B.9) $$

Hence, $g$ verifies the following differential equation

$$ \forall z \in \mathbb{R}, g'(z) + \frac{a}{g(z)} = b \qquad (B.10) $$

So $g$ is solution of (B.1) in the autonomous case.
C Shifting Method

In this section, we will see how to turn a problem with inhomogeneous boundary conditions into an equivalent problem with homogeneous boundary conditions. Assume we want to solve the following PDE problem where $\mathcal{L}$, $\mathcal{B}_0$ and $\mathcal{B}_L$ are linear partial derivative operators.

$$ \left\{ \begin{array}{ll} \mathcal{L}u(z,t) = f(z,t) & \text{on } (0,L) \times (0,T) \\ \mathcal{B}_0u(0,t) = g_0(t) & \text{on } (0,T) \\ \mathcal{B}_Lu(L,t) = g_L(t) & \text{on } (0,T) \\ u(z,0) = u_0(z) & \text{on } (0,L) \end{array} \right. \qquad (C.1) $$

The shifting method consists of finding a solution $u$ of the form $u(z,t) = v(z,t) + w(z,t)$, with $w$ that verifies the boundary conditions i.e. $\mathcal{B}_0w(0,t) = g_0(t)$ and $\mathcal{B}_Lw(L,t) = g_L(t)$. Injecting such a $u$ in (C.1) implies that $v$ must be solution of the following homogeneous boundary conditions problem.

$$ \left\{ \begin{array}{ll} \mathcal{L}v(z,t) = f(z,t) - \mathcal{L}w(z,t) & \text{on } (0,L) \times (0,T) \\ \mathcal{B}_0v(0,t) = 0 & \text{on } (0,T) \\ \mathcal{B}_Lv(L,t) = 0 & \text{on } (0,T) \\ v(z,0) = u_0(z) - w(z,0) & \text{on } (0,L) \end{array} \right. \qquad (C.2) $$

A typical choice for $w$ is $w(z,t) = a(t)z + b(t)$. $a$ and $b$ are to be determined using the boundary conditions.
References

[1] O. Marti, S. Nguyen, P. Braconnot, S. Valcke, F. Lemarié, and E. Blayo. A schwarz iterative method to evaluate ocean–atmosphere coupling schemes: implementation and diagnostics in ipsl-cm6-sw-vlr. *Geoscientific Model Development*, 14(5):2959–2975, 2021.

[2] Gary A. Maykut and Norbert Untersteiner. Some results from a time-dependent thermodynamic model of sea ice. *Journal of Geophysical Research (1896-1977)*, 76(6):1550–1575, 1971.

[3] Albert J. Semtner. A model for the thermodynamic growth of sea ice in numerical investigations of climate. *Journal of Physical Oceanography*, 6(3):379–389, 1976.

[4] A. S. Thorndike. A toy model linking atmospheric thermal radiation and sea ice growth. *Journal of Geophysical Research: Oceans*, 97(C6):9401–9410, 1992.