# A second-order correction method for loosely coupled discretizations applied to parabolic-parabolic interface problems

Erik Burman, Rebecca Durst, Miguel Angel Fernández, Johnny Guzmán, Sijing Liu

► To cite this version:

Erik Burman, Rebecca Durst, Miguel Angel Fernández, Johnny Guzmán, Sijing Liu. A second-order correction method for loosely coupled discretizations applied to parabolic-parabolic interface problems. 2024. [https://hal.science/hal-04530926](https://hal.science/hal-04530926)

HAL Id: hal-04530926

https://hal.science/hal-04530926

Preprint submitted on 3 Apr 2024

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.
# A second-order correction method for loosely coupled discretizations applied to parabolic-parabolic interface problems

Erik Burman, Rebecca Durst, Miguel Angel Fernández, Johnny Guzmán, Sijing Liu

► To cite this version:

Erik Burman, Rebecca Durst, Miguel Angel Fernández, Johnny Guzmán, Sijing Liu. A second-order correction method for loosely coupled discretizations applied to parabolic-parabolic interface problems. IMA Journal of Numerical Analysis, In press. hal-04530926

HAL Id: hal-04530926

https://hal.science/hal-04530926v1

Submitted on 3 Apr 2024

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.

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

Distributed under a Creative Commons Attribution 4.0 International License
A SECOND-ORDER CORRECTION METHOD FOR LOOSELY COUPLED
DISCRETIZATIONS APPLIED TO PARABOLIC-PARABOLIC INTERFACE
PROBLEMS

ERIK BURMAN, REBECCA DURST, MIGUEL A. FERNÁNDEZ, JOHNNY GUZMÁN, AND SIJING LIU

ABSTRACT. We consider a parabolic-parabolic interface problem and construct a loosely coupled prediction-correction scheme based on the Robin-Robin splitting method analyzed in [*J. Numer. Math.*, **31**(1):59–77, 2023]. We show that the errors of the correction step converge at $O((\Delta t)^2)$, under suitable convergence rate assumptions on the discrete time derivative of the prediction step, where $\Delta t$ stands for the time-step length. Numerical results are shown to support our analysis and the assumptions.

# 1. INTRODUCTION

A number of studies have been reported in the literature on loosely coupled schemes for incompressible fluid-structure interaction (FSI) with thick-walled solids (see, e.g., [19, 5, 20, 12, 26, 30, 13, 27, 28, 10]). In most of these works, either no rigorous error analysis is carried out or error estimates are provided which are not uniform with respect to the spatial discretization. Robin-Robin loosely coupled schemes, first proposed in [20] for FSI problems, were recently analyzed in [11, 18, 16, 14, 17, 25]. Both the parabolic-parabolic and the hyperbolic-parabolic couplings have been considered in [14]. The FSI case corresponds to the hyperbolic-parabolic coupling with the additional difficulty coming from the fluid incompressibility constraint [11, 16, 25, 17]. The method is shown to be unconditionally stable and, for the first time, with a splitting error that is spatially uniform. Sub-optimal convergence rates in time of order $O(\sqrt{\Delta t})$ were first proved in [11, 18, 16]. However, this was later improved to nearly first-order accuracy in [25, 14, 17].

By leveraging the good stability properties and spatially uniform splitting error of the method in [16, 17], our aim is to construct a higher-order scheme in time by building on the prediction-correction methodology [9], with the Robin-Robin loosely coupled method as the prediction step. Although our ultimate goal is to develop a higher-order method for FSI problems, in this paper we start by constructing and analyzing a method for a coupled parabolic-parabolic problem. Indeed, the parabolic-parabolic problem is a classical model problem for fluid-fluid interaction, but it also involves some of the difficulties present in linear FSI problems, making it a good first step. It should be mentioned that correction methods have been used widely for FSI problems; see for example, [3, 19, 5, 20, 16, 17]. Burman and Fernández [19] already argued that one should expect a $\Delta t$ improvement for every correction step. However, there does not seem to be an analysis of this for FSI problems (or for parabolic-parabolic problems). A method based on subcycling (several correction steps) was proposed in [7, 6], with linear convergence in the number of steps. In a domain decomposition framework, focusing on a multi-timestep approach, a Robin-Robin coupling for time-dependent advection–diffusion was introduced in [21] with numerical investigation of the stability.

In this paper, we aim to provide rigorous evidence of this by proving second-order convergence of a prediction-correction method with a single correction step, under an assumption of the convergence of the prediction step. The guiding principle of our correction method is to have discrete time differences of the prediction step error appearing in the right-hand side of the correction step error equations. Then, one uses the fact that the discrete time differences converge with second-order accuracy. We should mention that this strategy has also been used in [1] for a non-linear fluid problem, and the idea goes back to the work of Stetter and Böhmer et al. [32, 8, 9]. Typically, in the previous work on
A SECOND-ORDER CORRECTION METHOD FOR LOOSELY COUPLED
DISCRETIZATIONS APPLIED TO PARABOLIC-PARABOLIC INTERFACE
PROBLEMS

ERIK BURMAN, REBECCA DURST, MIGUEL A. FERNÁNDEZ, JOHNNY GUZMÁN, AND SIJING LIU

ABSTRACT. We consider a parabolic-parabolic interface problem and construct a loosely coupled prediction-correction scheme based on the Robin-Robin splitting method analyzed in [*J. Numer. Math.*, **31**(1):59–77, 2023]. We show that the errors of the correction step converge at $O((\Delta t)^2)$, under suitable convergence rate assumptions on the discrete time derivative of the prediction step, where $\Delta t$ stands for the time-step length. Numerical results are shown to support our analysis and the assumptions.

# 1. INTRODUCTION

A number of studies have been reported in the literature on loosely coupled schemes for incompressible fluid-structure interaction (FSI) with thick-walled solids (see, e.g., [19, 5, 20, 12, 26, 30, 13, 27, 28, 10]). In most of these works, either no rigorous error analysis is carried out or error estimates are provided which are not uniform with respect to the spatial discretization. Robin-Robin loosely coupled schemes, first proposed in [20] for FSI problems, were recently analyzed in [11, 18, 16, 14, 17, 25]. Both the parabolic-parabolic and the hyperbolic-parabolic couplings have been considered in [14]. The FSI case corresponds to the hyperbolic-parabolic coupling with the additional difficulty coming from the fluid incompressibility constraint [11, 16, 25, 17]. The method is shown to be unconditionally stable and, for the first time, with a splitting error that is spatially uniform. Sub-optimal convergence rates in time of order $O(\sqrt{\Delta t})$ were first proved in [11, 18, 16]. However, this was later improved to nearly first-order accuracy in [25, 14, 17].

By leveraging the good stability properties and spatially uniform splitting error of the method in [16, 17], our aim is to construct a higher-order scheme in time by building on the prediction-correction methodology [9], with the Robin-Robin loosely coupled method as the prediction step. Although our ultimate goal is to develop a higher-order method for FSI problems, in this paper we start by constructing and analyzing a method for a coupled parabolic-parabolic problem. Indeed, the parabolic-parabolic problem is a classical model problem for fluid-fluid interaction, but it also involves some of the difficulties present in linear FSI problems, making it a good first step. It should be mentioned that correction methods have been used widely for FSI problems; see for example, [3, 19, 5, 20, 16, 17]. Burman and Fernández [19] already argued that one should expect a $\Delta t$ improvement for every correction step. However, there does not seem to be an analysis of this for FSI problems (or for parabolic-parabolic problems). A method based on subcycling (several correction steps) was proposed in [7, 6], with linear convergence in the number of steps. In a domain decomposition framework, focusing on a multi-timestep approach, a Robin-Robin coupling for time-dependent advection–diffusion was introduced in [21] with numerical investigation of the stability.

In this paper, we aim to provide rigorous evidence of this by proving second-order convergence of a prediction-correction method with a single correction step, under an assumption of the convergence of the prediction step. The guiding principle of our correction method is to have discrete time differences of the prediction step error appearing in the right-hand side of the correction step error equations. Then, one uses the fact that the discrete time differences converge with second-order accuracy. We should mention that this strategy has also been used in [1] for a non-linear fluid problem, and the idea goes back to the work of Stetter and Böhmer et al. [32, 8, 9]. Typically, in the previous work on
defect-correction for interface problems [1], the physical coupling conditions in fluid-fluid interaction introduced some dissipation, which made it possible to control the splitting error. This is often the case for fluid-fluid interactions modelling atmosphere-ocean interaction, where the two systems are coupled through the drag force on the ocean surface [23, 24, 22, 34, 33, 35, 31, 29]. In our case, no such dissipative mechanism is present, and therefore the arguments of [1] fail. The conservative coupling presents an important challenge for the analysis, both for the low-order splitting scheme and the defect-correction scheme. Indeed, to the best of our knowledge, this is the first time second-order convergence is proved for a splitting scheme where the coupling conditions are conservative. A key milestone towards achieving second-order convergence of the correction step was reached in an accompanying paper [15], where we prove that many of the time difference errors of the prediction step are indeed second-order accurate in the $L^2$-norm. When the interface is perpendicular to the two sides of the domain, we also prove that the $H^2$-norm of the time difference errors of the prediction step also converges as $\mathcal{O}((\Delta t)^2)$. This $H^2$ estimate turns out to be a crucial ingredient in the analysis of the correction method. Although it is nontrivial to prove the $H^2$ estimate in a more general setting, numerical results indicate that the $H^2$ estimate also holds for more general cases.

Let us summarize the steps of the prediction-correction method at time $t_{n+1}$:

* Solve the prediction step to obtain the defect solutions $u_0^{n+1}$, $w_0^{n+1}$ and $\lambda_0^{n+1}$.

* Modify the right-hand side of the prediction step using the solutions $u_0^{n+1}$, $w_0^{n+1}$ and $\lambda_0^{n+1}$, while keeping the left-hand side of the prediction step unchanged.

* The modification of the right-hand sides is designed in such a way that the right-hand sides of the error equations of the correction step consists of the term $U_0^{n+1} - U_0^n$, $W_0^{n+1} - W_0^n$, $\Lambda_0^{n+1} - \Lambda_0^n$ and their gradients.

* Solve the new system with the modified right-hand side to obtain the correction solutions $u_1^{n+1}$, $w_1^{n+1}$ and $\lambda_1^{n+1}$.

Here, $\lambda_0$ denotes the Lagrange multiplier for the prediction step. The terms $U_0^n$, $W_0^n$ and $\Lambda_0^n$ are the corresponding error terms of the variables $u_0^n$, $w_0^n$ and $\lambda_0^n$ at time $t_n$. As mentioned above, it has been observed numerically that $U_0^{n+1} - U_0^n$, $W_0^{n+1} - W_0^n$, and $\Lambda_0^{n+1} - \Lambda_0^n$ converge with rate $\mathcal{O}((\Delta t)^2)$ in the $L^2$-norm in general. Moreover, we prove that the first two errors are second-order accurate in the accompanying paper [15]. We also prove that the $H^2$-norm of $U_0^{n+1} - U_0^n$ is second-order accurate for a special case. Note that this implies the second-order convergence of the Lagrange multiplier in $L^2$-norm. In Assumption 3.1 below, we assume the second-order accuracy of this error in order to carry out the analysis of the correction step. With this assumption we are able to prove second-order accuracy of the correction step. It should be noted that we perform the analysis for the time semi-discrete method. Extension to the fully discrete case are non-trivial and will be considered in the future.

The rest of the paper is organized as follows. In Section 2, we introduce a parabolic-parabolic interface problem and the corresponding prediction-correction method. In Section 4, we show that the correction method gives higher-order convergence in time under a crucial lemma when the interface is perpendicular to the two sides of $\Omega$. In Section 5, we present two numerical examples that illustrate the theoretical findings. In Section 6, we consider the interface problem with Dirichlet boundary conditions and discuss an inconsistency that arises and briefly describe how to remedy it. We conclude the paper by drawing some future directions in Section 7.

## 2. PREDICTION-CORRECTION METHOD FOR A PARABOLIC\PARABOLIC INTERFACE PROBLEM

Let $\Omega$ be a bounded domain that can be decomposed as $\Omega = \Omega_f \cup \Omega_s$ where the common interface $\Sigma = \partial\Omega_f \cap \partial\Omega_s$. The subset $\Gamma_D^f$ (resp., $\Gamma_D^s$) of the boundary $\partial\Omega_f$ (resp., $\partial\Omega_s$) is the bottom side (resp., top side) of the domain $\Omega_f$ (resp., $\Omega_s$). We then denote $\Gamma_N^f = \partial\Omega_f \setminus (\Gamma_D^f \cap \Sigma)$ and $\Gamma_N^s = \partial\Omega_s \setminus (\Gamma_D^s \cap \Sigma)$. See Figure 1 for an illustration.
defect-correction for interface problems [1], the physical coupling conditions in fluid-fluid interaction introduced some dissipation, which made it possible to control the splitting error. This is often the case for fluid-fluid interactions modelling atmosphere-ocean interaction, where the two systems are coupled through the drag force on the ocean surface [23, 24, 22, 34, 33, 35, 31, 29]. In our case, no such dissipative mechanism is present, and therefore the arguments of [1] fail. The conservative coupling presents an important challenge for the analysis, both for the low-order splitting scheme and the defect-correction scheme. Indeed, to the best of our knowledge, this is the first time second-order convergence is proved for a splitting scheme where the coupling conditions are conservative. A key milestone towards achieving second-order convergence of the correction step was reached in an accompanying paper [15], where we prove that many of the time difference errors of the prediction step are indeed second-order accurate in the $L^2$-norm. When the interface is perpendicular to the two sides of the domain, we also prove that the $H^2$-norm of the time difference errors of the prediction step also converges as $\mathcal{O}((\Delta t)^2)$. This $H^2$ estimate turns out to be a crucial ingredient in the analysis of the correction method. Although it is nontrivial to prove the $H^2$ estimate in a more general setting, numerical results indicate that the $H^2$ estimate also holds for more general cases.

Let us summarize the steps of the prediction-correction method at time $t_{n+1}$:

* Solve the prediction step to obtain the defect solutions $u_0^{n+1}$, $w_0^{n+1}$ and $\lambda_0^{n+1}$.

* Modify the right-hand side of the prediction step using the solutions $u_0^{n+1}$, $w_0^{n+1}$ and $\lambda_0^{n+1}$, while keeping the left-hand side of the prediction step unchanged.

* The modification of the right-hand sides is designed in such a way that the right-hand sides of the error equations of the correction step consists of the term $U_0^{n+1} - U_0^n$, $W_0^{n+1} - W_0^n$, $\Lambda_0^{n+1} - \Lambda_0^n$ and their gradients.

* Solve the new system with the modified right-hand side to obtain the correction solutions $u_1^{n+1}$, $w_1^{n+1}$ and $\lambda_1^{n+1}$.

Here, $\lambda_0$ denotes the Lagrange multiplier for the prediction step. The terms $U_0^n$, $W_0^n$ and $\Lambda_0^n$ are the corresponding error terms of the variables $u_0^n$, $w_0^n$ and $\lambda_0^n$ at time $t_n$. As mentioned above, it has been observed numerically that $U_0^{n+1} - U_0^n$, $W_0^{n+1} - W_0^n$, and $\Lambda_0^{n+1} - \Lambda_0^n$ converge with rate $\mathcal{O}((\Delta t)^2)$ in the $L^2$-norm in general. Moreover, we prove that the first two errors are second-order accurate in the accompanying paper [15]. We also prove that the $H^2$-norm of $U_0^{n+1} - U_0^n$ is second-order accurate for a special case. Note that this implies the second-order convergence of the Lagrange multiplier in $L^2$-norm. In Assumption 3.1 below, we assume the second-order accuracy of this error in order to carry out the analysis of the correction step. With this assumption we are able to prove second-order accuracy of the correction step. It should be noted that we perform the analysis for the time semi-discrete method. Extension to the fully discrete case are non-trivial and will be considered in the future.

The rest of the paper is organized as follows. In Section 2, we introduce a parabolic-parabolic interface problem and the corresponding prediction-correction method. In Section 4, we show that the correction method gives higher-order convergence in time under a crucial lemma when the interface is perpendicular to the two sides of $\Omega$. In Section 5, we present two numerical examples that illustrate the theoretical findings. In Section 6, we consider the interface problem with Dirichlet boundary conditions and discuss an inconsistency that arises and briefly describe how to remedy it. We conclude the paper by drawing some future directions in Section 7.

## 2. PREDICTION-CORRECTION METHOD FOR A PARABOLIC\PARABOLIC INTERFACE PROBLEM

Let $\Omega$ be a bounded domain that can be decomposed as $\Omega = \Omega_f \cup \Omega_s$ where the common interface $\Sigma = \partial\Omega_f \cap \partial\Omega_s$. The subset $\Gamma_D^f$ (resp., $\Gamma_D^s$) of the boundary $\partial\Omega_f$ (resp., $\partial\Omega_s$) is the bottom side (resp., top side) of the domain $\Omega_f$ (resp., $\Omega_s$). We then denote $\Gamma_N^f = \partial\Omega_f \setminus (\Gamma_D^f \cap \Sigma)$ and $\Gamma_N^s = \partial\Omega_s \setminus (\Gamma_D^s \cap \Sigma)$. See Figure 1 for an illustration.
We consider the following parabolic\parabolic interface problem.

$$
\begin{align}
(2.1a) \qquad & \partial_t u - \nu_f \Delta u = 0, && \text{in } [0, T] \times \Omega_f, \notag \\
& u(0, x) = u_0(x), && \text{on } \Omega_f, \notag \\
& u = 0, && \text{on } [0, T] \times \Gamma_D^f, \notag \\
& \frac{\partial u}{\partial n} = 0, && \text{on } [0, T] \times \Gamma_N^f \notag
\end{align}
$$

$$
\begin{align}
(2.1b) \qquad & \partial_t w - \nu_s \Delta w = 0, && \text{in } [0, T] \times \Omega_s, \notag \\
& w(0, x) = w_0(x), && \text{on } \Omega_s, \notag \\
& w = 0, && \text{on } [0, T] \times \Gamma_D^s, \notag \\
& \frac{\partial w}{\partial n} = 0, && \text{on } [0, T] \times \Gamma_N^s \notag
\end{align}
$$

$$
(2.1c) \qquad w - u = 0, \qquad \text{in } [0, T] \times \Sigma,
$$

$$
(2.1d) \qquad \nu_s \nabla w \cdot n_s + \nu_f \nabla u \cdot n_f = 0, \qquad \text{in } [0, T] \times \Sigma,
$$

where $\mathbf{n} = \mathbf{n}_f$ and $\mathbf{n}_s$ are the outward facing normal vectors of $\Sigma$ for $\Omega_f$ and $\Omega_s$, respectively. Here $\nu_f$ and $\nu_s$ are piece-wise constant functions. We also assume that the initial data is smooth and that $u$ and $w$ is smooth on $\Omega_f$ and $\Omega_s$ respectively.

Note that we consider Neumann boundary conditions on two sides of the domain instead of Dirichlet boundary conditions, this is different from our previous work [14]. There are two main reasons we consider the problem (2.1). First, the parabolic\parabolic interface problem (2.1) is more applicable in real-world scenarios, as numerous fluid-structure interaction problems involve a specific section within a broader domain where the fluid flows from left to right. Consequently, imposing Dirichlet boundary conditions to two vertical sides of the domain is unsuitable. Secondly, the numerical method for the Dirichlet problem has some inconsistencies which presents some difficulties in the analysis. We discuss this in Section 6.

<figure><img src="image_1.png" /><figcaption>FIGURE 1. The domains $\Omega_f$ and $\Omega_s$ with interface $\Sigma$ and Neumann boundaries.</figcaption></figure>
We consider the following parabolic\parabolic interface problem.

$$
\begin{align}
(2.1a) \qquad & \partial_t u - \nu_f \Delta u = 0, && \text{in } [0, T] \times \Omega_f, \notag \\
& u(0, x) = u_0(x), && \text{on } \Omega_f, \notag \\
& u = 0, && \text{on } [0, T] \times \Gamma_D^f, \notag \\
& \frac{\partial u}{\partial n} = 0, && \text{on } [0, T] \times \Gamma_N^f \notag
\end{align}
$$

$$
\begin{align}
(2.1b) \qquad & \partial_t w - \nu_s \Delta w = 0, && \text{in } [0, T] \times \Omega_s, \notag \\
& w(0, x) = w_0(x), && \text{on } \Omega_s, \notag \\
& w = 0, && \text{on } [0, T] \times \Gamma_D^s, \notag \\
& \frac{\partial w}{\partial n} = 0, && \text{on } [0, T] \times \Gamma_N^s \notag
\end{align}
$$

$$
(2.1c) \qquad w - u = 0, \qquad \text{in } [0, T] \times \Sigma,
$$

$$
(2.1d) \qquad \nu_s \nabla w \cdot n_s + \nu_f \nabla u \cdot n_f = 0, \qquad \text{in } [0, T] \times \Sigma,
$$

where $\mathbf{n} = \mathbf{n}_f$ and $\mathbf{n}_s$ are the outward facing normal vectors of $\Sigma$ for $\Omega_f$ and $\Omega_s$, respectively. Here $\nu_f$ and $\nu_s$ are piece-wise constant functions. We also assume that the initial data is smooth and that $u$ and $w$ is smooth on $\Omega_f$ and $\Omega_s$ respectively.

Note that we consider Neumann boundary conditions on two sides of the domain instead of Dirichlet boundary conditions, this is different from our previous work [14]. There are two main reasons we consider the problem (2.1). First, the parabolic\parabolic interface problem (2.1) is more applicable in real-world scenarios, as numerous fluid-structure interaction problems involve a specific section within a broader domain where the fluid flows from left to right. Consequently, imposing Dirichlet boundary conditions to two vertical sides of the domain is unsuitable. Secondly, the numerical method for the Dirichlet problem has some inconsistencies which presents some difficulties in the analysis. We discuss this in Section 6.

<figure><img src="image_2.png" /><figcaption>FIGURE 1. The domains $\Omega_f$ and $\Omega_s$ with interface $\Sigma$ and Neumann boundaries.</figcaption></figure>
2.1. **Variational form.** Let $(\cdot, \cdot)_i$ be the $L^2$-inner product on $\Omega_i$ for $i = f, s$. Moreover, let $\langle \cdot, \cdot \rangle$ be the $L^2$-inner product on $\Sigma$. Consider the spaces

$$ (2.2a) \qquad V_f := \{v \in H^1(\Omega_f) : v = 0 \text{ on } \Gamma_D^f\}, $$

$$ (2.2b) \qquad V_s := \{v \in H^1(\Omega_s) : v = 0 \text{ on } \Gamma_D^s\}, $$

$$ (2.2c) \qquad V_g := L^2(\Sigma). $$

We define the discrete time derivative as

$$ \partial_{\Delta t} v^{n+1} := \frac{v^{n+1} - v^n}{\Delta t}, $$

and the discrete average in time

$$ v^{n+\frac{1}{2}} := \frac{v^{n+1} + v^n}{2}. $$

Furthermore, define the time step $\Delta t$ and an integer $N$, such that $T = N\Delta t$, and let $u^n = u(t_n, \cdot)$, where $t_n = n\Delta t$ for $n \in \{0, 1, 2, \dots, N\}$. Define $l = \nu_f \nabla u \cdot n_f$ and assuming that $l \in L^2(\Sigma)$ we get that the solution to (2.1) solves, for all $t \in [0, T]$:

$$ (2.3) \qquad (\partial_t u, v)_f + \nu_f (\nabla u, \nabla v)_f - \langle l, v \rangle = 0, \quad \forall v \in V_f, $$

$$ (2.4) \qquad (\partial_t w, q)_s + \nu_s (\nabla w, \nabla q)_s + \langle l, q \rangle = 0, \quad \forall q \in V_s, $$

$$ (2.5) \qquad \langle u - w, \mu \rangle = 0, \quad \forall \mu \in V_g. $$

2.2. **Defect-correction method.** We propose the following defect-correction method. The method consists of two steps, namely, the prediction step and the correction step. The prediction step is a Robin-Robin coupling method [14]. The correction step is designed by modifying the right-hand side of the prediction step while keeping the left-hand side of the prediction step unchanged as we discussed in the introduction.

**Prediction step.** Find $w_0^{n+1} \in V_s$, $u_0^{n+1} \in V_f$, and $\lambda_0^{n+1} \in V_g$ such that for $n \ge 0$,

$$ (2.6a) \qquad (\partial_{\Delta t} w_0^{n+1}, z)_s + \nu_s (\nabla w_0^{n+1}, \nabla z)_s + \alpha \langle w_0^{n+1} - u_0^n, z \rangle + \langle \lambda_0^n, z \rangle = 0, \quad \forall z \in V_s, $$

$$ (2.6b) \qquad (\partial_{\Delta t} u_0^{n+1}, v)_f + \nu_f (\nabla u_0^{n+1}, \nabla v)_f - \langle \lambda_0^{n+1}, v \rangle = 0, \quad \forall v \in V_f, $$

$$ (2.6c) \qquad \langle (\lambda_0^{n+1} - \lambda_0^n) + \alpha(u_0^{n+1} - w_0^{n+1}), \mu \rangle = 0, \quad \forall \mu \in V_g. $$

**Correction step.** Find $w_1^{n+1} \in V_s$, $u_1^{n+1} \in V_f$, and $\lambda_1^{n+1} \in V_g$ such that for $n \ge 0$,

$$ (2.7a) \qquad (\partial_{\Delta t} w_1^{n+1}, z)_s + \nu_s (\nabla w_1^{n+1}, \nabla z)_s + \alpha \langle w_1^{n+1} - u_1^n, z \rangle + \langle \lambda_1^n, z \rangle \\ = \nu_s (\nabla w_0^{n+1}, \nabla z)_s + \alpha \langle w_0^{n+1} - w_0^n, z \rangle + \langle \lambda_0^n, z \rangle - \nu_s (\nabla w_0^{n+\frac{1}{2}}, \nabla z)_s - \langle \lambda_0^{n+\frac{1}{2}}, z \rangle, \quad \forall z \in V_s, $$

$$ (2.7b) \qquad (\partial_{\Delta t} u_1^{n+1}, v)_f + \nu_f (\nabla u_1^{n+1}, \nabla v)_f - \langle \lambda_1^{n+1}, v \rangle \\ = \nu_f (\nabla u_0^{n+1}, \nabla v)_f - \langle \lambda_0^{n+1}, v \rangle - \nu_f (\nabla u_0^{n+\frac{1}{2}}, \nabla v)_f + \langle \lambda_0^{n+\frac{1}{2}}, v \rangle, \quad \forall v \in V_f, $$

$$ (2.7c) \qquad \langle (\lambda_1^{n+1} - \lambda_1^n) + \alpha(u_1^{n+1} - w_1^{n+1}), \mu \rangle = \langle \lambda_0^{n+1} - \lambda_0^n, \mu \rangle, \quad \forall \mu \in V_g. $$

In practice, the splitting method should be implemented sequentially. Assuming the information from the previous time-step is known, we solve the prediction step (2.6) for $w_0^{n+1}$, $u_0^{n+1}$ and $\lambda_0^{n+1}$. We then compute the right-hand side of the correction step (2.7) using $w_0^{n+1}$, $u_0^{n+1}$ and $\lambda_0^{n+1}$ and solve for $w_1^{n+1}$, $u_1^{n+1}$ and $\lambda_1^{n+1}$ using (2.7). Moreover, by design, each problem (2.6), (2.7) can be split into two sub-problems.
2.1. **Variational form.** Let $(\cdot, \cdot)_i$ be the $L^2$-inner product on $\Omega_i$ for $i = f, s$. Moreover, let $\langle \cdot, \cdot \rangle$ be the $L^2$-inner product on $\Sigma$. Consider the spaces

$$
(2.2a) \qquad V_f := \{v \in H^1(\Omega_f) : v = 0 \text{ on } \Gamma_D^f\},
$$

$$
(2.2b) \qquad V_s := \{v \in H^1(\Omega_s) : v = 0 \text{ on } \Gamma_D^s\},
$$

$$
(2.2c) \qquad V_g := L^2(\Sigma).
$$

We define the discrete time derivative as

$$
\partial_{\Delta t} v^{n+1} := \frac{v^{n+1} - v^n}{\Delta t},
$$

and the discrete average in time

$$
v^{n+\frac{1}{2}} := \frac{v^{n+1} + v^n}{2}.
$$

Furthermore, define the time step $\Delta t$ and an integer $N$, such that $T = N\Delta t$, and let $u^n = u(t_n, \cdot)$, where $t_n = n\Delta t$ for $n \in \{0, 1, 2, \dots, N\}$. Define $l = \nu_f \nabla u \cdot n_f$ and assuming that $l \in L^2(\Sigma)$ we get that the solution to (2.1) solves, for all $t \in [0, T]$:

$$
(2.3) \qquad (\partial_t u, v)_f + \nu_f (\nabla u, \nabla v)_f - \langle l, v \rangle =0, \quad \forall v \in V_f,
$$

$$
(2.4) \qquad (\partial_t w, q)_s + \nu_s (\nabla w, \nabla q)_s + \langle l, q \rangle = 0, \quad \forall q \in V_s,
$$

$$
(2.5) \qquad \langle \mathbf{u} - \mathbf{w}, \mu \rangle = 0, \quad \forall \mu \in V_g.
$$

2.2. Defect-correction method. We propose the following defect-correction method. The method consists of two steps, namely, the prediction step and the correction step. The prediction step is a Robin-Robin coupling method [14]. The correction step is designed by modifying the right-hand side of the prediction step while keeping the left-hand side of the prediction step unchanged as we discussed in the introduction.

**Prediction step.** Find $w_0^{n+1} \in V_s$, $u_0^{n+1} \in V_f$, and $\lambda_0^{n+1} \in V_g$ such that for $n \ge 0$,

$$
(2.6a) \quad (\partial_{\Delta t} w_0^{n+1}, z)_s + \nu_s (\nabla w_0^{n+1}, \nabla z)_s + \alpha \langle w_0^{n+1} - u_0^n, z \rangle + \langle \lambda_0^n, z \rangle = 0, \quad \forall z \in V_s,
$$

$$
(2.6b) \qquad (\partial_{\Delta t} u_0^{n+1}, v)_f + \nu_f (\nabla u_0^{n+1}, \nabla v)_f - \langle \lambda_0^{n+1}, v \rangle = 0, \quad \forall v \in V_f,
$$

$$
(2.6c) \qquad \langle (\lambda_0^{n+1} - \lambda_0^n) + \alpha(u_0^{n+1} - w_0^{n+1}), \mu \rangle = 0, \quad \forall \mu \in V_g.
$$

**Correction step.** Find $w_1^{n+1} \in V_s$, $u_1^{n+1} \in V_f$, and $\lambda_1^{n+1} \in V_g$ such that for $n \ge 0$,

$$
(2.7a) \quad (\partial_{\Delta t} w_1^{n+1}, z)_s + \nu_s (\nabla w_1^{n+1}, \nabla z)_s + \alpha \langle w_1^{n+1} - u_1^n, z \rangle + \langle \lambda_1^n, z \rangle \\
\qquad = \nu_s (\nabla w_0^{n+1}, \nabla z)_s + \alpha \langle w_0^{n+1} - w_0^n, z \rangle + \langle \lambda_0^n, z \rangle - \nu_s (\nabla w_0^{n+\frac{1}{2}}, \nabla z)_s - \langle \lambda_0^{n+\frac{1}{2}}, z \rangle, \quad \forall z \in V_s,
$$

$$
(2.7b) \quad (\partial_{\Delta t} u_1^{n+1}, v)_f + \nu_f (\nabla u_1^{n+1}, \nabla v)_f - \langle \lambda_1^{n+1}, v \rangle \\
\qquad = \nu_f (\nabla u_0^{n+1}, \nabla v)_f - \langle \lambda_0^{n+1}, v \rangle - \nu_f (\nabla u_0^{n+\frac{1}{2}}, \nabla v)_f + \langle \lambda_0^{n+\frac{1}{2}}, v \rangle, \quad \forall v \in V_f,
$$

$$
(2.7c) \quad \langle (\lambda_1^{n+1} - \lambda_1^n) + \alpha(u_1^{n+1} - w_1^{n+1}), \mu \rangle = \langle \lambda_0^{n+1} - \lambda_0^n, \mu \rangle, \quad \forall \mu \in V_g.
$$

In practice, the splitting method should be implemented sequentially. Assuming the information from
the previous time-step is known, we solve the prediction step (2.6) for $w_0^{n+1}$, $u_0^{n+1}$ and $\lambda_0^{n+1}$. We
then compute the right-hand side of the correction step (2.7) using $w_0^{n+1}$, $u_0^{n+1}$ and $\lambda_0^{n+1}$ and solve
for $w_1^{n+1}$, $u_1^{n+1}$ and $\lambda_1^{n+1}$ using (2.7). Moreover, by design, each problem (2.6), (2.7) can be split into
two sub-problems.
### 3. ESTIMATES FOR THE PREDICTION STEP

The key to the analysis of the correction method is the error estimates of the time difference quantities of the prediction step (2.6). The proof of these results are established in [15]. In this section, we summarize these estimates. We define the errors

$$ (3.1) \qquad \begin{aligned} W_0^n &= w^n - w_0^n, \\ U_0^n &= u^n - u_0^n, \\ \Lambda_0^n &= l^n - \lambda_0^n. \end{aligned} $$

In [15, Corollary 4.2], we show the following estimates for the prediction step (2.6).

$$ (3.2) \qquad \begin{aligned} & \max_{0 \le n \le N-1} \|W_0^{n+1} - W_0^n\|_{L^2(\Omega_s)}^2 + \max_{0 \le n \le N-1} \|U_0^{n+1} - U_0^n\|_{L^2(\Omega_f)}^2 \\ & \quad + \Delta t \sum_{n=0}^{N-1} \nu_f \|\nabla(U_0^{n+1} - U_0^n)\|_{L^2(\Omega_f)}^2 + \Delta t \sum_{n=0}^{N-1} \nu_s \|\nabla(W_0^{n+1} - W_0^n)\|_{L^2(\Omega_s)}^2 \le C(\Delta t)^4 \mathfrak{Y}, \end{aligned} $$

Additionally, we obtain the following estimate in [15, Corollary 4.2],

$$ (3.3) \qquad \Delta t \sum_{n=0}^{M-1} \nu_f \left\| \frac{\nabla(U_0^{n+1} - 2U_0^n + U_0^{n-1})}{(\Delta t)^2} \right\|_{L^2(\Omega_f)}^2 \le C(\Delta t)^2 \mathfrak{Y}. $$

Here $\mathfrak{Y}$ and $\mathfrak{Z}$ are defined in Appendix A. At last, the following relation is established in [15] as well,

$$ (3.4) \qquad \Lambda_0^{n+1} = \nu_f \nabla U_0^{n+1} \cdot \mathbf{n}_f \quad \text{on } \Sigma. $$

We then make the following crucial assumption.

**Assumption 3.1.** We have,

$$ (3.5) \qquad \max_{0 \le n \le N-1} \|\nabla(U_0^{n+1} - U_0^n)\|_{L^2(\Omega_f)}^2 + \Delta t \sum_{n=0}^{N-1} \|D^2(U_0^{n+1} - U_0^n)\|_{L^2(\Omega_f)}^2 \le C(\Delta t)^4 \mathfrak{Y}, $$

where $\mathfrak{Y}$ contains the regularity of the solutions.

Due to Assumption 3.1, the following Corollary is valid.

**Corollary 3.2.** We have,

$$ (3.6) \qquad \Delta t \sum_{n=0}^{N-1} \|\Lambda_0^{n+1} - \Lambda_0^n\|_{L^2(\Sigma)}^2 \le C(\Delta t)^4 \mathfrak{Y} $$

for any $n = 0, 1, 2, \dots, N-1$.

*Proof.* By a standard trace inequality, the relation (3.4) and (3.5), the estimate (3.6) immediately follows. □

*Remark 3.3.* If we consider a simple setting where the interface is perpendicular to the two sides of the domain $\Omega$, as depicted in Figure 2, we can prove Assumption 3.1 with $\mathfrak{Y} = \mathfrak{Y} + \mathfrak{Z} + \mathfrak{Y} + \nu_f \|\partial_t^3 \mathbf{w}\|_{L^2((0,T),L^2(\Omega_s))}^2$. Here $\mathfrak{Y}$ and $\mathfrak{Z}$ are defined in Appendix A. We refer to [15] for the proof of Assumption 3.1 in this special case. However, it is nontrivial to prove Assumption 3.1 for a more general case. Nonetheless, numerical results in [15] indicate that Assumption 3.1 holds for general cases.

### 4. ERROR ANALYSIS FOR THE CORRECTION STEP

In this section we provide the error analysis of (2.7) under Assumption 3.1. As argued in [14], the prediction step is first-order accurate. In this section, we will prove that the correction method produces a second-order accurate method provided Assumption 3.1 is satisfied.
### 3. ESTIMATES FOR THE PREDICTION STEP

The key to the analysis of the correction method is the error estimates of the time difference quantities of the prediction step (2.6). The proof of these results are established in [15]. In this section, we summarize these estimates. We define the errors

$$ (3.1) \qquad \begin{aligned} W_0^n &= w^n - w_0^n, \\ U_0^n &= u^n - u_0^n, \\ \Lambda_0^n &= l^n - \lambda_0^n. \end{aligned} $$

In [15, Corollary 4.2], we show the following estimates for the prediction step (2.6).

$$ (3.2) \qquad \begin{aligned} & \max_{0 \le n \le N-1} \|W_0^{n+1} - W_0^n\|_{L^2(\Omega_s)}^2 + \max_{0 \le n \le N-1} \|U_0^{n+1} - U_0^n\|_{L^2(\Omega_f)}^2 \\ & \quad + \Delta t \sum_{n=0}^{N-1} \nu_f \|\nabla(U_0^{n+1} - U_0^n)\|_{L^2(\Omega_f)}^2 + \Delta t \sum_{n=0}^{N-1} \nu_s \|\nabla(W_0^{n+1} - W_0^n)\|_{L^2(\Omega_s)}^2 \le C(\Delta t)^4 \mathfrak{Y}, \end{aligned} $$

Additionally, we obtain the following estimate in [15, Corollary 4.2],

$$ (3.3) \qquad \Delta t \sum_{n=0}^{M-1} \nu_f \left\| \frac{\nabla(U_0^{n+1} - 2U_0^n + U_0^{n-1})}{(\Delta t)^2} \right\|_{L^2(\Omega_f)}^2 \le C(\Delta t)^2 \mathfrak{Y}. $$

Here $\mathfrak{Y}$ and $\mathfrak{Z}$ are defined in Appendix A. At last, the following relation is established in [15] as well,

$$ (3.4) \qquad \Lambda_0^{n+1} = \nu_f \nabla U_0^{n+1} \cdot n_f \quad \text{on } \Sigma. $$

We then make the following crucial assumption.

**Assumption 3.1.** We have,

$$ (3.5) \qquad \max_{0 \le n \le N-1} \|\nabla(U_0^{n+1} - U_0^n)\|_{L^2(\Omega_f)}^2 + \Delta t \sum_{n=0}^{N-1} \|D^2(U_0^{n+1} - U_0^n)\|_{L^2(\Omega_f)}^2 \le C(\Delta t)^4 \mathfrak{Y}, $$

where $\mathfrak{Y}$ contains the regularity of the solutions.

Due to Assumption 3.1, the following Corollary is valid.

**Corollary 3.2.** We have,

$$ (3.6) \qquad \Delta t \sum_{n=0}^{N-1} \|\Lambda_0^{n+1} - \Lambda_0^n\|_{L^2(\Sigma)}^2 \le C(\Delta t)^4 \mathfrak{Y} $$

for any $n = 0, 1, 2, \dots, N-1$.

*Proof.* By a standard trace inequality, the relation (3.4) and (3.5), the estimate (3.6) immediately follows. □

*Remark 3.3.* If we consider a simple setting where the interface is perpendicular to the two sides of the domain $\Omega$, as depicted in Figure 2, we can prove Assumption 3.1 with $\mathfrak{Y} = \mathfrak{Y} + \mathfrak{Z} + \mathfrak{Y} + \nu_f \|\partial_t^3 w\|_{L^2((0,T),L^2(\Omega_s))}^2$. Here $\mathfrak{Y}$ and $\mathfrak{Z}$ are defined in Appendix A. We refer to [15] for the proof of Assumption 3.1 in this special case. However, it is nontrivial to prove Assumption 3.1 for a more general case. Nonetheless, numerical results in [15] indicate that Assumption 3.1 holds for general cases.

### 4. ERROR ANALYSIS FOR THE CORRECTION STEP

In this section we provide the error analysis of (2.7) under Assumption 3.1. As argued in [14], the prediction step is first-order accurate. In this section, we will prove that the correction method produces a second-order accurate method provided Assumption 3.1 is satisfied.
<figure><img src="image_2.png" /><figcaption>FIGURE 2. The domains $\Omega_f$ and $\Omega_s$ with horizontal interface $\Sigma$</figcaption></figure>

4.1. **Error equations.** To obtain error equations, we write the exact solution in a way that can be compared to the numerical method.

$$ (4.1a) \quad (\partial_{\Delta t} \mathbf{w}^{n+1}, z)_s + \nu_s(\nabla \mathbf{w}^{n+1}, \nabla z)_s + \alpha \langle \mathbf{w}^{n+1} - \mathbf{u}^n, z \rangle + \langle |^n, z \rangle \\ = \nu_s(\nabla \mathbf{w}^{n+1}, \nabla z)_s + \alpha \langle \mathbf{w}^{n+1} - \mathbf{u}^n, z \rangle + \langle |^n, z \rangle \\ - \nu_s(\nabla \mathbf{w}^{n+\frac{1}{2}}, \nabla z)_s - \langle |^{n+\frac{1}{2}}, z \rangle + T_1(z), \quad \forall z \in V_s, $$

$$ (4.1b) \quad (\partial_{\Delta t} \mathbf{u}^{n+1}, v)_f + \nu_f(\nabla \mathbf{u}^{n+1}, \nabla v)_f - \langle |^{n+1}, v \rangle \\ = \nu_f(\nabla \mathbf{u}^{n+1}, \nabla v)_f - \langle |^{n+1}, v \rangle - \nu_f(\nabla \mathbf{u}^{n+\frac{1}{2}}, \nabla v)_f + \langle |^{n+\frac{1}{2}}, v \rangle + T_2(v), \quad \forall v \in V_f, $$

$$ (4.1c) \quad \langle (|^{n+1} - |^n) + \alpha(\mathbf{u}^{n+1} - \mathbf{w}^{n+1}), \mu \rangle = \langle |^{n+1} - |^n, \mu \rangle, \quad \forall \mu \in V_g, $$

where

$$ (4.2a) \quad T_1(z) := (\partial_{\Delta t} \mathbf{w}^{n+1}, z)_s + \nu_s(\nabla \mathbf{w}^{n+\frac{1}{2}}, \nabla z)_s + \langle |^{n+\frac{1}{2}}, z \rangle, $$

$$ (4.2b) \quad T_2(v) := (\partial_{\Delta t} \mathbf{u}^{n+1}, v)_f + \nu_f(\nabla \mathbf{u}^{n+\frac{1}{2}}, \nabla v)_f - \langle |^{n+\frac{1}{2}}, v \rangle. $$

Now we define the errors

$$ W_1^n = \mathbf{w}^n - w_1^n, $$

$$ (4.3) \quad U_1^n = \mathbf{u}^n - u_1^n, $$

$$ \Lambda_1^n = |^n - \lambda_1^n. $$

Subtract (2.7) and (4.1), we obtain the error equations for $n \ge 0$,

$$ (4.4a) \quad (\partial_{\Delta t} W_1^{n+1}, z)_s + \nu_s(\nabla W_1^{n+1}, \nabla z)_s + \alpha \langle W_1^{n+1} - U_1^n, z \rangle + \langle \Lambda_1^n, z \rangle = R_1(z) + T_1(z) \quad \forall z \in V_s, $$

$$ (4.4b) \quad (\partial_{\Delta t} U_1^{n+1}, v)_f + \nu_f(\nabla U_1^{n+1}, \nabla v)_f - \langle \Lambda_1^{n+1}, v \rangle = R_2(v) + T_2(v) \quad \forall v \in V_f, $$

$$ (4.4c) \quad \langle \Lambda_1^{n+1} - \Lambda_1^n + \alpha(U_1^{n+1} - W_1^{n+1}), \mu \rangle = \langle \Lambda_0^{n+1} - \Lambda_0^n, \mu \rangle \quad \forall \mu \in V_g. $$

Here

$$ (4.5a) \quad R_1(z) := \frac{1}{2} \nu_s (\nabla (W_0^{n+1} - W_0^n), \nabla z)_s + \alpha \langle (W_0^{n+1} - W_0^n), z \rangle - \frac{1}{2} \langle \Lambda_0^{n+1} - \Lambda_0^n, z \rangle, $$

$$ (4.5b) \quad R_2(v) := \frac{1}{2} \nu_f (\nabla (U_0^{n+1} - U_0^n), \nabla v)_f - \frac{1}{2} \langle \Lambda_0^{n+1} - \Lambda_0^n, v \rangle. $$

It follows from (4.4c) that,

$$ (4.6) \quad \Lambda_1^{n+1} - \Lambda_1^n + \alpha(U_1^{n+1} - W_1^{n+1}) = \Lambda_0^{n+1} - \Lambda_0^n \quad \text{on } \Sigma. $$
<figure><img src="image_3.png" /><figcaption>FIGURE 2. The domains $\Omega_f$ and $\Omega_s$ with horizontal interface $\Sigma$</figcaption></figure>

4.1. **Error equations.** To obtain error equations, we write the exact solution in a way that can be compared to the numerical method.

$$
\begin{align}
(4.1\text{a}) \quad & (\partial_{\Delta t} \mathbf{w}^{n+1}, z)_s + \nu_s (\nabla \mathbf{w}^{n+1}, \nabla z)_s + \alpha \langle \mathbf{w}^{n+1} - \mathbf{u}^n, z \rangle + \langle |^n, z \rangle \nonumber \\
& = \nu_s (\nabla \mathbf{w}^{n+1}, \nabla z)_s + \alpha \langle \mathbf{w}^{n+1} - \mathbf{u}^n, z \rangle + \langle |^n, z \rangle \nonumber \\
& \phantom{=} - \nu_s (\nabla \mathbf{w}^{n+\frac{1}{2}}, \nabla z)_s - \langle |^{n+\frac{1}{2}}, z \rangle + T_1(z), \quad \forall z \in V_s, \nonumber
\end{align}
$$

$$
\begin{align}
(4.1\text{b}) \quad & (\partial_{\Delta t} \mathbf{u}^{n+1}, v)_f + \nu_f (\nabla \mathbf{u}^{n+1}, \nabla v)_f - \langle |^{n+1}, v \rangle \nonumber \\
& = \nu_f (\nabla \mathbf{u}^{n+1}, \nabla v)_f - \langle |^{n+1}, v \rangle - \nu_f (\nabla \mathbf{u}^{n+\frac{1}{2}}, \nabla v)_f + \langle |^{n+\frac{1}{2}}, v \rangle + T_2(v), \quad \forall v \in V_f, \nonumber
\end{align}
$$

$$
(4.1c) \quad \langle (|^{n+1} - |^n) + \alpha(\mathbf{u}^{n+1} - \mathbf{w}^{n+1}), \mu \rangle = \langle |^{n+1} - |^n, \mu \rangle, \quad \forall \mu \in V_g,
$$

where

$$
(4.2\text{a}) \quad T_1(z) := (\partial_{\Delta t} \mathbf{w}^{n+1}, z)_s + \nu_s (\nabla \mathbf{w}^{n+\frac{1}{2}}, \nabla z)_s + \langle |^{n+\frac{1}{2}}, z \rangle,
$$

$$
(4.2\text{b}) \quad T_2(v) := (\partial_{\Delta t} \mathbf{u}^{n+1}, v)_f + \nu_f (\nabla \mathbf{u}^{n+\frac{1}{2}}, \nabla v)_f - \langle |^{n+\frac{1}{2}}, v \rangle.
$$

Now we define the errors

$$
W_1^n = \mathbf{w}^n - w_1^n,
$$

$$
(4.3) \qquad U_1^n = \mathbf{u}^n - u_1^n,
$$

$$
\Lambda_1^n = |^n - \lambda_1^n.
$$

Subtract (2.7) and (4.1), we obtain the error equations for $n \ge 0$,

$$
(4.4\text{a}) \quad (\partial_{\Delta t} W_1^{n+1}, z)_s + \nu_s (\nabla W_1^{n+1}, \nabla z)_s + \alpha \langle W_1^{n+1} - U_1^n, z \rangle + \langle \Lambda_1^n, z \rangle = R_1(z) + T_1(z) \quad \forall z \in V_s,
$$

$$
(4.4\text{b}) \quad (\partial_{\Delta t} U_1^{n+1}, v)_f + \nu_f (\nabla U_1^{n+1}, \nabla v)_f - \langle \Lambda_1^{n+1}, v \rangle = R_2(v) + T_2(v) \quad \forall v \in V_f,
$$

$$
(4.4c) \quad \langle \Lambda_1^{n+1} - \Lambda_1^n + \alpha(U_1^{n+1} - W_1^{n+1}), \mu \rangle = \langle \Lambda_0^{n+1} - \Lambda_0^n, \mu \rangle \quad \forall \mu \in V_g.
$$

Here

$$
(4.5\text{a}) \quad R_1(z) := \frac{1}{2} \nu_s (\nabla (W_0^{n+1} - W_0^n), \nabla z)_s + \alpha \langle (W_0^{n+1} - W_0^n), z \rangle - \frac{1}{2} \langle \Lambda_0^{n+1} - \Lambda_0^n, z \rangle,
$$

$$
(4.5\text{b}) \quad R_2(v) := \frac{1}{2} \nu_f (\nabla (U_0^{n+1} - U_0^n), \nabla v)_f - \frac{1}{2} \langle \Lambda_0^{n+1} - \Lambda_0^n, v \rangle.
$$

It follows from (4.4c) that,

$$
(4.6) \qquad \Lambda_1^{n+1} - \Lambda_1^n + \alpha(U_1^{n+1} - W_1^{n+1}) = \Lambda_0^{n+1} - \Lambda_0^n \quad \text{on } \Sigma.
$$
We use the convention $u_1^i = u^i = u_0$ and $w_1^i = w^i = w_0$ for $i < 0$.

4.2. **Preliminary estimate for the correction step.** In this section, we establish some preliminary estimates that will be used in the convergence analysis. The techniques we use here is similar but slightly different to those in [15]. First, we define the following quantities for $n \ge 0$,

$$Z_1^{n+1} := \frac{1}{2} \|W_1^{n+1}\|_{L^2(\Omega_s)}^2 + \frac{1}{2} \|U_1^{n+1}\|_{L^2(\Omega_f)}^2 + \frac{\Delta t \alpha}{2} \|U_1^{n+1}\|_{L^2(\Sigma)}^2 + \frac{\Delta t}{2\alpha} \|\Lambda_1^{n+1}\|_{L^2(\Sigma)}^2,$$

$$S_1^{n+1} := \Delta t (\nu_f \|\nabla U_1^{n+1}\|_{L^2(\Omega_f)}^2 + \nu_s \|\nabla W_1^{n+1}\|_{L^2(\Omega_s)}^2) + \frac{1}{2} (\|W_1^{n+1} - W_1^n\|_{L^2(\Omega_s)}^2 + \|U_1^{n+1} - U_1^n\|_{L^2(\Omega_f)}^2) \\ + \frac{\alpha \Delta t}{2} \|(U_1^n - U_1^{n+1}) + \frac{1}{\alpha}(\Lambda_1^n - \Lambda_1^{n+1})\|_{L^2(\Sigma)}^2.$$

We will use the following identity where $(\cdot, \cdot)$ is an inner product and $\|\cdot\|$ is the corresponding norm:

$$(4.7) \qquad (\varphi - \psi, \vartheta) = \frac{1}{2} (\|\varphi\|^2 - \|\psi\|^2 + \|\vartheta - \psi\|^2 - \|\vartheta - \varphi\|^2).$$

We start with the following lemma.

**Lemma 4.1.** It holds,

$$(4.8) \qquad Z_1^{n+1} + S_1^{n+1} = Z_1^n + \Delta t F_1^{n+1} + \Delta t R^{n+1} + \Delta t T^{n+1} + \frac{\Delta t}{\alpha} \langle \Lambda_1^{n+1}, \Lambda_0^{n+1} - \Lambda_0^n \rangle,$$

where

$$F_1^{n+1} := - \langle W_1^{n+1}, \Lambda_0^n - \Lambda_0^{n+1} \rangle + \langle U_1^n - U_1^{n+1}, \Lambda_0^n - \Lambda_0^{n+1} \rangle,$$

$$R^{n+1} := R_1(W_1^{n+1}) + R_2(U_1^{n+1}),$$

$$T^{n+1} := T_1(W_1^{n+1}) + T_2(U_1^{n+1}).$$

*Proof.* To begin, we set $z = \Delta t W_1^{n+1}$ in (4.4a) and $v = \Delta t U_1^{n+1}$ in (4.4b) and we obtain, by (4.7),

$$\frac{1}{2} \|W_1^{n+1}\|_{L^2(\Omega_s)}^2 + \frac{1}{2} \|U_1^{n+1}\|_{L^2(\Omega_f)}^2 + \frac{1}{2} \|W_1^{n+1} - W_1^n\|_{L^2(\Omega_s)}^2 + \frac{1}{2} \|U_1^{n+1} - U_1^n\|_{L^2(\Omega_f)}^2 \\ + \nu_s \Delta t \|\nabla W_1^{n+1}\|_{L^2(\Omega_s)}^2 + \nu_f \Delta t \|\nabla U_1^{n+1}\|_{L^2(\Omega_f)}^2$$

$$(4.9) \qquad = \frac{1}{2} \|W_1^n\|_{L^2(\Omega_s)}^2 + \frac{1}{2} \|U_1^n\|_{L^2(\Omega_f)}^2 + \Delta t (J^{n+1} + R^{n+1} + T^{n+1}),$$

where

$$(4.10) \qquad J^{n+1} := - \alpha \langle W_1^{n+1} - U_1^n, W_1^{n+1} \rangle - \langle \Lambda_1^n, W_1^{n+1} \rangle + \langle \Lambda_1^{n+1}, U_1^{n+1} \rangle.$$

After some manipulations and use (4.6), we have

$$J^{n+1} = \mathcal{J}^{n+1} - \langle W_1^{n+1}, \Lambda_0^n - \Lambda_0^{n+1} \rangle - \frac{1}{\alpha} \langle \Lambda_1^{n+1}, \Lambda_0^n - \Lambda_0^{n+1} \rangle + \langle U_1^n - U_1^{n+1}, \Lambda_0^n - \Lambda_0^{n+1} \rangle,$$

where

$$J^{n+1} := \alpha \langle U_1^n - U_1^{n+1}, U_1^{n+1} \rangle + \frac{1}{\alpha} \langle \Lambda_1^n - \Lambda_1^{n+1}, \Lambda_1^{n+1} \rangle - \langle U_1^n - U_1^{n+1}, \Lambda_1^n - \Lambda_1^{n+1} \rangle.$$

One can easily show that

$$(4.11) \qquad J^{n+1} = \frac{\alpha}{2} (\|U_1^n\|_{L^2(\Sigma)}^2 - \|U_1^{n+1}\|_{L^2(\Sigma)}^2) + \frac{1}{2\alpha} (\|\Lambda_1^n\|_{L^2(\Sigma)}^2 - \|\Lambda_1^{n+1}\|_{L^2(\Sigma)}^2) \\ - \frac{\alpha}{2} \|(U_1^n - U_1^{n+1}) + \frac{1}{\alpha}(\Lambda_1^n - \Lambda_1^{n+1})\|_{L^2(\Sigma)}^2.$$

Combining (4.11), (4.10) and (4.9) we arrive at the identity (4.8). □

In the following subsections, we give preliminary estimates of all the terms on the right-hand side of (4.8). To this purpose, we will need the following trace inequality, which follows from the continuity of the trace operator and from Poincaré's inequality.
We use the convention $u_1^i = u^i = u_0$ and $w_1^i = w^i = w_0$ for $i < 0$.

4.2. **Preliminary estimate for the correction step.** In this section, we establish some preliminary estimates that will be used in the convergence analysis. The techniques we use here is similar but slightly different to those in [15]. First, we define the following quantities for $n \ge 0$,

$$Z_1^{n+1} := \frac{1}{2} \|W_1^{n+1}\|_{L^2(\Omega_s)}^2 + \frac{1}{2} \|U_1^{n+1}\|_{L^2(\Omega_f)}^2 + \frac{\Delta t \alpha}{2} \|U_1^{n+1}\|_{L^2(\Sigma)}^2 + \frac{\Delta t}{2\alpha} \|\Lambda_1^{n+1}\|_{L^2(\Sigma)}^2,$$

$$S_1^{n+1} := \Delta t (\nu_f \|\nabla U_1^{n+1}\|_{L^2(\Omega_f)}^2 + \nu_s \|\nabla W_1^{n+1}\|_{L^2(\Omega_s)}^2) + \frac{1}{2} (\|W_1^{n+1} - W_1^n\|_{L^2(\Omega_s)}^2 + \|U_1^{n+1} - U_1^n\|_{L^2(\Omega_f)}^2) \\ + \frac{\alpha \Delta t}{2} \|(U_1^n - U_1^{n+1}) + \frac{1}{\alpha}(\Lambda_1^n - \Lambda_1^{n+1})\|_{L^2(\Sigma)}^2.$$

We will use the following identity where $(\cdot, \cdot)$ is an inner product and $\|\cdot\|$ is the corresponding norm:

$$(4.7) \qquad (\varphi - \psi, \vartheta) = \frac{1}{2} (\|\varphi\|^2 - \|\psi\|^2 + \|\vartheta - \psi\|^2 - \|\vartheta - \varphi\|^2).$$

We start with the following lemma.

**Lemma 4.1.** It holds,

$$(4.8) \qquad Z_1^{n+1} + S_1^{n+1} = Z_1^n + \Delta t F_1^{n+1} + \Delta t R^{n+1} + \Delta t T^{n+1} + \frac{\Delta t}{\alpha} \langle \Lambda_1^{n+1}, \Lambda_0^{n+1} - \Lambda_0^n \rangle,$$

where

$$F_1^{n+1} := - \langle W_1^{n+1}, \Lambda_0^n - \Lambda_0^{n+1} \rangle + \langle U_1^n - U_1^{n+1}, \Lambda_0^n - \Lambda_0^{n+1} \rangle,$$

$$R^{n+1} := R_1(W_1^{n+1}) + R_2(U_1^{n+1}),$$

$$T^{n+1} := T_1(W_1^{n+1}) + T_2(U_1^{n+1}).$$

*Proof.* To begin, we set $z = \Delta t W_1^{n+1}$ in (4.4a) and $v = \Delta t U_1^{n+1}$ in (4.4b) and we obtain, by (4.7),

$$\frac{1}{2} \|W_1^{n+1}\|_{L^2(\Omega_s)}^2 + \frac{1}{2} \|U_1^{n+1}\|_{L^2(\Omega_f)}^2 + \frac{1}{2} \|W_1^{n+1} - W_1^n\|_{L^2(\Omega_s)}^2 + \frac{1}{2} \|U_1^{n+1} - U_1^n\|_{L^2(\Omega_f)}^2 \\ + \nu_s \Delta t \|\nabla W_1^{n+1}\|_{L^2(\Omega_s)}^2 + \nu_f \Delta t \|\nabla U_1^{n+1}\|_{L^2(\Omega_f)}^2$$

$$(4.9) \qquad = \frac{1}{2} \|W_1^n\|_{L^2(\Omega_s)}^2 + \frac{1}{2} \|U_1^n\|_{L^2(\Omega_f)}^2 + \Delta t (J^{n+1} + R^{n+1} + T^{n+1}),$$

where

$$(4.10) \qquad J^{n+1} := -\alpha \langle W_1^{n+1} - U_1^n, W_1^{n+1} \rangle - \langle \Lambda_1^n, W_1^{n+1} \rangle + \langle \Lambda_1^{n+1}, U_1^{n+1} \rangle.$$

After some manipulations and use (4.6), we have

$$J^{n+1} = \mathcal{J}^{n+1} - \langle W_1^{n+1}, \Lambda_0^n - \Lambda_0^{n+1} \rangle - \frac{1}{\alpha} \langle \Lambda_1^{n+1}, \Lambda_0^n - \Lambda_0^{n+1} \rangle + \langle U_1^n - U_1^{n+1}, \Lambda_0^n - \Lambda_0^{n+1} \rangle,$$

where

$$J^{n+1} := \alpha \langle U_1^n - U_1^{n+1}, U_1^{n+1} \rangle + \frac{1}{\alpha} \langle \Lambda_1^n - \Lambda_1^{n+1}, \Lambda_1^{n+1} \rangle - \langle U_1^n - U_1^{n+1}, \Lambda_1^n - \Lambda_1^{n+1} \rangle.$$

One can easily show that

$$(4.11) \qquad J^{n+1} = \frac{\alpha}{2} (\|U_1^n\|_{L^2(\Sigma)}^2 - \|U_1^{n+1}\|_{L^2(\Sigma)}^2) + \frac{1}{2\alpha} (\|\Lambda_1^n\|_{L^2(\Sigma)}^2 - \|\Lambda_1^{n+1}\|_{L^2(\Sigma)}^2) \\ - \frac{\alpha}{2} \|(U_1^n - U_1^{n+1}) + \frac{1}{\alpha}(\Lambda_1^n - \Lambda_1^{n+1})\|_{L^2(\Sigma)}^2.$$

Combining (4.11), (4.10) and (4.9) we arrive at the identity (4.8). □

In the following subsections, we give preliminary estimates of all the terms on the right-hand side of (4.8). To this purpose, we will need the following trace inequality, which follows from the continuity of the trace operator and from Poincaré's inequality.
**Proposition 4.2.** There exists constants $C_{tr}$ such that

$$ (4.12) \qquad \|v\|_{L^2(\Sigma)} \le C_{tr} \|\nabla v\|_{L^2(\Omega_i)} \quad \forall v \in V_i \quad i = s, f. $$

4.3. Estimate for $F_1^{n+1}$.

**Lemma 4.3.** Let $1 \le M \le N$, then

$$ (4.13) \qquad \Delta t \sum_{n=0}^{M-1} F_1^{n+1} \le \frac{1}{8} \sum_{n=0}^{M-1} S_1^{n+1} + CD_1(M), $$

where

$$ D_1(M) := \Delta t \sum_{n=0}^{M-1} \left( \frac{C_{tr}^2}{\nu_f} + \frac{C_{tr}^2}{\nu_s} \right) \| \Lambda_0^n - \Lambda_0^{n+1} \|_{L^2(\Sigma)}^2. $$

*Proof.* It follows from (4.12) and the definition of $S_1^{n+1}$ that,

$$ \begin{aligned} \Delta t \sum_{n=0}^{M-1} & (-\langle W_1^{n+1}, \Lambda_0^n - \Lambda_0^{n+1} \rangle + \langle U_1^n - U_1^{n+1}, \Lambda_0^n - \Lambda_0^{n+1} \rangle) \\ & \le \frac{\Delta t}{16} \sum_{n=0}^{M-1} \nu_s \|\nabla W_1^{n+1}\|_{L^2(\Omega_s)}^2 + \frac{\Delta t}{16} \sum_{n=0}^{M-1} \nu_f \|\nabla U_1^{n+1}\|_{L^2(\Omega_f)}^2 \\ & \qquad + C \Delta t \sum_{n=0}^{M-1} \left( \frac{C_{tr}^2}{\nu_f} + \frac{C_{tr}^2}{\nu_s} \right) \|\Lambda_0^n - \Lambda_0^{n+1}\|_{L^2(\Sigma)}^2 \le \frac{1}{8} \sum_{n=0}^{M-1} S_1^{n+1} + CD_1(M), \end{aligned} $$

which completes the proof. $\square$

4.4. Estimate for $R^{n+1}$.

**Lemma 4.4.** Let $1 \le M \le N$, then

$$ (4.14) \qquad \Delta t \sum_{n=0}^{M-1} R^{n+1} \le \frac{1}{8} \sum_{n=0}^{M-1} S_1^{n+1} + CD_2(M), $$

where

$$ D_2(M) := \Delta t \sum_{n=0}^{M-1} \left( \frac{C_{tr}^2}{\nu_f} + \frac{C_{tr}^2}{\nu_s} \right) \| \Lambda_0^n - \Lambda_0^{n+1} \|_{L^2(\Sigma)}^2 + \Delta t \sum_{n=0}^{M-1} \nu_s \| \nabla (W_0^{n+1} - W_0^n) \|_{L^2(\Omega_s)}^2 \\ + \Delta t \sum_{n=0}^{M-1} \frac{C_{tr}^2}{\nu_s} \alpha^2 \| W_0^{n+1} - W_0^n \|_{L^2(\Sigma)}^2 + \Delta t \sum_{n=0}^{M-1} \nu_f \| \nabla (U_0^{n+1} - U_0^n) \|_{L^2(\Omega_f)}^2. $$
**Proposition 4.2.** There exists constants $C_{tr}$ such that

$$ (4.12) \qquad \|v\|_{L^2(\Sigma)} \le C_{tr} \|\nabla v\|_{L^2(\Omega_i)} \quad \forall v \in V_i \quad i = s, f. $$

### 4.3. Estimate for $F_1^{n+1}$.

**Lemma 4.3.** Let $1 \le M \le N$, then

$$ (4.13) \qquad \Delta t \sum_{n=0}^{M-1} F_1^{n+1} \le \frac{1}{8} \sum_{n=0}^{M-1} S_1^{n+1} + CD_1(M), $$

where

$$ D_1(M) := \Delta t \sum_{n=0}^{M-1} \left( \frac{C_{tr}^2}{\nu_f} + \frac{C_{tr}^2}{\nu_s} \right) \|\Lambda_0^n - \Lambda_0^{n+1}\|_{L^2(\Sigma)}^2. $$

*Proof.* It follows from (4.12) and the definition of $S_1^{n+1}$ that,

$$ \begin{aligned} \Delta t \sum_{n=0}^{M-1} & (-\langle W_1^{n+1}, \Lambda_0^n - \Lambda_0^{n+1} \rangle + \langle U_1^n - U_1^{n+1}, \Lambda_0^n - \Lambda_0^{n+1} \rangle) \\ & \le \frac{\Delta t}{16} \sum_{n=0}^{M-1} \nu_s \|\nabla W_1^{n+1}\|_{L^2(\Omega_s)}^2 + \frac{\Delta t}{16} \sum_{n=0}^{M-1} \nu_f \|\nabla U_1^{n+1}\|_{L^2(\Omega_f)}^2 \\ & \quad + C \Delta t \sum_{n=0}^{M-1} \left( \frac{C_{tr}^2}{\nu_f} + \frac{C_{tr}^2}{\nu_s} \right) \|\Lambda_0^n - \Lambda_0^{n+1}\|_{L^2(\Sigma)}^2 \le \frac{1}{8} \sum_{n=0}^{M-1} S_1^{n+1} + CD_1(M), \end{aligned} $$

which completes the proof. $\square$

### 4.4. Estimate for $R^{n+1}$.

**Lemma 4.4.** Let $1 \le M \le N$, then

$$ (4.14) \qquad \Delta t \sum_{n=0}^{M-1} R^{n+1} \le \frac{1}{8} \sum_{n=0}^{M-1} S_1^{n+1} + CD_2(M), $$

where

$$ D_2(M) := \Delta t \sum_{n=0}^{M-1} \left( \frac{C_{tr}^2}{\nu_f} + \frac{C_{tr}^2}{\nu_s} \right) \|\Lambda_0^n - \Lambda_0^{n+1}\|_{L^2(\Sigma)}^2 + \Delta t \sum_{n=0}^{M-1} \nu_s \|\nabla(W_0^{n+1} - W_0^n)\|_{L^2(\Omega_s)}^2 \\ + \Delta t \sum_{n=0}^{M-1} \frac{C_{tr}^2}{\nu_s} \alpha^2 \|W_0^{n+1} - W_0^n\|_{L^2(\Sigma)}^2 + \Delta t \sum_{n=0}^{M-1} \nu_f \|\nabla(U_0^{n+1} - U_0^n)\|_{L^2(\Omega_f)}^2. $$
*Proof.* From (4.12) and the definition of $S_1^{n+1}$, it follows that

$$
\begin{align*}
& \Delta t \sum_{n=0}^{M-1} R_1(W_1^{n+1}) + R_2(U^{n+1}) \\
&= \Delta t \sum_{n=0}^{M-1} \frac{1}{2} \nu_s (\nabla(W_0^{n+1} - W_0^n), \nabla W_1^{n+1})_s + \alpha \langle (W_0^{n+1} - W_0^n), W_1^{n+1} \rangle - \frac{1}{2} \langle \Lambda_0^{n+1} - \Lambda_0^n, W_1^{n+1} \rangle \\
&\quad + \Delta t \sum_{n=0}^{M-1} \frac{1}{2} \nu_f (\nabla(U_0^{n+1} - U_0^n), \nabla U_1^{n+1})_f - \frac{1}{2} \langle \Lambda_0^{n+1} - \Lambda_0^n, U_1^{n+1} \rangle \\
&\le \frac{\Delta t}{8} \sum_{n=0}^{M-1} \nu_s \|\nabla W_1^{n+1}\|_{L^2(\Omega_s)}^2 + C \Delta t \sum_{n=0}^{M-1} \nu_s \|\nabla(W_0^{n+1} - W_0^n)\|_{L^2(\Omega_s)}^2 \\
&\quad + C \Delta t \frac{C_{tr}^2}{\nu_s} \alpha^2 \sum_{n=0}^{M-1} \|W_0^{n+1} - W_0^n\|_{L^2(\Sigma)}^2 + C \Delta t \frac{C_{tr}^2}{\nu_s} \sum_{n=0}^{M-1} \|\Lambda_0^{n+1} - \Lambda_0^n\|_{L^2(\Sigma)}^2 \\
&\quad + \frac{\Delta t}{8} \sum_{n=0}^{M-1} \nu_f \|\nabla U_1^{n+1}\|_{L^2(\Omega_f)}^2 + C \Delta t \sum_{n=0}^{M-1} \nu_f \|\nabla(U_0^{n+1} - U_0^n)\|_{L^2(\Omega_f)}^2 \\
&\quad + C \Delta t \frac{C_{tr}^2}{\nu_f} \sum_{n=0}^{M-1} \|\Lambda_0^{n+1} - \Lambda_0^n\|_{L^2(\Sigma)}^2 \le \frac{1}{8} \sum_{n=0}^{M-1} S_1^{n+1} + CD_2(M),
\end{align*}
$$

which yields (4.14) and completes the proof.

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

### 4.5. Estimate for $T^{n+1}$.

**Lemma 4.5.** Let $1 \le M \le N$, then

$$ (4.15) \qquad \Delta t \sum_{n=0}^{M-1} T^{n+1} \le \frac{1}{4} \max_{1 \le n \le M} Z_1^n + \frac{1}{8} \sum_{n=0}^{M-1} S_1^{n+1} + CD_3(M), $$

where

$$
\begin{align*}
D_3(M) := \Delta t T \sum_{n=0}^{M-1} & \left( \|\partial_{\Delta t} w^{n+1} - \partial_t w(t_{n+\frac{1}{2}})\|_{L^2(\Omega_s)}^2 + \|\partial_{\Delta t} u^{n+1} - \partial_t u(t_{n+\frac{1}{2}})\|_{L^2(\Omega_f)}^2 \right) \\
& + \Delta t \sum_{n=0}^{M-1} \left( \nu_s \|\nabla w^{n+\frac{1}{2}} - \nabla w(t_{n+\frac{1}{2}})\|_{L^2(\Omega_s)}^2 + \nu_f \|\nabla u^{n+\frac{1}{2}} - \nabla u(t_{n+\frac{1}{2}})\|_{L^2(\Omega_f)}^2 \right) \\
& + \Delta t \sum_{n=0}^{M-1} \left( \left( \frac{C_{tr}^2}{\nu_f} + \frac{C_{tr}^2}{\nu_s} \right) \|l^{n+\frac{1}{2}} - l(t_{n+\frac{1}{2}})\|_{L^2(\Sigma)}^2 \right).
\end{align*}
$$

*Proof.* We have, by (4.2),

$$
\begin{align}
\Delta t \sum_{n=0}^{M-1} T^{n+1} &= \Delta t \sum_{n=0}^{M-1} \left[ T_1^{n+1}(W_1^{n+1}) + T_2^{n+1}(U_1^{n+1}) \right] \nonumber \\
(4.16) \qquad &= \Delta t \sum_{n=0}^{M-1} \left[ \begin{aligned}
& (\partial_{\Delta t} w^{n+1}, W_1^{n+1})_s + \nu_s (\nabla w^{n+\frac{1}{2}}, \nabla W_1^{n+1})_s + \langle l^{n+\frac{1}{2}}, W_1^{n+1} \rangle \\
& + (\partial_{\Delta t} u^{n+1}, U_1^{n+1})_f + \nu_s (\nabla u^{n+\frac{1}{2}}, \nabla U_1^{n+1})_f - \langle l^{n+\frac{1}{2}}, U_1^{n+1} \rangle
\end{aligned} \right]. \nonumber
\end{align}
$$

It follows from (2.4) at $t_{n+\frac{1}{2}} := \frac{t_{n+1}+t_n}{2}$ that,

$$ (4.17) \qquad (\partial_t w(t_{n+\frac{1}{2}}), W_1^{n+1})_s + \nu_s (\nabla w(t_{n+\frac{1}{2}}), \nabla W_1^{n+1})_s + \langle l(t_{n+\frac{1}{2}}), W_1^{n+1} \rangle = 0. $$
Similarly, it follows from (2.3) that,

$$ (4.18) \qquad (\partial_t u(t_{n+\frac{1}{2}}), U_1^{n+1})_f + \nu_f (\nabla u(t_{n+\frac{1}{2}}), \nabla U_1^{n+1})_f - \langle l(t_{n+\frac{1}{2}}), U_1^{n+1} \rangle = 0. $$

Subtract (4.17) and (4.18) from (4.16), we obtain

$$ (4.19) \qquad \begin{aligned} \Delta t \sum_{n=0}^{M-1} T^{n+1} &= \Delta t \sum_{n=0}^{M-1} \left[ T_1^{n+1}(W_1^{n+1}) + T_2^{n+1}(U_1^{n+1}) \right] \\ &= \Delta t \sum_{n=0}^{M-1} \left[ (\partial_{\Delta t} w^{n+1} - \partial_t w(t_{n+\frac{1}{2}}), W_1^{n+1})_s + \nu_s (\nabla w^{n+\frac{1}{2}} - \nabla w(t_{n+\frac{1}{2}}), \nabla W_1^{n+1})_s \right. \\ & \qquad \left. + \langle l^{n+\frac{1}{2}} - l(t_{n+\frac{1}{2}}), W_1^{n+1} \rangle \right. \\ & \qquad \left. + (\partial_{\Delta t} u^{n+1} - \partial_t u(t_{n+\frac{1}{2}}), U_1^{n+1})_f + \nu_s (\nabla u^{n+\frac{1}{2}} - \nabla u(t_{n+\frac{1}{2}}), \nabla U_1^{n+1})_f \right. \\ & \qquad \left. - \langle l^{n+\frac{1}{2}} - l(t_{n+\frac{1}{2}}), U_1^{n+1} \rangle \right]. \end{aligned} $$

Consequently, we estimate (4.19) by Cauchy-Schwarz inequality and Young's inequality as follows,

$$ \begin{aligned} \Delta t \sum_{n=0}^{M-1} T^{n+1} &\le \frac{\Delta t}{8T} \sum_{n=0}^{M-1} \left( \|W_1^{n+1}\|_{L^2(\Omega_s)}^2 + \|U_1^{n+1}\|_{L^2(\Omega_f)}^2 \right) \\ &\quad + C \Delta t T \sum_{n=0}^{M-1} \left( \|\partial_{\Delta t} w^{n+1} - \partial_t w(t_{n+\frac{1}{2}})\|_{L^2(\Omega_s)}^2 + \|\partial_{\Delta t} u^{n+1} - \partial_t u(t_{n+\frac{1}{2}})\|_{L^2(\Omega_f)}^2 \right) \\ &\quad + \frac{\Delta t}{8} \sum_{n=0}^{M-1} \left( \nu_s \|\nabla W_1^{n+1}\|_{L^2(\Omega_s)}^2 + \nu_f \|\nabla U_1^{n+1}\|_{L^2(\Omega_f)}^2 \right) \\ &\quad + C \Delta t \sum_{n=0}^{M-1} \left( \nu_s \|\nabla w^{n+\frac{1}{2}} - \nabla w(t_{n+\frac{1}{2}})\|_{L^2(\Omega_s)}^2 + \nu_f \|\nabla u^{n+\frac{1}{2}} - \nabla u(t_{n+\frac{1}{2}})\|_{L^2(\Omega_f)}^2 \right) \\ &\quad + C \Delta t \sum_{n=0}^{M-1} \left( \left( \frac{C_{\text{tr}}^2}{\nu_f} + \frac{C_{\text{tr}}^2}{\nu_s} \right) \|l^{n+\frac{1}{2}} - l(t_{n+\frac{1}{2}})\|_{L^2(\Sigma)}^2 \right), \end{aligned} $$

which gives (4.15) and hence completes the proof. □

4.6. **Estimate for** $\frac{\Delta t}{\alpha} \langle \Lambda_1^{n+1}, \Lambda_0^{n+1} - \Lambda_0^n \rangle$. In this section we estimate the term

$$ (4.20) \qquad \frac{\Delta t}{\alpha} \sum_{n=0}^{M-1} \langle \Lambda_1^{n+1}, \Lambda_0^{n+1} - \Lambda_0^n \rangle. $$

Note that $\nabla U_0^{n+1} \cdot n_f$ is well-defined on $\Omega_f$. Hence we may extend $\nabla U_0^{n+1} \cdot n_f$ into the domain $\Omega_f$ in a similar fashion as in [15]. Let $\phi : \Omega_f \to \mathbb{R}$ be an affine function that is 1 on $\Sigma$ and vanishes on $\Gamma_D^f$. Therefore, the following properties are satisfied:

$$ (4.21a) \qquad \phi \in V_f, $$

$$ (4.21b) \qquad 0 \le \phi \le 1, $$

$$ (4.21c) \qquad \phi = 1 \text{ on } \Sigma, $$

$$ (4.21d) \qquad \|\nabla \phi\|_{L^2(\Omega_f)} \le C. $$

Define the function

$$ (4.22) \qquad \tilde{\Lambda}_0^{n+1} := \phi(\nu_f \nabla U_0^{n+1} \cdot n_f). $$
Similarly, it follows from (2.3) that,

$$ (4.18) \qquad (\partial_t u(t_{n+\frac{1}{2}}), U_1^{n+1})_f + \nu_f (\nabla u(t_{n+\frac{1}{2}}), \nabla U_1^{n+1})_f - \langle l(t_{n+\frac{1}{2}}), U_1^{n+1} \rangle = 0. $$

Subtract (4.17) and (4.18) from (4.16), we obtain

$$ (4.19) \qquad \begin{aligned} \Delta t \sum_{n=0}^{M-1} T^{n+1} &= \Delta t \sum_{n=0}^{M-1} \left[ T_1^{n+1}(W_1^{n+1}) + T_2^{n+1}(U_1^{n+1}) \right] \\ &= \Delta t \sum_{n=0}^{M-1} \left[ (\partial_{\Delta t} w^{n+1} - \partial_t w(t_{n+\frac{1}{2}}), W_1^{n+1})_s + \nu_s (\nabla w^{n+\frac{1}{2}} - \nabla w(t_{n+\frac{1}{2}}), \nabla W_1^{n+1})_s \right. \\ & \qquad \left. + \langle |^{n+\frac{1}{2}} - l(t_{n+\frac{1}{2}}), W_1^{n+1} \rangle \right. \\ & \qquad \left. + (\partial_{\Delta t} u^{n+1} - \partial_t u(t_{n+\frac{1}{2}}), U_1^{n+1})_f + \nu_s (\nabla u^{n+\frac{1}{2}} - \nabla u(t_{n+\frac{1}{2}}), \nabla U_1^{n+1})_f \right. \\ & \qquad \left. - \langle |^{n+\frac{1}{2}} - l(t_{n+\frac{1}{2}}), U_1^{n+1} \rangle \right]. \end{aligned} $$

Consequently, we estimate (4.19) by Cauchy-Schwarz inequality and Young's inequality as follows,

$$ \begin{aligned} \Delta t \sum_{n=0}^{M-1} T^{n+1} &\le \frac{\Delta t}{8T} \sum_{n=0}^{M-1} \left( \|W_1^{n+1}\|_{L^2(\Omega_s)}^2 + \|U_1^{n+1}\|_{L^2(\Omega_f)}^2 \right) \\ &\quad + C \Delta t T \sum_{n=0}^{M-1} \left( \|\partial_{\Delta t} w^{n+1} - \partial_t w(t_{n+\frac{1}{2}})\|_{L^2(\Omega_s)}^2 + \|\partial_{\Delta t} u^{n+1} - \partial_t u(t_{n+\frac{1}{2}})\|_{L^2(\Omega_f)}^2 \right) \\ &\quad + \frac{\Delta t}{8} \sum_{n=0}^{M-1} \left( \nu_s \|\nabla W_1^{n+1}\|_{L^2(\Omega_s)}^2 + \nu_f \|\nabla U_1^{n+1}\|_{L^2(\Omega_f)}^2 \right) \\ &\quad + C \Delta t \sum_{n=0}^{M-1} \left( \nu_s \|\nabla w^{n+\frac{1}{2}} - \nabla w(t_{n+\frac{1}{2}})\|_{L^2(\Omega_s)}^2 + \nu_f \|\nabla u^{n+\frac{1}{2}} - \nabla u(t_{n+\frac{1}{2}})\|_{L^2(\Omega_f)}^2 \right) \\ &\quad + C \Delta t \sum_{n=0}^{M-1} \left( \left( \frac{C_{\text{tr}}^2}{\nu_f} + \frac{C_{\text{tr}}^2}{\nu_s} \right) \| |^{n+\frac{1}{2}} - l(t_{n+\frac{1}{2}}) \|_{L^2(\Sigma)}^2 \right), \end{aligned} $$

which gives (4.15) and hence completes the proof. □

4.6. **Estimate for** $\frac{\Delta t}{\alpha} \langle \Lambda_1^{n+1}, \Lambda_0^{n+1} - \Lambda_0^n \rangle$. In this section we estimate the term

$$ (4.20) \qquad \frac{\Delta t}{\alpha} \sum_{n=0}^{M-1} \langle \Lambda_1^{n+1}, \Lambda_0^{n+1} - \Lambda_0^n \rangle. $$

Note that $\nabla U_0^{n+1} \cdot n_f$ is well-defined on $\Omega_f$. Hence we may extend $\nabla U_0^{n+1} \cdot n_f$ into the domain $\Omega_f$ in a similar fashion as in [15]. Let $\phi : \Omega_f \to \mathbb{R}$ be an affine function that is 1 on $\Sigma$ and vanishes on $\Gamma_D^f$. Therefore, the following properties are satisfied:

$$ (4.21a) \qquad \phi \in V_f, $$

$$ (4.21b) \qquad 0 \le \phi \le 1, $$

$$ (4.21c) \qquad \phi = 1 \text{ on } \Sigma, $$

$$ (4.21d) \qquad \|\nabla \phi\|_{L^2(\Omega_f)} \le C. $$

Define the function

$$ (4.22) \qquad \tilde{\Lambda}_0^{n+1} := \phi(\nu_f \nabla U_0^{n+1} \cdot n_f). $$
Note that $\tilde{\Lambda}_0^{n+1} = \Lambda_0^{n+1}$ on $\Sigma$ and $\tilde{\Lambda}_0^{n+1} \in V_f$, so it is indeed an extension. It then follows that

$$ (4.23) \qquad \langle \Lambda_1^{n+1}, \Lambda_0^{n+1} - \Lambda_0^n \rangle = \langle \Lambda_1^{n+1}, \tilde{\Lambda}_0^{n+1} - \tilde{\Lambda}_0^n \rangle. $$

Denote $\mathcal{L}^{n+1} := \Lambda_0^{n+1} - \Lambda_0^n$ and $\tilde{\mathcal{L}}^{n+1} := \tilde{\Lambda}_0^{n+1} - \tilde{\Lambda}_0^n$. Furthermore, owing to (4.4b) with $v = \tilde{\mathcal{L}}^{n+1}$, we have

$$ (4.24) \qquad \langle \Lambda_1^{n+1}, \tilde{\mathcal{L}}^{n+1} \rangle = (\partial_{\Delta t} U_1^{n+1}, \tilde{\mathcal{L}}^{n+1})_f + \nu_f (\nabla U_1^{n+1}, \nabla \tilde{\mathcal{L}}^{n+1})_f - R_2(\tilde{\mathcal{L}}^{n+1}) - T_2(\tilde{\mathcal{L}}^{n+1}). $$

Therefore, we obtain

$$ (4.25) \qquad \frac{\Delta t}{\alpha} \sum_{n=0}^{M-1} \langle \Lambda_1^{n+1}, \Lambda_0^{n+1} - \Lambda_0^n \rangle = \frac{1}{\alpha} \sum_{n=0}^{M-1} (U_1^{n+1} - U_1^n, \tilde{\mathcal{L}}^{n+1})_f + \frac{\Delta t}{\alpha} \sum_{n=0}^{M-1} \nu_f (\nabla U_1^{n+1}, \nabla \tilde{\mathcal{L}}^{n+1})_f \\ - \frac{\Delta t}{\alpha} \sum_{n=0}^{M-1} R_2(\tilde{\mathcal{L}}^{n+1}) - \frac{\Delta t}{\alpha} \sum_{n=0}^{M-1} T_2(\tilde{\mathcal{L}}^{n+1}). $$

The following lemma states the estimate for (4.20), based on the previous identity.

**Lemma 4.6.** Let $1 \le M \le N$, then

$$ (4.26) \qquad \frac{\Delta t}{\alpha} \sum_{n=0}^{M-1} \langle \Lambda_1^{n+1}, \Lambda_0^{n+1} - \Lambda_0^n \rangle \le \frac{1}{8} \sum_{n=0}^{M-1} S^{n+1} + \frac{1}{4} \max_{1 \le n \le M} Z^n + CD_4(M), $$

where

$$ \begin{aligned} D_4(M) := & \Delta t \frac{\nu_f}{\alpha^2} \sum_{n=0}^{M-1} \|\nabla \tilde{\mathcal{L}}^{n+1}\|_{L^2(\Omega_f)}^2 + \frac{\Delta t}{\alpha} \sum_{n=0}^{M-1} \|\tilde{\mathcal{L}}^{n+1}\|_{L^2(\Omega_f)}^2 \\ & + \frac{T(\Delta t)^3}{\alpha^2} \sum_{n=1}^{M-1} \|\frac{\tilde{\mathcal{L}}^{n+1} - \tilde{\mathcal{L}}^n}{(\Delta t)^2}\|_{L^2(\Omega_f)}^2 + \frac{1}{\alpha^2} \|\tilde{\mathcal{L}}^M\|_{L^2(\Omega_f)}^2 \\ & + \Delta t \frac{C_{tr}^2}{\nu_f} \sum_{n=0}^{M-1} \|\Lambda_0^{n+1} - \Lambda_0^n\|_{L^2(\Sigma)}^2 + \Delta t \sum_{n=0}^{M-1} \nu_f \|\nabla(U_0^{n+1} - U_0^n)\|_{L^2(\Omega_f)}^2 \\ & + \frac{\Delta t}{\alpha} \sum_{n=0}^{M-1} \|\partial_{\Delta t} u^{n+1} - \partial_t u(t_{n+\frac{1}{2}})\|_{L^2(\Omega_f)}^2 + \Delta t \sum_{n=0}^{M-1} \nu_f \|\nabla u^{n+\frac{1}{2}} - \nabla u(t_{n+\frac{1}{2}})\|_{L^2(\Omega_f)}^2 \\ & + \Delta t \sum_{n=0}^{M-1} \frac{C_{tr}^2}{\nu_f} \|u^{n+\frac{1}{2}} - l(t_{n+\frac{1}{2}})\|_{L^2(\Sigma)}^2. \end{aligned} \qquad (4.27) $$

*Proof.* From (4.25), (4.5b) and (4.2b), it follows that

$$ \begin{aligned} & \frac{\Delta t}{\alpha} \sum_{n=0}^{M-1} \langle \Lambda_1^{n+1}, \Lambda_0^{n+1} - \Lambda_0^n \rangle \\ = & \frac{1}{\alpha} \sum_{n=0}^{M-1} (U_1^{n+1} - U_1^n, \tilde{\mathcal{L}}^{n+1})_f + \frac{\Delta t}{\alpha} \sum_{n=0}^{M-1} \nu_f (\nabla U_1^{n+1}, \nabla \tilde{\mathcal{L}}^{n+1})_f \\ & - \frac{\Delta t}{\alpha} \sum_{n=0}^{M-1} \left[ \frac{1}{2} \nu_f (\nabla (U_0^{n+1} - U_0^n), \nabla \tilde{\mathcal{L}}^{n+1})_f - \frac{1}{2} \langle \Lambda_0^{n+1} - \Lambda_0^n, \tilde{\mathcal{L}}^{n+1} \rangle \right] \\ & - \frac{\Delta t}{\alpha} \sum_{n=0}^{M-1} \left[ (\partial_{\Delta t} u^{n+1}, \tilde{\mathcal{L}}^{n+1})_f + \nu_f (\nabla u^{n+\frac{1}{2}}, \nabla \tilde{\mathcal{L}}^{n+1})_f - \langle u^{n+\frac{1}{2}}, \tilde{\mathcal{L}}^{n+1} \rangle \right]. \end{aligned} \qquad (4.28) $$
Note that $\tilde{\Lambda}_0^{n+1} = \Lambda_0^{n+1}$ on $\Sigma$ and $\tilde{\Lambda}_0^{n+1} \in V_f$, so it is indeed an extension. It then follows that

$$ (4.23) \qquad \langle \Lambda_1^{n+1}, \Lambda_0^{n+1} - \Lambda_0^n \rangle = \langle \Lambda_1^{n+1}, \tilde{\Lambda}_0^{n+1} - \tilde{\Lambda}_0^n \rangle. $$

Denote $\mathcal{L}^{n+1} := \Lambda_0^{n+1} - \Lambda_0^n$ and $\tilde{\mathcal{L}}^{n+1} := \tilde{\Lambda}_0^{n+1} - \tilde{\Lambda}_0^n$. Furthermore, owing to (4.4b) with $v = \tilde{\mathcal{L}}^{n+1}$, we have

$$ (4.24) \qquad \langle \Lambda_1^{n+1}, \tilde{\mathcal{L}}^{n+1} \rangle = (\partial_{\Delta t} U_1^{n+1}, \tilde{\mathcal{L}}^{n+1})_f + \nu_f (\nabla U_1^{n+1}, \nabla \tilde{\mathcal{L}}^{n+1})_f - R_2(\tilde{\mathcal{L}}^{n+1}) - T_2(\tilde{\mathcal{L}}^{n+1}). $$

Therefore, we obtain

$$ (4.25) \qquad \frac{\Delta t}{\alpha} \sum_{n=0}^{M-1} \langle \Lambda_1^{n+1}, \Lambda_0^{n+1} - \Lambda_0^n \rangle = \frac{1}{\alpha} \sum_{n=0}^{M-1} (U_1^{n+1} - U_1^n, \tilde{\mathcal{L}}^{n+1})_f + \frac{\Delta t}{\alpha} \sum_{n=0}^{M-1} \nu_f (\nabla U_1^{n+1}, \nabla \tilde{\mathcal{L}}^{n+1})_f \\ - \frac{\Delta t}{\alpha} \sum_{n=0}^{M-1} R_2(\tilde{\mathcal{L}}^{n+1}) - \frac{\Delta t}{\alpha} \sum_{n=0}^{M-1} T_2(\tilde{\mathcal{L}}^{n+1}). $$

The following lemma states the estimate for (4.20), based on the previous identity.

**Lemma 4.6.** Let $1 \le M \le N$, then

$$ (4.26) \qquad \frac{\Delta t}{\alpha} \sum_{n=0}^{M-1} \langle \Lambda_1^{n+1}, \Lambda_0^{n+1} - \Lambda_0^n \rangle \le \frac{1}{8} \sum_{n=0}^{M-1} S^{n+1} + \frac{1}{4} \max_{1 \le n \le M} Z^n + CD_4(M), $$

where

$$ \begin{aligned} D_4(M) := & \Delta t \frac{\nu_f}{\alpha^2} \sum_{n=0}^{M-1} \|\nabla \tilde{\mathcal{L}}^{n+1}\|_{L^2(\Omega_f)}^2 + \frac{\Delta t}{\alpha} \sum_{n=0}^{M-1} \|\tilde{\mathcal{L}}^{n+1}\|_{L^2(\Omega_f)}^2 \\ & + \frac{T(\Delta t)^3}{\alpha^2} \sum_{n=1}^{M-1} \|\frac{\tilde{\mathcal{L}}^{n+1} - \tilde{\mathcal{L}}^n}{(\Delta t)^2}\|_{L^2(\Omega_f)}^2 + \frac{1}{\alpha^2} \|\tilde{\mathcal{L}}^M\|_{L^2(\Omega_f)}^2 \\ & + \Delta t \frac{C_{tr}^2}{\nu_f} \sum_{n=0}^{M-1} \|\Lambda_0^{n+1} - \Lambda_0^n\|_{L^2(\Sigma)}^2 + \Delta t \sum_{n=0}^{M-1} \nu_f \|\nabla(U_0^{n+1} - U_0^n)\|_{L^2(\Omega_f)}^2 \\ & + \frac{\Delta t}{\alpha} \sum_{n=0}^{M-1} \|\partial_{\Delta t} \mathbf{u}^{n+1} - \partial_t \mathbf{u}(t_{n+\frac{1}{2}})\|_{L^2(\Omega_f)}^2 + \Delta t \sum_{n=0}^{M-1} \nu_f \|\nabla \mathbf{u}^{n+\frac{1}{2}} - \nabla \mathbf{u}(t_{n+\frac{1}{2}})\|_{L^2(\Omega_f)}^2 \\ & + \Delta t \sum_{n=0}^{M-1} \frac{C_{tr}^2}{\nu_f} \|\mathbf{u}^{n+\frac{1}{2}} - \mathbf{u}(t_{n+\frac{1}{2}})\|_{L^2(\Sigma)}^2. \end{aligned} \quad (4.27) $$

*Proof.* From (4.25), (4.5b) and (4.2b), it follows that

$$ \begin{aligned} & \frac{\Delta t}{\alpha} \sum_{n=0}^{M-1} \langle \Lambda_1^{n+1}, \Lambda_0^{n+1} - \Lambda_0^n \rangle \\ = & \frac{1}{\alpha} \sum_{n=0}^{M-1} (U_1^{n+1} - U_1^n, \tilde{\mathcal{L}}^{n+1})_f + \frac{\Delta t}{\alpha} \sum_{n=0}^{M-1} \nu_f (\nabla U_1^{n+1}, \nabla \tilde{\mathcal{L}}^{n+1})_f \\ & - \frac{\Delta t}{\alpha} \sum_{n=0}^{M-1} \left[ \frac{1}{2} \nu_f (\nabla (U_0^{n+1} - U_0^n), \nabla \tilde{\mathcal{L}}^{n+1})_f - \frac{1}{2} \langle \Lambda_0^{n+1} - \Lambda_0^n, \tilde{\mathcal{L}}^{n+1} \rangle \right] \\ & - \frac{\Delta t}{\alpha} \sum_{n=0}^{M-1} \left[ (\partial_{\Delta t} \mathbf{u}^{n+1}, \tilde{\mathcal{L}}^{n+1})_f + \nu_f (\nabla \mathbf{u}^{n+\frac{1}{2}}, \nabla \tilde{\mathcal{L}}^{n+1})_f - \langle \mathbf{u}^{n+\frac{1}{2}}, \tilde{\mathcal{L}}^{n+1} \rangle \right]. \end{aligned} \quad (4.28) $$
By using Cauchy-Schwarz inequality, Young's inequality and a similar argument as in Lemma 4.5, we obtain

$$
\begin{align*}
& \frac{\Delta t}{\alpha} \sum_{n=0}^{M-1} (\Lambda_1^{n+1}, \Lambda_0^{n+1} - \Lambda_0^n) \\
\le & \frac{1}{\alpha} \sum_{n=0}^{M-1} (U_1^{n+1} - U_1^n, \tilde{L}^{n+1})_f + \frac{\Delta t}{8} \sum_{n=0}^{M-1} \nu_f \|\nabla U_1^{n+1}\|_{L^2(\Omega_f)}^2 + C \Delta t \frac{\nu_f}{\alpha^2} \sum_{n=0}^{M-1} \|\nabla \tilde{L}^{n+1}\|_{L^2(\Omega_f)}^2 \\
& + \Delta t \sum_{n=0}^{M-1} \nu_f \|\nabla (U_0^{n+1} - U_0^n)\|_{L^2(\Omega_f)}^2 + \Delta t \frac{C_{tr}^2}{\nu_f} \sum_{n=0}^{M-1} \|\Lambda_0^{n+1} - \Lambda_0^n\|_{L^2(\Sigma)}^2 \\
& + C \frac{\Delta t}{\alpha} \sum_{n=0}^{M-1} \|\tilde{L}^{n+1}\|_{L^2(\Omega_f)}^2 + C \frac{\Delta t}{\alpha} \sum_{n=0}^{M-1} \|\partial_{\Delta t} \mathbf{u}^{n+1} - \partial_t \mathbf{u}(t_{n+\frac{1}{2}})\|_{L^2(\Omega_f)}^2 \\
& + C \Delta t \sum_{n=0}^{M-1} \nu_f \|\nabla \mathbf{u}^{n+\frac{1}{2}} - \nabla \mathbf{u}(t_{n+\frac{1}{2}})\|_{L^2(\Omega_f)}^2 + C \Delta t \sum_{n=0}^{M-1} \frac{C_{tr}^2}{\nu_f} \|l^{n+\frac{1}{2}} - l(t_{n+\frac{1}{2}})\|_{L^2(\Sigma)}^2 \\
\le & \frac{1}{\alpha} \sum_{n=0}^{M-1} (U_1^{n+1} - U_1^n, \tilde{L}^{n+1})_f + \frac{1}{8} \sum_{n=0}^{M-1} S_1^{n+1} + CD_4(M).
\end{align*}
$$

By a summation by parts formula, and using $U_1^0 = 0$, we have

$$
\begin{align*}
\frac{1}{\alpha} \sum_{n=0}^{M-1} (U_1^{n+1} - U_1^n, \tilde{L}^{n+1})_f &= \frac{1}{\alpha} \sum_{n=1}^{M-1} (U_1^n, \tilde{L}^n - \tilde{L}^{n+1})_f + \frac{1}{\alpha} (U_1^M, \tilde{L}^M)_f - \frac{1}{\alpha} (U_1^0, \tilde{L}^1)_f \\
&= -\frac{1}{\alpha} \sum_{n=1}^{M-1} (U_1^n, \tilde{L}^{n+1} - \tilde{L}^n)_f + \frac{1}{\alpha} (U_1^M, \tilde{L}^M)_f \\
&\le \frac{\Delta t}{16T} \sum_{n=1}^{M-1} \|U_1^n\|_{L^2(\Omega_f)}^2 + \frac{CT(\Delta t)^3}{\alpha^2} \sum_{n=1}^{M-1} \left\| \frac{\tilde{L}^{n+1} - \tilde{L}^n}{(\Delta t)^2} \right\|_{L^2(\Omega_f)}^2 \\
&\quad + \frac{1}{16} \|U_1^M\|_{L^2(\Omega_f)}^2 + \frac{C}{\alpha^2} \|\tilde{L}^M\|_{L^2(\Omega_f)}^2 \le \frac{1}{4} \max_{1 \le n \le M} \mathbf{Z}^n + CD_4(M).
\end{align*}
$$

The estimate (4.26) then follows by combining the above two inequalities, which completes the proof. $\square$

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

**4.7. Main result: error estimate.** In this section, we state and prove our main result under Assumption 3.1 and Corollary 3.2. To this purpose, we first need the following lemma.

**Lemma 4.7.** We have, for $1 \le M \le N$,

$$
(4.29) \quad D_1(M) \le C(\Delta t)^4 \left( \frac{C_{tr}^2}{\nu_f} + \frac{C_{tr}^2}{\nu_s} \right) Y,
$$

$$
(4.30) \quad D_2(M) \le C(\Delta t)^4 \left( \frac{C_{tr}^2}{\nu_f} + \frac{C_{tr}^2}{\nu_s} \right) Y + \left( \frac{C_{tr}^4}{\nu_s^2} \alpha^2 + 1 \right) y,
$$

$$
(4.31) \quad D_3(M) \le (\Delta t)^4 \left( \begin{aligned}[t]
& T \|\partial_t^3 \mathbf{w}\|_{L^2((0,T),L^2(\Omega_s))}^2 + T \|\partial_t^3 \mathbf{u}\|_{L^2((0,T),L^2(\Omega_f))}^2 \\
& + \nu_s \|\partial_t^2 \mathbf{w}\|_{L^2((0,T),H^1(\Omega_s))}^2 + \nu_f \|\partial_t^2 \mathbf{u}\|_{L^2((0,T),H^1(\Omega_f))}^2 \\
& + \left( \frac{C_{tr}^2}{\nu_f} + \frac{C_{tr}^2}{\nu_s} \right) \|\partial_t^2 \mathbf{l}\|_{L^2((0,T),L^2(\Sigma))}^2
\end{aligned} \right),
$$
By using Cauchy-Schwarz inequality, Young's inequality and a similar argument as in Lemma 4.5, we obtain

$$
\begin{align*}
& \frac{\Delta t}{\alpha} \sum_{n=0}^{M-1} (\Lambda_1^{n+1}, \Lambda_0^{n+1} - \Lambda_0^n) \\
\le & \frac{1}{\alpha} \sum_{n=0}^{M-1} (U_1^{n+1} - U_1^n, \tilde{L}^{n+1})_f + \frac{\Delta t}{8} \sum_{n=0}^{M-1} \nu_f \|\nabla U_1^{n+1}\|_{L^2(\Omega_f)}^2 + C \Delta t \frac{\nu_f}{\alpha^2} \sum_{n=0}^{M-1} \|\nabla \tilde{L}^{n+1}\|_{L^2(\Omega_f)}^2 \\
& + \Delta t \sum_{n=0}^{M-1} \nu_f \|\nabla (U_0^{n+1} - U_0^n)\|_{L^2(\Omega_f)}^2 + \Delta t \frac{C_{tr}^2}{\nu_f} \sum_{n=0}^{M-1} \|\Lambda_0^{n+1} - \Lambda_0^n\|_{L^2(\Sigma)}^2 \\
& + C \frac{\Delta t}{\alpha} \sum_{n=0}^{M-1} \|\tilde{L}^{n+1}\|_{L^2(\Omega_f)}^2 + C \frac{\Delta t}{\alpha} \sum_{n=0}^{M-1} \|\partial_{\Delta t} u^{n+1} - \partial_t u(t_{n+\frac{1}{2}})\|_{L^2(\Omega_f)}^2 \\
& + C \Delta t \sum_{n=0}^{M-1} \nu_f \|\nabla u^{n+\frac{1}{2}} - \nabla u(t_{n+\frac{1}{2}})\|_{L^2(\Omega_f)}^2 + C \Delta t \sum_{n=0}^{M-1} \frac{C_{tr}^2}{\nu_f} \|l^{n+\frac{1}{2}} - l(t_{n+\frac{1}{2}})\|_{L^2(\Sigma)}^2 \\
\le & \frac{1}{\alpha} \sum_{n=0}^{M-1} (U_1^{n+1} - U_1^n, \tilde{L}^{n+1})_f + \frac{1}{8} \sum_{n=0}^{M-1} S_1^{n+1} + CD_4(M).
\end{align*}
$$

By a summation by parts formula, and using $U_1^0 = 0$, we have

$$
\begin{align*}
\frac{1}{\alpha} \sum_{n=0}^{M-1} (U_1^{n+1} - U_1^n, \tilde{L}^{n+1})_f &= \frac{1}{\alpha} \sum_{n=1}^{M-1} (U_1^n, \tilde{L}^n - \tilde{L}^{n+1})_f + \frac{1}{\alpha} (U_1^M, \tilde{L}^M)_f - \frac{1}{\alpha} (U_1^0, \tilde{L}^1)_f \\
&= -\frac{1}{\alpha} \sum_{n=1}^{M-1} (U_1^n, \tilde{L}^{n+1} - \tilde{L}^n)_f + \frac{1}{\alpha} (U_1^M, \tilde{L}^M)_f \\
&\le \frac{\Delta t}{16T} \sum_{n=1}^{M-1} \|U_1^n\|_{L^2(\Omega_f)}^2 + \frac{CT(\Delta t)^3}{\alpha^2} \sum_{n=1}^{M-1} \left\| \frac{\tilde{L}^{n+1} - \tilde{L}^n}{(\Delta t)^2} \right\|_{L^2(\Omega_f)}^2 \\
&\quad + \frac{1}{16} \|U_1^M\|_{L^2(\Omega_f)}^2 + \frac{C}{\alpha^2} \|\tilde{L}^M\|_{L^2(\Omega_f)}^2 \le \frac{1}{4} \max_{1 \le n \le M} Z^n + CD_4(M).
\end{align*}
$$

The estimate (4.26) then follows by combining the above two inequalities, which completes the proof. $\square$

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

4.7. Main result: error estimate. In this section, we state and prove our main result under Assumption 3.1 and Corollary 3.2. To this purpose, we first need the following lemma.

**Lemma 4.7.** We have, for $1 \le M \le N$,

$$
(4.29) \quad D_1(M) \le C(\Delta t)^4 \left( \frac{C_{tr}^2}{\nu_f} + \frac{C_{tr}^2}{\nu_s} \right) Y,
$$

$$
(4.30) \quad D_2(M) \le C(\Delta t)^4 \left( \frac{C_{tr}^2}{\nu_f} + \frac{C_{tr}^2}{\nu_s} \right) Y + \left( \frac{C_{tr}^4}{\nu_s^2} \alpha^2 + 1 \right) y,
$$

$$
(4.31) \quad D_3(M) \le (\Delta t)^4 \left( \begin{aligned}[t]
& T \|\partial_t^3 \mathbf{w}\|_{L^2((0,T),L^2(\Omega_s))}^2 + T \|\partial_t^3 \mathbf{u}\|_{L^2((0,T),L^2(\Omega_f))}^2 \\
& + \nu_s \|\partial_t^2 \mathbf{w}\|_{L^2((0,T),H^1(\Omega_s))}^2 + \nu_f \|\partial_t^2 \mathbf{u}\|_{L^2((0,T),H^1(\Omega_f))}^2 \\
& + \left( \frac{C_{tr}^2}{\nu_f} + \frac{C_{tr}^2}{\nu_s} \right) \|\partial_t^2 \mathbf{l}\|_{L^2((0,T),L^2(\Sigma))}^2
\end{aligned} \right),
$$
$$ (4.32) \quad D_4(M) \le C(\Delta t)^4 \left( \left( \frac{\nu_f}{\alpha^2} + \frac{1}{\alpha} + 1 \right) \mathcal{Y} + \left( \frac{\nu_f^2}{\alpha^2} + \frac{\nu_f}{\alpha^2} + \frac{C_{tr}^2}{\nu_f} \right) \mathcal{Y} + \frac{T}{\alpha^2} \mathcal{Y} \right. \\ \left. + \frac{1}{\alpha} \|\partial_t^3 \mathbf{u}\|_{L^2((0,T),L^2(\Omega_s))}^2 + \nu_f \|\partial_t^2 \mathbf{u}\|_{L^2((0,T),H^1(\Omega_s))}^2 + \frac{C_{tr}^2}{\nu_f} \|\partial_t^2 \mathbf{I}\|_{L^2((0,T),L^2(\Sigma))}^2 \right), $$

where $\mathcal{Y}$ is defined in Remark 3.3 and $\mathcal{Y}$, $\mathcal{Z}$ are defined in Appendix A.

*Proof.* For $D_1(M)$, using Corollary 3.2, we have

$$ (4.33) \quad D_1(M) = \Delta t \sum_{n=0}^{M-1} \left( \frac{C_{tr}^2}{\nu_f} + \frac{C_{tr}^2}{\nu_s} \right) \| \Lambda_0^n - \Lambda_0^{n+1} \|_{L^2(\Sigma)}^2 \le \left( \frac{C_{tr}^2}{\nu_f} + \frac{C_{tr}^2}{\nu_s} \right) (\Delta t)^4 \mathcal{Y}. $$

From (3.2) and Corollary 3.2, it follows that

$$ \begin{aligned} D_2(M) = & \Delta t \sum_{n=0}^{M-1} \left( \frac{C_{tr}^2}{\nu_f} + \frac{C_{tr}^2}{\nu_s} \right) \| \Lambda_0^n - \Lambda_0^{n+1} \|_{L^2(\Sigma)}^2 + \Delta t \sum_{n=0}^{M-1} \nu_s \| \nabla (W_0^{n+1} - W_0^n) \|_{L^2(\Omega_s)}^2 \\ & + \Delta t \sum_{n=0}^{M-1} \frac{C_{tr}^2}{\nu_s} \alpha^2 \| W_0^{n+1} - W_0^n \|_{L^2(\Sigma)}^2 + \Delta t \sum_{n=0}^{M-1} \nu_f \| \nabla (U_0^{n+1} - U_0^n) \|_{L^2(\Omega_f)}^2 \\ & \le \left( \frac{C_{tr}^2}{\nu_f} + \frac{C_{tr}^2}{\nu_s} \right) (\Delta t)^4 \mathcal{Y} + C (\Delta t)^4 \mathcal{Y} + C \frac{C_{tr}^4}{\nu_s^2} \alpha^2 (\Delta t)^4 \mathcal{Y} + C (\Delta t)^4 \mathcal{Y}. \end{aligned} $$

For $D_3(M)$, we first notice that (see, e.g., [14] for similar bounds)

$$ (4.34a) \quad \sum_{n=0}^{M-1} \| \partial_t \mathbf{w}^{n+1} - \partial_t \mathbf{w}(t_{n+\frac{1}{2}}) \|_{L^2(\Omega_s)}^2 \le \Delta t^3 \| \partial_t^3 \mathbf{w} \|_{L^2((0,T),L^2(\Omega_s))}^2, $$

$$ (4.34b) \quad \sum_{n=0}^{M-1} \| \nabla (\mathbf{w}^{n+1} - \mathbf{w}(t_{n+\frac{1}{2}})) \|_{L^2(\Omega_s)}^2 \le \Delta t^3 \| \partial_t^2 \mathbf{w} \|_{L^2((0,T),H^1(\Omega_s))}^2, $$

$$ (4.34c) \quad \sum_{n=0}^{M-1} \| \mathbf{I}^{n+1} - \mathbf{I}(t_{n+\frac{1}{2}}) \|_{L^2(\Sigma)}^2 \le \Delta t^3 \| \partial_t^2 \mathbf{I} \|_{L^2((0,T),L^2(\Sigma))}^2. $$

Hence, from these bounds, we get

$$ \begin{aligned} D_3(M) = \Delta t T \sum_{n=0}^{M-1} & \left( \| \partial_{\Delta t} \mathbf{w}^{n+1} - \partial_t \mathbf{w}(t_{n+\frac{1}{2}}) \|_{L^2(\Omega_s)}^2 + \| \partial_{\Delta t} \mathbf{u}^{n+1} - \partial_t \mathbf{u}(t_{n+\frac{1}{2}}) \|_{L^2(\Omega_f)}^2 \right) \\ & + \Delta t \sum_{n=0}^{M-1} \left( \nu_s \| \nabla \mathbf{w}^{n+\frac{1}{2}} - \nabla \mathbf{w}(t_{n+\frac{1}{2}}) \|_{L^2(\Omega_s)}^2 + \nu_f \| \nabla \mathbf{u}^{n+\frac{1}{2}} - \nabla \mathbf{u}(t_{n+\frac{1}{2}}) \|_{L^2(\Omega_f)}^2 \right) \\ & + \Delta t \sum_{n=0}^{M-1} \left( \left( \frac{C_{tr}^2}{\nu_f} + \frac{C_{tr}^2}{\nu_s} \right) \| \mathbf{I}^{n+\frac{1}{2}} - \mathbf{I}(t_{n+\frac{1}{2}}) \|_{L^2(\Sigma)}^2 \right) \\ & \le T(\Delta t)^4 ( \| \partial_t^3 \mathbf{w} \|_{L^2((0,T),L^2(\Omega_s))}^2 + \| \partial_t^3 \mathbf{u} \|_{L^2((0,T),L^2(\Omega_f))}^2 ) \\ & + (\Delta t)^4 ( \nu_s \| \partial_t^2 \mathbf{w} \|_{L^2((0,T),H^1(\Omega_s))}^2 + \nu_f \| \partial_t^2 \mathbf{u} \|_{L^2((0,T),H^1(\Omega_f))}^2 ) \\ & + (\Delta t)^4 \left( \frac{C_{tr}^2}{\nu_f} + \frac{C_{tr}^2}{\nu_s} \right) \| \partial_t^2 \mathbf{I} \|_{L^2((0,T),L^2(\Sigma))}^2. \end{aligned} $$
$$ (4.32) \quad D_4(M) \le C(\Delta t)^4 \left( \left( \frac{\nu_f}{\alpha^2} + \frac{1}{\alpha} + 1 \right) Y + \left( \frac{\nu_f^2}{\alpha^2} + \frac{\nu_f}{\alpha^2} + \frac{C_{tr}^2}{\nu_f} \right) Y + \frac{T}{\alpha^2} \mathfrak{Z} \right. \\ \left. + \frac{1}{\alpha} \|\partial_t^3 u\|_{L^2((0,T),L^2(\Omega_s))}^2 + \nu_f \|\partial_t^2 u\|_{L^2((0,T),H^1(\Omega_s))}^2 + \frac{C_{tr}^2}{\nu_f} \|\partial_t^2 l\|_{L^2((0,T),L^2(\Sigma))}^2 \right), $$

where $Y$ is defined in Remark 3.3 and $Y$, $\mathfrak{Z}$ are defined in Appendix A.

*Proof.* For $D_1(M)$, using Corollary 3.2, we have

$$ (4.33) \quad D_1(M) = \Delta t \sum_{n=0}^{M-1} \left( \frac{C_{tr}^2}{\nu_f} + \frac{C_{tr}^2}{\nu_s} \right) \|\Lambda_0^n - \Lambda_0^{n+1}\|_{L^2(\Sigma)}^2 \le \left( \frac{C_{tr}^2}{\nu_f} + \frac{C_{tr}^2}{\nu_s} \right) (\Delta t)^4 Y. $$

From (3.2) and Corollary 3.2, it follows that

$$ \begin{aligned} D_2(M) ={}& \Delta t \sum_{n=0}^{M-1} \left( \frac{C_{tr}^2}{\nu_f} + \frac{C_{tr}^2}{\nu_s} \right) \|\Lambda_0^n - \Lambda_0^{n+1}\|_{L^2(\Sigma)}^2 + \Delta t \sum_{n=0}^{M-1} \nu_s \|\nabla(W_0^{n+1} - W_0^n)\|_{L^2(\Omega_s)}^2 \\ & + \Delta t \sum_{n=0}^{M-1} \frac{C_{tr}^2}{\nu_s} \alpha^2 \|W_0^{n+1} - W_0^n\|_{L^2(\Sigma)}^2 + \Delta t \sum_{n=0}^{M-1} \nu_f \|\nabla(U_0^{n+1} - U_0^n)\|_{L^2(\Omega_f)}^2 \\ & \le \left( \frac{C_{tr}^2}{\nu_f} + \frac{C_{tr}^2}{\nu_s} \right) (\Delta t)^4 Y + C(\Delta t)^4 Y + C \frac{C_{tr}^4}{\nu_s^2} \alpha^2 (\Delta t)^4 Y + C(\Delta t)^4 Y. \end{aligned} $$

For $D_3(M)$, we first notice that (see, e.g., [14] for similar bounds)

$$ (4.34a) \quad \sum_{n=0}^{M-1} \|\partial_t w^{n+1} - \partial_t w(t_{n+\frac{1}{2}})\|_{L^2(\Omega_s)}^2 \le \Delta t^3 \|\partial_t^3 w\|_{L^2((0,T),L^2(\Omega_s))}^2, $$

$$ (4.34b) \quad \sum_{n=0}^{M-1} \|\nabla(w^{n+1} - w(t_{n+\frac{1}{2}}))\|_{L^2(\Omega_s)}^2 \le \Delta t^3 \|\partial_t^2 w\|_{L^2((0,T),H^1(\Omega_s))}^2, $$

$$ (4.34c) \quad \sum_{n=0}^{M-1} \|l^{n+1} - l(t_{n+\frac{1}{2}})\|_{L^2(\Sigma)}^2 \le \Delta t^3 \|\partial_t^2 l\|_{L^2((0,T),L^2(\Sigma))}^2. $$

Hence, from these bounds, we get

$$ \begin{aligned} D_3(M) ={}& \Delta t T \sum_{n=0}^{M-1} \left( \|\partial_{\Delta t} w^{n+1} - \partial_t w(t_{n+\frac{1}{2}})\|_{L^2(\Omega_s)}^2 + \|\partial_{\Delta t} u^{n+1} - \partial_t u(t_{n+\frac{1}{2}})\|_{L^2(\Omega_f)}^2 \right) \\ & + \Delta t \sum_{n=0}^{M-1} \left( \nu_s \|\nabla w^{n+\frac{1}{2}} - \nabla w(t_{n+\frac{1}{2}})\|_{L^2(\Omega_s)}^2 + \nu_f \|\nabla u^{n+\frac{1}{2}} - \nabla u(t_{n+\frac{1}{2}})\|_{L^2(\Omega_f)}^2 \right) \\ & + \Delta t \sum_{n=0}^{M-1} \left( \left( \frac{C_{tr}^2}{\nu_f} + \frac{C_{tr}^2}{\nu_s} \right) \|l^{n+\frac{1}{2}} - l(t_{n+\frac{1}{2}})\|_{L^2(\Sigma)}^2 \right) \\ & \le T(\Delta t)^4 (\|\partial_t^3 w\|_{L^2((0,T),L^2(\Omega_s))}^2 + \|\partial_t^3 u\|_{L^2((0,T),L^2(\Omega_f))}^2) \\ & + (\Delta t)^4 (\nu_s \|\partial_t^2 w\|_{L^2((0,T),H^1(\Omega_s))}^2 + \nu_f \|\partial_t^2 u\|_{L^2((0,T),H^1(\Omega_f))}^2) \\ & + (\Delta t)^4 \left( \frac{C_{tr}^2}{\nu_f} + \frac{C_{tr}^2}{\nu_s} \right) \|\partial_t^2 l\|_{L^2((0,T),L^2(\Sigma))}^2. \end{aligned} $$
For $D_4(M)$, we first estimate the terms involving $\tilde{L}^{n+1}$. From (4.21b), (4.21d), (3.2) and Assumption 3.1, it follows that

$$ (4.36) \qquad \begin{aligned} & \Delta t \frac{\nu_f}{\alpha^2} \sum_{n=0}^{M-1} \|\nabla \tilde{L}^{n+1}\|_{L^2(\Omega_f)}^2 \\ & \le \Delta t \frac{\nu_f}{\alpha^2} \sum_{n=0}^{M-1} \nu_f \|\nabla (U_0^{n+1} - U_0^n)\|_{L^2(\Omega_f)}^2 + \Delta t \frac{\nu_f}{\alpha^2} \sum_{n=0}^{M-1} \nu_f \|D^2(U_0^{n+1} - U_0^n)\|_{L^2(\Omega_f)}^2 \\ & \le C(\Delta t)^4 \left( \frac{\nu_f}{\alpha^2} y + \frac{\nu_f^2}{\alpha^2} y \right). \end{aligned} $$

Similarly, from (3.2), we have

$$ (4.37) \qquad \frac{\Delta t}{\alpha} \sum_{n=0}^{M-1} \|\tilde{L}^{n+1}\|_{L^2(\Omega_f)}^2 \le C \frac{\Delta t}{\alpha} \sum_{n=0}^{M-1} \nu_f \|\nabla(U_0^{n+1} - U_0^n)\|_{L^2(\Omega_f)}^2 \le \frac{C}{\alpha} (\Delta t)^4 y. $$

On the other hand, from (3.3), we have

$$ (4.38) \qquad \begin{aligned} & \frac{T(\Delta t)^3}{\alpha^2} \sum_{n=1}^{M-1} \|\frac{\tilde{L}^{n+1} - \tilde{L}^n}{(\Delta t)^2}\|_{L^2(\Omega_f)}^2 \le C \frac{T(\Delta t)^3}{\alpha^2} \sum_{n=1}^{M-1} \nu_f \left\| \frac{\nabla(U_0^{n+1} - 2U_0^n + U_0^{n-1})}{(\Delta t)^2} \right\|_{L^2(\Omega_f)}^2 \\ & \le C \frac{T}{\alpha^2} (\Delta t)^4 y. \end{aligned} $$

At last, from Assumption 3.1, we get

$$ (4.39) \qquad \frac{1}{\alpha^2} \|\tilde{L}^M\|_{L^2(\Omega_f)}^2 \le C \frac{\nu_f}{\alpha^2} \|\nabla(U_0^{M+1} - U_0^M)\|_{L^2(\Omega_f)}^2 \le C \frac{\nu_f}{\alpha^2} (\Delta t)^4 y. $$

The remaining terms in $D_4(M)$ can be estimated similarly to those of $D_1(M)$ to $D_3(M)$ as follows, using (3.2), (3.6) and (4.34a)-(4.34c), yielding

$$ (4.40) \qquad \begin{aligned} & \Delta t \frac{C_{tr}^2}{\nu_f} \sum_{n=0}^{M-1} \|\Lambda_0^{n+1} - \Lambda_0^n\|_{L^2(\Sigma)}^2 + \Delta t \sum_{n=0}^{M-1} \nu_f \|\nabla(U_0^{n+1} - U_0^n)\|_{L^2(\Omega_f)}^2 \\ & + \frac{\Delta t}{\alpha} \sum_{n=0}^{M-1} \|\partial_{\Delta t} u^{n+1} - \partial_t u(t_{n+\frac{1}{2}})\|_{L^2(\Omega_f)}^2 + \Delta t \sum_{n=0}^{M-1} \nu_f \|\nabla u^{n+\frac{1}{2}} - \nabla u(t_{n+\frac{1}{2}})\|_{L^2(\Omega_f)}^2 \\ & + \Delta t \sum_{n=0}^{M-1} \frac{C_{tr}^2}{\nu_f} \|l^{n+\frac{1}{2}} - l(t_{n+\frac{1}{2}})\|_{L^2(\Sigma)}^2 \\ & \le \frac{C_{tr}^2}{\nu_f} (\Delta t)^4 y + C(\Delta t)^4 y + \frac{1}{\alpha} \Delta t^4 \|\partial_t^3 u\|_{L^2((0,T),L^2(\Omega_s))}^2 \\ & + \nu_f \Delta t^4 \|\partial_t^2 u\|_{L^2((0,T),H^1(\Omega_s))}^2 + \frac{C_{tr}^2}{\nu_f} \Delta t^4 \|\partial_t^2 l\|_{L^2((0,T),L^2(\Sigma))}^2. \end{aligned} $$

We conclude the proof by combining all the above inequalities. $\square$

We now state and prove the main result of this paper.

**Theorem 4.8.** The following error estimate holds

$$ (4.41) \qquad \max_{1 \le M \le N} Z_1^M + \sum_{n=0}^{N-1} S_1^{n+1} \le C(\Delta t)^4 Y_1, $$
For $D_4(M)$, we first estimate the terms involving $\tilde{L}^{n+1}$. From (4.21b), (4.21d), (3.2) and Assumption 3.1, it follows that

$$ (4.36) \qquad \begin{aligned} & \Delta t \frac{\nu_f}{\alpha^2} \sum_{n=0}^{M-1} \|\nabla \tilde{L}^{n+1}\|_{L^2(\Omega_f)}^2 \\ & \le \Delta t \frac{\nu_f}{\alpha^2} \sum_{n=0}^{M-1} \nu_f \|\nabla (U_0^{n+1} - U_0^n)\|_{L^2(\Omega_f)}^2 + \Delta t \frac{\nu_f}{\alpha^2} \sum_{n=0}^{M-1} \nu_f \|D^2(U_0^{n+1} - U_0^n)\|_{L^2(\Omega_f)}^2 \\ & \le C(\Delta t)^4 \left( \frac{\nu_f}{\alpha^2} y + \frac{\nu_f^2}{\alpha^2} y \right). \end{aligned} $$

Similarly, from (3.2), we have

$$ (4.37) \qquad \frac{\Delta t}{\alpha} \sum_{n=0}^{M-1} \|\tilde{L}^{n+1}\|_{L^2(\Omega_f)}^2 \le C \frac{\Delta t}{\alpha} \sum_{n=0}^{M-1} \nu_f \|\nabla(U_0^{n+1} - U_0^n)\|_{L^2(\Omega_f)}^2 \le \frac{C}{\alpha} (\Delta t)^4 y. $$

On the other hand, from (3.3), we have

$$ (4.38) \qquad \begin{aligned} & \frac{T(\Delta t)^3}{\alpha^2} \sum_{n=1}^{M-1} \|\frac{\tilde{L}^{n+1} - \tilde{L}^n}{(\Delta t)^2}\|_{L^2(\Omega_f)}^2 \le C \frac{T(\Delta t)^3}{\alpha^2} \sum_{n=1}^{M-1} \nu_f \left\| \frac{\nabla(U_0^{n+1} - 2U_0^n + U_0^{n-1})}{(\Delta t)^2} \right\|_{L^2(\Omega_f)}^2 \\ & \le C \frac{T}{\alpha^2} (\Delta t)^4 y. \end{aligned} $$

At last, from Assumption 3.1, we get

$$ (4.39) \qquad \frac{1}{\alpha^2} \|\tilde{L}^M\|_{L^2(\Omega_f)}^2 \le C \frac{\nu_f}{\alpha^2} \|\nabla(U_0^{M+1} - U_0^M)\|_{L^2(\Omega_f)}^2 \le C \frac{\nu_f}{\alpha^2} (\Delta t)^4 y. $$

The remaining terms in $D_4(M)$ can be estimated similarly to those of $D_1(M)$ to $D_3(M)$ as follows, using (3.2), (3.6) and (4.34a)-(4.34c), yielding

$$ (4.40) \qquad \begin{aligned} & \Delta t \frac{C_{tr}^2}{\nu_f} \sum_{n=0}^{M-1} \|\Lambda_0^{n+1} - \Lambda_0^n\|_{L^2(\Sigma)}^2 + \Delta t \sum_{n=0}^{M-1} \nu_f \|\nabla(U_0^{n+1} - U_0^n)\|_{L^2(\Omega_f)}^2 \\ & \qquad + \frac{\Delta t}{\alpha} \sum_{n=0}^{M-1} \|\partial_{\Delta t} u^{n+1} - \partial_t u(t_{n+\frac{1}{2}})\|_{L^2(\Omega_f)}^2 + \Delta t \sum_{n=0}^{M-1} \nu_f \|\nabla u^{n+\frac{1}{2}} - \nabla u(t_{n+\frac{1}{2}})\|_{L^2(\Omega_f)}^2 \\ & \qquad + \Delta t \sum_{n=0}^{M-1} \frac{C_{tr}^2}{\nu_f} \|u^{n+\frac{1}{2}} - l(t_{n+\frac{1}{2}})\|_{L^2(\Sigma)}^2 \\ & \le \frac{C_{tr}^2}{\nu_f} (\Delta t)^4 y + C(\Delta t)^4 y + \frac{1}{\alpha} \Delta t^4 \|\partial_t^3 u\|_{L^2((0,T),L^2(\Omega_s))}^2 \\ & \qquad + \nu_f \Delta t^4 \|\partial_t^2 u\|_{L^2((0,T),H^1(\Omega_s))}^2 + \frac{C_{tr}^2}{\nu_f} \Delta t^4 \|\partial_t^2 l\|_{L^2((0,T),L^2(\Sigma))}^2. \end{aligned} $$

We conclude the proof by combining all the above inequalities. □

We now state and prove the main result of this paper.

**Theorem 4.8.** The following error estimate holds

$$ (4.41) \qquad \max_{1 \le M \le N} Z_1^M + \sum_{n=0}^{N-1} S_1^{n+1} \le C(\Delta t)^4 Y_1, $$
where

$$
\begin{aligned}
Y_1 := & \left( \frac{C_{tr}^2}{\nu_f} + \frac{C_{tr}^2}{\nu_s} + \frac{C_{tr}^2}{\nu_f^2} \alpha^2 + \frac{\nu_f + \alpha}{\alpha^2} + 1 \right) Y + \left( \frac{C_{tr}^2}{\nu_f} + \frac{C_{tr}^2}{\nu_s} + \frac{\nu_f + \nu_f^2}{\alpha^2} \right) Y \\
& + \frac{T}{\alpha^2} \mathfrak{Z} + T \| \partial_t^3 \mathbf{w} \|_{L^2((0,T),L^2(\Omega_s))}^2 + \left( T + \frac{1}{\alpha} \right) \| \partial_t^3 \mathbf{u} \|_{L^2((0,T),L^2(\Omega_f))}^2 \\
& + \nu_s \| \partial_t^2 \mathbf{w} \|_{L^2((0,T),H^1(\Omega_s))}^2 + \nu_f \| \partial_t^2 \mathbf{u} \|_{L^2((0,T),H^1(\Omega_f))}^2 + \left( \frac{C_{tr}^2}{\nu_f} + \frac{C_{tr}^2}{\nu_s} \right) \| \partial_t^2 \mathbf{l} \|_{L^2((0,T),L^2(\Sigma))}^2.
\end{aligned}
$$

Here, the notation $Y$ is defined in Assumption 3.1 and $Y, \mathfrak{Z}$ are defined in Appendix A.

*Proof.* For $1 \le M \le N$, by taking the sum over $n = 0, \dots, M-1$ in Lemma 4.1 and using Lemmas 4.3 and 4.5, we have

$$Z_1^M + \frac{1}{2} \sum_{n=0}^{M-1} S_1^{n+1} \le \frac{1}{2} \max_{1 \le n \le M} Z_1^n + C \sum_{i=1}^4 D_i(M).$$

Taking the maximum over $1 \le M \le N$, yields

$$\frac{1}{2} \max_{1 \le M \le N} Z^M \le C \max_{1 \le M \le N} \sum_{i=1}^4 D_i(M).$$

Hence, we also get

$$\sum_{n=0}^{N-1} S^{n+1} \le C \max_{1 \le M \le N} \sum_{i=1}^{4} D_i(M).$$

The estimate (4.41) then follows from Lemma 4.7, which concludes the proof. $\square$

## 5. NUMERICAL EXPERIMENTS

In this section, we present two examples which numerically illustrate the analysis and observations of the previous sections. We use the following notations to denote errors of the correction step,

$$e_{u1} := \|u_1^N - u^N\|_{L^2(\Omega_f)},$$

$$e_{w1} := \|w_1^N - w^N\|_{L^2(\Omega_s)},$$

$$e_{du1} := \|\nabla(u_1^N - u^N)\|_{L^2(\Omega_f)},$$

$$e_{\lambda} := \|\lambda_0^N - I^N\|_{L^2(\Sigma)},$$

$$e_{1,\lambda} := \|(\lambda_0^N - I^N) - (\lambda_0^{N-1} - I^{N-1})\|_{L^2(\Sigma)}.$$

All numerical experiments are performed using FEniCS and multiphenics [2, 4].

*Example 5.1 (Non-horizontal interface).* In this example, we test the defect-correction scheme (2.7) in a problem with non-horizontal interface. We let $\Omega = (0, 1)^2$ and we take the solution to be

$$\mathbf{w} = \mathbf{u} = e^{-2\pi^2 t} \cos(\pi x_1) \sin(\pi x_2).$$

We let $\Sigma$ be defined as the straight line connecting $(0, 0.25)$ and $(1, 0.75)$. We then define $\Omega_s$ as the region above $\Sigma$ and $\Omega_f$ as the region below $\Sigma$ (see Figure 1 for an illustration). The diffusion parameters are both set to be one, $\nu_s = \nu_f = 1$. We choose $h = \Delta t, T = 0.25, \alpha = 4$ and use the piece-wise linear finite element method for the spatial discretization. The results are reported in Table 1. As we can observe, second-order convergence $\mathcal{O}((\Delta t)^2)$ is obtained for $e_{u1}$ and $e_{w1}$, which is in agreement with the theoretical findings provided by Theorem 4.8. Finally, we also observe that $e_{1,\lambda}$ converges as $\mathcal{O}((\Delta t)^2)$. This matches with Corollary 3.2.
where

$$Y_1 := \left( \frac{C_{tr}^2}{\nu_f} + \frac{C_{tr}^2}{\nu_s} + \frac{C_{tr}^2}{\nu_f^2} \alpha^2 + \frac{\nu_f + \alpha}{\alpha^2} + 1 \right) Y + \left( \frac{C_{tr}^2}{\nu_f} + \frac{C_{tr}^2}{\nu_s} + \frac{\nu_f + \nu_f^2}{\alpha^2} \right) Y \\ + \frac{T}{\alpha^2} \mathfrak{Z} + T \| \partial_t^3 w \|_{L^2((0,T),L^2(\Omega_s))}^2 + \left( T + \frac{1}{\alpha} \right) \| \partial_t^3 u \|_{L^2((0,T),L^2(\Omega_f))}^2 \\ + \nu_s \| \partial_t^2 w \|_{L^2((0,T),H^1(\Omega_s))}^2 + \nu_f \| \partial_t^2 u \|_{L^2((0,T),H^1(\Omega_f))}^2 + \left( \frac{C_{tr}^2}{\nu_f} + \frac{C_{tr}^2}{\nu_s} \right) \| \partial_t^2 l \|_{L^2((0,T),L^2(\Sigma))}^2.$$

Here, the notation $Y$ is defined in Assumption 3.1 and $Y, \mathfrak{Z}$ are defined in Appendix A.

*Proof.* For $1 \le M \le N$, by taking the sum over $n = 0, \dots, M-1$ in Lemma 4.1 and using Lemmas 4.3 and 4.5, we have

$$Z_1^M + \frac{1}{2} \sum_{n=0}^{M-1} S_1^{n+1} \le \frac{1}{2} \max_{1 \le n \le M} Z_1^n + C \sum_{i=1}^4 D_i(M).$$

Taking the maximum over $1 \le M \le N$, yields

$$\frac{1}{2} \max_{1 \le M \le N} Z^M \le C \max_{1 \le M \le N} \sum_{i=1}^4 D_i(M).$$

Hence, we also get

$$\sum_{n=0}^{N-1} S^{n+1} \le C \max_{1 \le M \le N} \sum_{i=1}^{4} D_i(M).$$

The estimate (4.41) then follows from Lemma 4.7, which concludes the proof. $\square$

## 5. NUMERICAL EXPERIMENTS

In this section, we present two examples which numerically illustrate the analysis and observations of the previous sections. We use the following notations to denote errors of the correction step,

$$e_{u1} := \|u_1^N - u^N\|_{L^2(\Omega_f)},$$

$$e_{w1} := \|w_1^N - w^N\|_{L^2(\Omega_s)},$$

$$e_{du1} := \|\nabla(u_1^N - u^N)\|_{L^2(\Omega_f)},$$

$$e_{\lambda} := \|\lambda_0^N - l^N\|_{L^2(\Sigma)},$$

$$e_{1,\lambda} := \|(\lambda_0^N - l^N) - (\lambda_0^{N-1} - l^{N-1})\|_{L^2(\Sigma)}.$$

All numerical experiments are performed using FEniCS and multiphenics [2, 4].

*Example 5.1 (Non-horizontal interface).* In this example, we test the defect-correction scheme (2.7) in a problem with non-horizontal interface. We let $\Omega = (0, 1)^2$ and we take the solution to be

$$w = u = e^{-2\pi^2 t} \cos(\pi x_1) \sin(\pi x_2).$$

We let $\Sigma$ be defined as the straight line connecting $(0, 0.25)$ and $(1, 0.75)$. We then define $\Omega_s$ as the region above $\Sigma$ and $\Omega_f$ as the region below $\Sigma$ (see Figure 1 for an illustration). The diffusion parameters are both set to be one, $\nu_s = \nu_f = 1$. We choose $h = \Delta t, T = 0.25, \alpha = 4$ and use the piece-wise linear finite element method for the spatial discretization. The results are reported in Table 1. As we can observe, second-order convergence $\mathcal{O}((\Delta t)^2)$ is obtained for $e_{u1}$ and $e_{w1}$, which is in agreement with the theoretical findings provided by Theorem 4.8. Finally, we also observe that $e_{1,\lambda}$ converges as $\mathcal{O}((\Delta t)^2)$. This matches with Corollary 3.2.
<table><caption>TABLE 1. Errors and convergence rates for Example 5.1.</caption><thead><tr><th>Δt</th><th>e<sub>u1</sub></th><th>rates</th><th>e<sub>w1</sub></th><th>rates</th><th>e<sub>λ</sub></th><th>rates</th><th>e<sub>1,λ</sub></th><th>rates</th><th>e<sub>du1</sub></th><th>rates</th></tr></thead><tbody><tr><td>(1/2)<sup>2</sup></td><td>2.47e-02</td><td>-</td><td>1.14e-01</td><td>-</td><td>8.51e-01</td><td>-</td><td>8.54e-01</td><td>-</td><td>2.29e-01</td><td>-</td></tr><tr><td>(1/2)<sup>3</sup></td><td>1.44e-02</td><td>0.78</td><td>5.68e-02</td><td>1.01</td><td>5.01e-01</td><td>0.76</td><td>3.81e-01</td><td>1.16</td><td>7.43e-02</td><td>1.63</td></tr><tr><td>(1/2)<sup>4</sup></td><td>1.36e-02</td><td>0.08</td><td>1.35e-02</td><td>2.08</td><td>1.94e-01</td><td>1.37</td><td>1.45e-01</td><td>1.39</td><td>7.31e-02</td><td>0.02</td></tr><tr><td>(1/2)<sup>5</sup></td><td>5.89e-03</td><td>1.21</td><td>3.20e-03</td><td>2.07</td><td>5.34e-02</td><td>1.86</td><td>2.44e-02</td><td>2.58</td><td>2.82e-02</td><td>1.37</td></tr><tr><td>(1/2)<sup>6</sup></td><td>1.50e-03</td><td>1.98</td><td>1.07e-03</td><td>1.59</td><td>1.84e-02</td><td>1.54</td><td>4.38e-03</td><td>2.48</td><td>6.69e-03</td><td>2.08</td></tr><tr><td>(1/2)<sup>7</sup></td><td>3.59e-04</td><td>2.06</td><td>2.75e-04</td><td>1.95</td><td>7.49e-03</td><td>1.30</td><td>9.10e-04</td><td>2.27</td><td>1.56e-03</td><td>2.10</td></tr><tr><td>(1/2)<sup>8</sup></td><td>8.69e-05</td><td>2.05</td><td>6.84e-05</td><td>2.01</td><td>3.39e-03</td><td>1.14</td><td>2.07e-04</td><td>2.13</td><td>3.71e-04</td><td>2.07</td></tr><tr><td>(1/2)<sup>9</sup></td><td>2.13e-05</td><td>2.03</td><td>1.69e-05</td><td>2.01</td><td>1.61e-03</td><td>1.07</td><td>4.95e-05</td><td>2.07</td><td>9.01e-05</td><td>2.04</td></tr></tbody></table>

Example 5.2 (Different Viscosity). In this example, the defect-correction scheme (2.7) is tested with an exact solution involving different diffusion parameters in $\Omega_f$ and $\Omega_s$. We let $\Omega = (0,1)^2$, $\Sigma$ to be the straight line $x_2 = 0.75$ and we take the solution to be

$$ (5.1) \qquad \mathbf{w} = e^{-2\pi^2 t} \cos(\pi x_1) \sin\left(4\pi \frac{\nu_f}{\nu_s}(x_2 - 0.75)\right), \quad \mathbf{u} = e^{-2\pi^2 t} \cos(\pi x_1) \sin(4\pi(x_2 - 0.75)). $$

Notice that $\mathbf{w} = \mathbf{u}$ and $\nu_s \nabla \mathbf{w} \cdot \mathbf{n}_s + \nu_f \nabla \mathbf{u} \cdot \mathbf{n}_f = 0$ on the interface $\Sigma$. We also have $\mathbf{u} \in V_f$ and $\mathbf{w} \in V_s$ if $\frac{\nu_f}{\nu_s}$ is an integer. We calculate $g_1 = \partial_t \mathbf{u} - \nu_f \Delta \mathbf{u}$, $g_2 = \partial_t \mathbf{w} - \nu_s \Delta \mathbf{w}$ and solve the parabolic-parabolic interface problem (2.1) with nonhomegeneous right-hand sides $g_1$ and $g_2$. We choose $\nu_f = 2$, $\nu_s = 1$, $h = \Delta t$, $T = 0.25$, $\alpha = 4$ and use the piece-wise linear finite element method for spatial discretization. We use the same notations as that of Example 5.1.

<table><caption>TABLE 2. Error and convergence rates for Example 5.2</caption><thead><tr><th>Δt</th><th>e<sub>u1</sub></th><th>rates</th><th>e<sub>w1</sub></th><th>rates</th><th>e<sub>λ</sub></th><th>rates</th><th>e<sub>1,λ</sub></th><th>rates</th><th>e<sub>du</sub></th><th>rates</th></tr></thead><tbody><tr><td>(1/2)<sup>2</sup></td><td>5.21e-01</td><td>-</td><td>3.49e-01</td><td>-</td><td>1.35e+00</td><td>-</td><td>7.90e-01</td><td>-</td><td>2.63e+00</td><td>-</td></tr><tr><td>(1/2)<sup>3</sup></td><td>1.38e-01</td><td>1.92</td><td>3.07e-01</td><td>0.19</td><td>1.23e+00</td><td>0.13</td><td>4.11e-01</td><td>0.94</td><td>1.70e+00</td><td>0.63</td></tr><tr><td>(1/2)<sup>4</sup></td><td>3.81e-02</td><td>1.86</td><td>7.14e-02</td><td>2.10</td><td>8.42e-01</td><td>0.55</td><td>9.23e-02</td><td>2.15</td><td>4.89e-01</td><td>1.79</td></tr><tr><td>(1/2)<sup>5</sup></td><td>9.85e-03</td><td>1.95</td><td>1.92e-02</td><td>1.90</td><td>4.65e-01</td><td>0.86</td><td>1.69e-02</td><td>2.45</td><td>1.27e-01</td><td>1.95</td></tr><tr><td>(1/2)<sup>6</sup></td><td>2.48e-03</td><td>1.99</td><td>4.87e-03</td><td>1.98</td><td>2.40e-01</td><td>0.96</td><td>3.89e-03</td><td>2.12</td><td>3.21e-02</td><td>1.98</td></tr><tr><td>(1/2)<sup>7</sup></td><td>6.22e-04</td><td>2.00</td><td>1.22e-03</td><td>1.99</td><td>1.22e-01</td><td>0.98</td><td>9.54e-04</td><td>2.03</td><td>8.06e-03</td><td>1.99</td></tr><tr><td>(1/2)<sup>8</sup></td><td>1.55e-04</td><td>2.00</td><td>3.06e-04</td><td>2.00</td><td>6.15e-02</td><td>0.99</td><td>2.40e-04</td><td>1.99</td><td>2.02e-03</td><td>2.00</td></tr><tr><td>(1/2)<sup>9</sup></td><td>3.89e-05</td><td>2.00</td><td>7.65e-05</td><td>2.00</td><td>3.09e-02</td><td>0.99</td><td>6.02e-05</td><td>1.99</td><td>5.05e-04</td><td>2.00</td></tr></tbody></table>

We can clearly see $O((\Delta t)^2)$ convergence rates for $e_{u1}$, $e_{w1}$ and $e_{du}$, which agrees with our analysis. One also observes the optimal convergence rates for $e_{1,\lambda}$ which match with Corollary 3.2.

## 6. THE DIRICHLET INTERFACE PROBLEM

In this section, we briefly discuss the defect-correction methods for the following parabolic\parabolic interface problem with Dirichlet boundary conditions on all sides (cf. [14]). As we shall demonstrate there is an inconsistency in the formulation at the the end points of $\Sigma$ when $\Sigma$ is not perpendicular to the two sides of $\partial\Omega$. Then, this inconsistency makes the second-order rates of $\Lambda_0^{n+1} - \Lambda_0^n$ deteriorate, which means that the Corollary 3.2 is not valid in this case. This in turn destroys the second-order convergence rate for the correction method measured in the $H^1$ norm. Nonetheless, when measured in the $L^2$ norm the approximations of the correction method seem to be second-order accurate.
<table><caption>TABLE 1. Errors and convergence rates for Example 5.1.</caption><thead><tr><th>Δt</th><th>e<sub>u1</sub></th><th>rates</th><th>e<sub>w1</sub></th><th>rates</th><th>e<sub>λ</sub></th><th>rates</th><th>e<sub>1,λ</sub></th><th>rates</th><th>e<sub>du1</sub></th><th>rates</th></tr></thead><tbody><tr><td>(1/2)<sup>2</sup></td><td>2.47e-02</td><td>-</td><td>1.14e-01</td><td>-</td><td>8.51e-01</td><td>-</td><td>8.54e-01</td><td>-</td><td>2.29e-01</td><td>-</td></tr><tr><td>(1/2)<sup>3</sup></td><td>1.44e-02</td><td>0.78</td><td>5.68e-02</td><td>1.01</td><td>5.01e-01</td><td>0.76</td><td>3.81e-01</td><td>1.16</td><td>7.43e-02</td><td>1.63</td></tr><tr><td>(1/2)<sup>4</sup></td><td>1.36e-02</td><td>0.08</td><td>1.35e-02</td><td>2.08</td><td>1.94e-01</td><td>1.37</td><td>1.45e-01</td><td>1.39</td><td>7.31e-02</td><td>0.02</td></tr><tr><td>(1/2)<sup>5</sup></td><td>5.89e-03</td><td>1.21</td><td>3.20e-03</td><td>2.07</td><td>5.34e-02</td><td>1.86</td><td>2.44e-02</td><td>2.58</td><td>2.82e-02</td><td>1.37</td></tr><tr><td>(1/2)<sup>6</sup></td><td>1.50e-03</td><td>1.98</td><td>1.07e-03</td><td>1.59</td><td>1.84e-02</td><td>1.54</td><td>4.38e-03</td><td>2.48</td><td>6.69e-03</td><td>2.08</td></tr><tr><td>(1/2)<sup>7</sup></td><td>3.59e-04</td><td>2.06</td><td>2.75e-04</td><td>1.95</td><td>7.49e-03</td><td>1.30</td><td>9.10e-04</td><td>2.27</td><td>1.56e-03</td><td>2.10</td></tr><tr><td>(1/2)<sup>8</sup></td><td>8.69e-05</td><td>2.05</td><td>6.84e-05</td><td>2.01</td><td>3.39e-03</td><td>1.14</td><td>2.07e-04</td><td>2.13</td><td>3.71e-04</td><td>2.07</td></tr><tr><td>(1/2)<sup>9</sup></td><td>2.13e-05</td><td>2.03</td><td>1.69e-05</td><td>2.01</td><td>1.61e-03</td><td>1.07</td><td>4.95e-05</td><td>2.07</td><td>9.01e-05</td><td>2.04</td></tr></tbody></table>

Example 5.2 (Different Viscosity). In this example, the defect-correction scheme (2.7) is tested with an exact solution involving different diffusion parameters in $\Omega_f$ and $\Omega_s$. We let $\Omega = (0,1)^2$, $\Sigma$ to be the straight line $x_2 = 0.75$ and we take the solution to be

$$ (5.1) \qquad w = e^{-2\pi^2 t} \cos(\pi x_1) \sin\left(4\pi \frac{\nu_f}{\nu_s}(x_2 - 0.75)\right), \quad u = e^{-2\pi^2 t} \cos(\pi x_1) \sin\left(4\pi (x_2 - 0.75)\right). $$

Notice that $w = u$ and $\nu_s \nabla w \cdot n_s + \nu_f \nabla u \cdot n_f = 0$ on the interface $\Sigma$. We also have $u \in V_f$ and $w \in V_s$ if $\nu_f/\nu_s$ is an integer. We calculate $g_1 = \partial_t u - \nu_f \Delta u$, $g_2 = \partial_t w - \nu_s \Delta w$ and solve the parabolic-parabolic interface problem (2.1) with nonhomegeneous right-hand sides $g_1$ and $g_2$. We choose $\nu_f = 2$, $\nu_s = 1$, $h = \Delta t$, $T = 0.25$, $\alpha = 4$ and use the piece-wise linear finite element method for spatial discretization. We use the same notations as that of Example 5.1.

<table><caption>TABLE 2. Error and convergence rates for Example 5.2</caption><thead><tr><th>Δt</th><th>e<sub>u1</sub></th><th>rates</th><th>e<sub>w1</sub></th><th>rates</th><th>e<sub>λ</sub></th><th>rates</th><th>e<sub>1,λ</sub></th><th>rates</th><th>e<sub>du</sub></th><th>rates</th></tr></thead><tbody><tr><td>(1/2)<sup>2</sup></td><td>5.21e-01</td><td>-</td><td>3.49e-01</td><td>-</td><td>1.35e+00</td><td>-</td><td>7.90e-01</td><td>-</td><td>2.63e+00</td><td>-</td></tr><tr><td>(1/2)<sup>3</sup></td><td>1.38e-01</td><td>1.92</td><td>3.07e-01</td><td>0.19</td><td>1.23e+00</td><td>0.13</td><td>4.11e-01</td><td>0.94</td><td>1.70e+00</td><td>0.63</td></tr><tr><td>(1/2)<sup>4</sup></td><td>3.81e-02</td><td>1.86</td><td>7.14e-02</td><td>2.10</td><td>8.42e-01</td><td>0.55</td><td>9.23e-02</td><td>2.15</td><td>4.89e-01</td><td>1.79</td></tr><tr><td>(1/2)<sup>5</sup></td><td>9.85e-03</td><td>1.95</td><td>1.92e-02</td><td>1.90</td><td>4.65e-01</td><td>0.86</td><td>1.69e-02</td><td>2.45</td><td>1.27e-01</td><td>1.95</td></tr><tr><td>(1/2)<sup>6</sup></td><td>2.48e-03</td><td>1.99</td><td>4.87e-03</td><td>1.98</td><td>2.40e-01</td><td>0.96</td><td>3.89e-03</td><td>2.12</td><td>3.21e-02</td><td>1.98</td></tr><tr><td>(1/2)<sup>7</sup></td><td>6.22e-04</td><td>2.00</td><td>1.22e-03</td><td>1.99</td><td>1.22e-01</td><td>0.98</td><td>9.54e-04</td><td>2.03</td><td>8.06e-03</td><td>1.99</td></tr><tr><td>(1/2)<sup>8</sup></td><td>1.55e-04</td><td>2.00</td><td>3.06e-04</td><td>2.00</td><td>6.15e-02</td><td>0.99</td><td>2.40e-04</td><td>1.99</td><td>2.02e-03</td><td>2.00</td></tr><tr><td>(1/2)<sup>9</sup></td><td>3.89e-05</td><td>2.00</td><td>7.65e-05</td><td>2.00</td><td>3.09e-02</td><td>0.99</td><td>6.02e-05</td><td>1.99</td><td>5.05e-04</td><td>2.00</td></tr></tbody></table>

We can clearly see $O((\Delta t)^2)$ convergence rates for $e_{u1}$, $e_{w1}$ and $e_{du}$, which agrees with our analysis. One also observes the optimal convergence rates for $e_{1,\lambda}$ which match with Corollary 3.2.

## 6. THE DIRICHLET INTERFACE PROBLEM

In this section, we briefly discuss the defect-correction methods for the following parabolic\parabolic interface problem with Dirichlet boundary conditions on all sides (cf. [14]). As we shall demonstrate there is an inconsistency in the formulation at the the end points of $\Sigma$ when $\Sigma$ is not perpendicular to the two sides of $\partial\Omega$. Then, this inconsistency makes the second-order rates of $\Lambda_0^{n+1} - \Lambda_0^n$ deteriorate, which means that the Corollary 3.2 is not valid in this case. This in turn destroys the second-order convergence rate for the correction method measured in the $H^1$ norm. Nonetheless, when measured in the $L^2$ norm the approximations of the correction method seem to be second-order accurate.
The Dirichlet problem is as follows.

$$ (6.1a) \qquad \partial_t u - \nu_f \Delta u = 0, \quad \text{in } (0, T) \times \Omega_f, $$

$$ (6.1b) \qquad u(0, x) = u_0(x), \quad \text{in } \Omega_f, $$

$$ (6.1c) \qquad u = 0, \quad \text{on } (0, T) \times \partial \Omega_f \setminus \Sigma, $$

$$ (6.2a) \qquad \partial_t w - \nu_s \Delta w = 0, \quad \text{in } (0, T) \times \Omega_s, $$

$$ (6.2b) \qquad w(0, x) = w_0(x), \quad \text{in } \Omega_s, $$

$$ (6.2c) \qquad w = 0, \quad \text{on } (0, T) \times \partial \Omega_s \setminus \Sigma, $$

and

$$ (6.3a) \qquad w - u = 0, \quad \text{in } (0, T) \times \Sigma, $$

$$ (6.3b) \qquad \nu_s \nabla w \cdot n_s + \nu_f \nabla u \cdot n_f = 0 \quad \text{in } (0, T) \times \Sigma, $$

where $n_f$ and $n_s$ are the outward facing normal vectors for $\Omega_f$ and $\Omega_s$, respectively.

We define the following spaces:

$$ (6.4a) \qquad \tilde{V}_f = \{v \in H^1(\Omega_f) : v = 0 \text{ on } \partial\Omega_f \setminus \Sigma\}, $$

$$ (6.4b) \qquad \tilde{V}_s = \{v \in H^1(\Omega_s) : v = 0 \text{ on } \partial\Omega_s \setminus \Sigma\}, $$

$$ (6.4c) \qquad \tilde{V}_g = L^2(\Sigma). $$

Note that the variational form for (6.1)-(6.3) is the same as that of (2.1). The defect-correction method for (6.1)-(6.3) is almost identical to that of (2.1), the only difference is that we replace the $V$ spaces with the $\tilde{V}$ spaces.

There is an inconsistency of $\lambda_0$ at the end points of $\Sigma$. As in (2.6c) we have $\lambda_0^{n+1} - \lambda_0^n = \alpha(w_0^{n+1} - u_0^{n+1})$ on $\Sigma$. Moreover, $w_0^{n+1}$ and $u_0^{n+1}$ have zero Dirichlet boundary conditions at two end points of $\Sigma$ which force $\lambda_0^{n+1} = \lambda_0^n$ at those endpoints. Therefore, $\lambda_0^{n+1}$ does not change in time hence it does not approximate $|n+1|$ at those endpoints and consequently affect the convergence rates on $\Sigma$.

*Example 6.1* (Non-horizontal Interface for the Dirichlet parabolic-parabolic problems). We now present a numerical example which shows this phenomena. We let $\Omega = (0, 1)^2$ and we take the solution to be

$$ w = u = e^{-2\pi^2 t} \sin(\pi x_1) \sin(\pi x_2). $$

We use the same notations and parameters as those in Example 5.1. The convergence results are shown in Table 3.

<table><caption>TABLE 3. Errors and convergence rates for Example 6.1</caption><thead><tr><th>Δt</th><th>e<sub>u1</sub></th><th>rates</th><th>e<sub>w1</sub></th><th>rates</th><th>e<sub>λ</sub></th><th>rates</th><th>e<sub>1,λ</sub></th><th>rates</th><th>e<sub>du1</sub></th><th>rates</th></tr></thead><tbody><tr><td>(1/2)<sup>2</sup></td><td>1.99e-02</td><td>-</td><td>4.18e-02</td><td>-</td><td>7.04e-01</td><td>-</td><td>6.99e-01</td><td>-</td><td>2.21e-01</td><td>-</td></tr><tr><td>(1/2)<sup>3</sup></td><td>7.58e-03</td><td>1.39</td><td>5.82e-02</td><td>-0.48</td><td>5.04e-01</td><td>0.48</td><td>1.18e-01</td><td>2.57</td><td>4.08e-02</td><td>2.43</td></tr><tr><td>(1/2)<sup>4</sup></td><td>1.40e-02</td><td>-0.88</td><td>1.18e-02</td><td>2.30</td><td>2.75e-01</td><td>0.88</td><td>1.13e-01</td><td>0.06</td><td>7.61e-02</td><td>-0.90</td></tr><tr><td>(1/2)<sup>5</sup></td><td>7.05e-03</td><td>0.99</td><td>3.55e-03</td><td>1.73</td><td>1.54e-01</td><td>0.83</td><td>2.61e-02</td><td>2.11</td><td>3.48e-02</td><td>1.13</td></tr><tr><td>(1/2)<sup>6</sup></td><td>1.83e-03</td><td>1.94</td><td>1.41e-03</td><td>1.33</td><td>1.02e-01</td><td>0.59</td><td>6.20e-03</td><td>2.07</td><td>9.32e-03</td><td>1.90</td></tr><tr><td>(1/2)<sup>7</sup></td><td>4.47e-04</td><td>2.04</td><td>3.80e-04</td><td>1.89</td><td>7.03e-02</td><td>0.54</td><td>1.84e-03</td><td>1.75</td><td>2.74e-03</td><td>1.76</td></tr><tr><td>(1/2)<sup>8</sup></td><td>1.08e-04</td><td>2.05</td><td>9.63e-05</td><td>1.98</td><td>4.94e-02</td><td>0.51</td><td>5.95e-04</td><td>1.63</td><td>1.01e-03</td><td>1.45</td></tr><tr><td>(1/2)<sup>9</sup></td><td>2.65e-05</td><td>2.03</td><td>2.41e-05</td><td>2.00</td><td>3.47e-02</td><td>0.51</td><td>2.00e-04</td><td>1.57</td><td>4.43e-04</td><td>1.19</td></tr></tbody></table>

We can clearly see that $e_{1,\lambda}$ is not second-order accurate and hence Corollary 3.2 is not valid in this case. Moreover, the gradient error of $u$ of the correction method, $e_{du1}$, is not second-order. Finally, we
The Dirichlet problem is as follows.

$$ (6.1a) \qquad \partial_t u - \nu_f \Delta u = 0, \quad \text{in } (0, T) \times \Omega_f, $$

$$ (6.1b) \qquad u(0, x) = u_0(x), \quad \text{in } \Omega_f, $$

$$ (6.1c) \qquad u = 0, \quad \text{on } (0, T) \times \partial \Omega_f \setminus \Sigma, $$

$$ (6.2a) \qquad \partial_t w - \nu_s \Delta w = 0, \quad \text{in } (0, T) \times \Omega_s, $$

$$ (6.2b) \qquad w(0, x) = w_0(x), \quad \text{in } \Omega_s, $$

$$ (6.2c) \qquad w = 0, \quad \text{on } (0, T) \times \partial \Omega_s \setminus \Sigma, $$

and

$$ (6.3a) \qquad w - u = 0, \quad \text{in } (0, T) \times \Sigma, $$

$$ (6.3b) \qquad \nu_s \nabla w \cdot n_s + \nu_f \nabla u \cdot n_f = 0 \quad \text{in } (0, T) \times \Sigma, $$

where $n_f$ and $n_s$ are the outward facing normal vectors for $\Omega_f$ and $\Omega_s$, respectively.

We define the following spaces:

$$ (6.4a) \qquad \tilde{V}_f = \{v \in H^1(\Omega_f) : v = 0 \text{ on } \partial\Omega_f \setminus \Sigma\}, $$

$$ (6.4b) \qquad \tilde{V}_s = \{v \in H^1(\Omega_s) : v = 0 \text{ on } \partial\Omega_s \setminus \Sigma\}, $$

$$ (6.4c) \qquad \tilde{V}_g = L^2(\Sigma). $$

Note that the variational form for (6.1)-(6.3) is the same as that of (2.1). The defect-correction method for (6.1)-(6.3) is almost identical to that of (2.1), the only difference is that we replace the $V$ spaces with the $\tilde{V}$ spaces.

There is an inconsistency of $\lambda_0$ at the end points of $\Sigma$. As in (2.6c) we have $\lambda_0^{n+1} - \lambda_0^n = \alpha(w_0^{n+1} - u_0^{n+1})$ on $\Sigma$. Moreover, $w_0^{n+1}$ and $u_0^{n+1}$ have zero Dirichlet boundary conditions at two end points of $\Sigma$ which force $\lambda_0^{n+1} = \lambda_0^n$ at those endpoints. Therefore, $\lambda_0^{n+1}$ does not change in time hence it does not approximate $|n+1|$ at those endpoints and consequently affect the convergence rates on $\Sigma$.

*Example 6.1* (Non-horizontal Interface for the Dirichlet parabolic-parabolic problems). We now present a numerical example which shows this phenomena. We let $\Omega = (0, 1)^2$ and we take the solution to be

$$ w = u = e^{-2\pi^2 t} \sin(\pi x_1) \sin(\pi x_2). $$

We use the same notations and parameters as those in Example 5.1. The convergence results are shown in Table 3.

<table><caption>TABLE 3. Errors and convergence rates for Example 6.1</caption><thead><tr><th>Δt</th><th>e<sub>u1</sub></th><th>rates</th><th>e<sub>w1</sub></th><th>rates</th><th>e<sub>λ</sub></th><th>rates</th><th>e<sub>1,λ</sub></th><th>rates</th><th>e<sub>du1</sub></th><th>rates</th></tr></thead><tbody><tr><td>(1/2)<sup>2</sup></td><td>1.99e-02</td><td>-</td><td>4.18e-02</td><td>-</td><td>7.04e-01</td><td>-</td><td>6.99e-01</td><td>-</td><td>2.21e-01</td><td>-</td></tr><tr><td>(1/2)<sup>3</sup></td><td>7.58e-03</td><td>1.39</td><td>5.82e-02</td><td>-0.48</td><td>5.04e-01</td><td>0.48</td><td>1.18e-01</td><td>2.57</td><td>4.08e-02</td><td>2.43</td></tr><tr><td>(1/2)<sup>4</sup></td><td>1.40e-02</td><td>-0.88</td><td>1.18e-02</td><td>2.30</td><td>2.75e-01</td><td>0.88</td><td>1.13e-01</td><td>0.06</td><td>7.61e-02</td><td>-0.90</td></tr><tr><td>(1/2)<sup>5</sup></td><td>7.05e-03</td><td>0.99</td><td>3.55e-03</td><td>1.73</td><td>1.54e-01</td><td>0.83</td><td>2.61e-02</td><td>2.11</td><td>3.48e-02</td><td>1.13</td></tr><tr><td>(1/2)<sup>6</sup></td><td>1.83e-03</td><td>1.94</td><td>1.41e-03</td><td>1.33</td><td>1.02e-01</td><td>0.59</td><td>6.20e-03</td><td>2.07</td><td>9.32e-03</td><td>1.90</td></tr><tr><td>(1/2)<sup>7</sup></td><td>4.47e-04</td><td>2.04</td><td>3.80e-04</td><td>1.89</td><td>7.03e-02</td><td>0.54</td><td>1.84e-03</td><td>1.75</td><td>2.74e-03</td><td>1.76</td></tr><tr><td>(1/2)<sup>8</sup></td><td>1.08e-04</td><td>2.05</td><td>9.63e-05</td><td>1.98</td><td>4.94e-02</td><td>0.51</td><td>5.95e-04</td><td>1.63</td><td>1.01e-03</td><td>1.45</td></tr><tr><td>(1/2)<sup>9</sup></td><td>2.65e-05</td><td>2.03</td><td>2.41e-05</td><td>2.00</td><td>3.47e-02</td><td>0.51</td><td>2.00e-04</td><td>1.57</td><td>4.43e-04</td><td>1.19</td></tr></tbody></table>

We can clearly see that $e_{1,\lambda}$ is not second-order accurate and hence Corollary 3.2 is not valid in this case. Moreover, the gradient error of $u$ of the correction method, $e_{du1}$, is not second-order. Finally, we
<figure><img src="image_5.png" /><figcaption>FIGURE 3. The domains $\Omega_f$ and $\Omega_s$ with non-horizontal interface $\Sigma$.</figcaption></figure>

observe that $e_{u1}$ and $e_{w1}$ are second-order accurate. It is worthwhile mentioning that if the interface is perpendicular to two of the sides of $\partial\Omega$ then there is no inconsistency for the Dirichlet problem and, numerically, one regains second-order convergence for all the variables above. We omit the numerical experiments here.

We finish this section by stating that there is a remedy for the Dirichlet problem, when $\nu_s = \nu_f = \nu$, in the case $\Sigma$ is straight line and that is not parallel to the two sides of $\Omega$. To be specific, we denote $\tau$ as the tangent vector of $\Sigma$ and $s$ as the tangent vector of the side of $\Omega$. It is easy to see that there exists constants $a$ and $b$ such that $\mathbf{n}_f = a\tau + bs$, and hence $\nabla u \cdot \mathbf{n}_f = a\nabla u \cdot \tau + b\nabla u \cdot s$. Note that $\mathbf{n}_s = -a\tau - bs$.

**Remark 6.2.** We take $\Omega = (0, 1)^2$ and let the interface $\Sigma$ be the straight line connecting the point $(0, 0.25)$ with the point $(1, 0.75)$. In this case, the unit tangent to $\Sigma$ is $\tau = [2/\sqrt{5}, 1/\sqrt{5}]^t$. The normal $\mathbf{n}_f$ is $\mathbf{n}_f = [-1/\sqrt{5}, 2/\sqrt{5}]^t$. We let $\mathbf{s} = [0, 1]^t$, hence we can write $\mathbf{n}_f = \frac{1}{2}\tau + \frac{\sqrt{5}}{2}\mathbf{s}$. See Figure 3 for an example.

We then define the following modified defect-correction method.

**Prediction step.** Find $w_0^{n+1} \in \tilde{V}_s$, $u_0^{n+1} \in \tilde{V}_f$, and $\lambda_0^{n+1} \in \tilde{V}_g$ such that for $n \ge 0$,

$$ (6.5a) \qquad (\partial_{\Delta t} w_0^{n+1}, z)_s + \nu(\nabla w_0^{n+1}, \nabla z)_s + b\alpha \langle w_0^{n+1} - u_0^n, z \rangle + a\nu \langle \nabla w_0^{n+1} \cdot \tau, z \rangle + \langle \lambda^n, z \rangle = 0, \qquad \forall z \in \tilde{V}_s, $$

$$ (6.5b) \qquad (\partial_{\Delta t} u_0^{n+1}, v)_f + \nu(\nabla u_0^{n+1}, \nabla v)_f - \langle \lambda_0^{n+1}, v \rangle - a\nu \langle \nabla u_0^{n+1} \cdot \tau, v \rangle = 0, \quad \forall v \in \tilde{V}_f, $$

$$ (6.5c) \qquad \langle (\lambda_0^{n+1} - \lambda_0^n) + \alpha(u_0^{n+1} - w_0^{n+1}), \mu \rangle = 0, \qquad \forall \mu \in \tilde{V}_g. $$

**Correction step.** Find $w_1^{n+1} \in \tilde{V}_s$, $u_1^{n+1} \in \tilde{V}_f$, and $\lambda_1^{n+1} \in \tilde{V}_g$ such that for $n \ge 0$,

$$ (\partial_{\Delta t} w_1^{n+1}, z)_s + \nu(\nabla w_1^{n+1}, \nabla z)_s + b\alpha \langle w_1^{n+1} - u_1^n, z \rangle + a\nu \langle \nabla w_1^{n+1} \cdot \tau, z \rangle + \langle \lambda_1^n, z \rangle $$

$$ (6.6a) \qquad = \nu(\nabla w_0^{n+1}, \nabla z)_s + b\alpha \langle w_0^{n+1} - w_0^n, z \rangle + \langle \lambda_0^n, z \rangle + a\nu \langle \nabla w_0^{n+1} \cdot \tau, z \rangle - \nu(\nabla w_0^{n+\frac{1}{2}}, \nabla z)_s - \langle \lambda_0^{n+\frac{1}{2}}, z \rangle - a\nu \langle \nabla w_0^{n+\frac{1}{2}} \cdot \tau, z \rangle, \qquad \forall z \in \tilde{V}_s, $$

$$ (6.6b) \qquad = \nu(\nabla u_0^{n+1}, \nabla v)_f - \langle \lambda_0^{n+1}, v \rangle - a\nu \langle \nabla u_0^{n+1} \cdot \tau, v \rangle - \nu(\nabla u_0^{n+\frac{1}{2}}, \nabla v)_f + \langle \lambda_0^{n+\frac{1}{2}}, v \rangle + a\nu \langle \nabla u_0^{n+\frac{1}{2}} \cdot \tau, v \rangle, \qquad \forall v \in \tilde{V}_f, $$

$$ (6.6c) \qquad \langle (\lambda_0^{n+1} - \lambda_0^n) + \alpha(u_1^{n+1} - w_1^{n+1}), \mu \rangle = \langle \lambda_0^{n+1} - \lambda_0^n, \mu \rangle \qquad \forall \mu \in \tilde{V}_g. $$

**Remark 6.3.** The difference between the original correction method and the modified correction method (6.5)-(6.6) is the definition of the Lagrange multiplier $\lambda^{n+1}$. The Lagrange multiplier $\lambda^{n+1}$ in (2.7) approximates $\nu\nabla u^{n+1} \cdot \mathbf{n}_f$ while the Lagrange multiplier $\lambda^{n+1}$ in (6.6) approximates $b\nu\nabla u^{n+1} \cdot \mathbf{s}$.
<figure><img src="image_5.png" /><figcaption>FIGURE 3. The domains $\Omega_f$ and $\Omega_s$ with non-horizontal interface $\Sigma$.</figcaption></figure>

observe that $e_{u1}$ and $e_{w1}$ are second-order accurate. It is worthwhile mentioning that if the interface is perpendicular to two of the sides of $\partial\Omega$ then there is no inconsistency for the Dirichlet problem and, numerically, one regains second-order convergence for all the variables above. We omit the numerical experiments here.

We finish this section by stating that there is a remedy for the Dirichlet problem, when $\nu_s = \nu_f = \nu$, in the case $\Sigma$ is straight line and that is not parallel to the two sides of $\Omega$. To be specific, we denote $\tau$ as the tangent vector of $\Sigma$ and $s$ as the tangent vector of the side of $\Omega$. It is easy to see that there exists constants $a$ and $b$ such that $\mathbf{n}_f = a\tau + bs$, and hence $\nabla u \cdot \mathbf{n}_f = a\nabla u \cdot \tau + b\nabla u \cdot s$. Note that $\mathbf{n}_s = -a\tau - bs$.

**Remark 6.2.** We take $\Omega = (0, 1)^2$ and let the interface $\Sigma$ be the straight line connecting the point $(0, 0.25)$ with the point $(1, 0.75)$. In this case, the unit tangent to $\Sigma$ is $\tau = [2/\sqrt{5}, 1/\sqrt{5}]^t$. The normal $\mathbf{n}_f$ is $\mathbf{n}_f = [-1/\sqrt{5}, 2/\sqrt{5}]^t$. We let $\mathbf{s} = [0, 1]^t$, hence we can write $\mathbf{n}_f = \frac{1}{2}\tau + \frac{\sqrt{5}}{2}\mathbf{s}$. See Figure 3 for an example.

We then define the following modified defect-correction method.

**Prediction step.** Find $w_0^{n+1} \in \tilde{V}_s$, $u_0^{n+1} \in \tilde{V}_f$, and $\lambda_0^{n+1} \in \tilde{V}_g$ such that for $n \ge 0$,

$$ (6.5a) \qquad (\partial_{\Delta t} w_0^{n+1}, z)_s + \nu(\nabla w_0^{n+1}, \nabla z)_s + b\alpha \langle w_0^{n+1} - u_0^n, z \rangle \\ \qquad \qquad + a\nu \langle \nabla w_0^{n+1} \cdot \tau, z \rangle + \langle \lambda^n, z \rangle = 0, \qquad \forall z \in \tilde{V}_s, $$

$$ (6.5b) \qquad (\partial_{\Delta t} u_0^{n+1}, v)_f + \nu(\nabla u_0^{n+1}, \nabla v)_f - \langle \lambda_0^{n+1}, v \rangle - a\nu \langle \nabla u_0^{n+1} \cdot \tau, v \rangle = 0, \quad \forall v \in \tilde{V}_f, $$

$$ (6.5c) \qquad \langle (\lambda_0^{n+1} - \lambda_0^n) + \alpha(u_0^{n+1} - w_0^{n+1}), \mu \rangle = 0, \qquad \forall \mu \in \tilde{V}_g. $$

**Correction step.** Find $w_1^{n+1} \in \tilde{V}_s$, $u_1^{n+1} \in \tilde{V}_f$, and $\lambda_1^{n+1} \in \tilde{V}_g$ such that for $n \ge 0$,

$$ (\partial_{\Delta t} w_1^{n+1}, z)_s + \nu(\nabla w_1^{n+1}, \nabla z)_s + b\alpha \langle w_1^{n+1} - u_1^n, z \rangle \\ + a\nu \langle \nabla w_1^{n+1} \cdot \tau, z \rangle + \langle \lambda_1^n, z \rangle $$

$$ (6.6a) \qquad = \nu(\nabla w_0^{n+1}, \nabla z)_s + b\alpha \langle w_0^{n+1} - w_0^n, z \rangle + \langle \lambda_0^n, z \rangle + a\nu \langle \nabla w_0^{n+1} \cdot \tau, z \rangle \\ \qquad \qquad - \nu(\nabla w_0^{n+\frac{1}{2}}, \nabla z)_s - \langle \lambda_0^{n+\frac{1}{2}}, z \rangle - a\nu \langle \nabla w_0^{n+\frac{1}{2}} \cdot \tau, z \rangle, \qquad \forall z \in \tilde{V}_s, \\ (\partial_{\Delta t} u_1^{n+1}, v)_f + \nu(\nabla u_1^{n+1}, \nabla v)_f - \langle \lambda_1^{n+1}, v \rangle - a\nu \langle \nabla u_1^{n+1} \cdot \tau, v \rangle $$

$$ (6.6b) \qquad = \nu(\nabla u_0^{n+1}, \nabla v)_f - \langle \lambda_0^{n+1}, v \rangle - a\nu \langle \nabla u_0^{n+1} \cdot \tau, v \rangle - \nu(\nabla u_0^{n+\frac{1}{2}}, \nabla v)_f \\ \qquad \qquad + \langle \lambda_0^{n+\frac{1}{2}}, v \rangle + a\nu \langle \nabla u_0^{n+\frac{1}{2}} \cdot \tau, v \rangle, \qquad \forall v \in \tilde{V}_f, $$

$$ (6.6c) \qquad \langle (\lambda_1^{n+1} - \lambda_1^n) + \alpha(u_1^{n+1} - w_1^{n+1}), \mu \rangle = \langle \lambda_0^{n+1} - \lambda_0^n, \mu \rangle \qquad \forall \mu \in \tilde{V}_g. $$

**Remark 6.3.** The difference between the original correction method and the modified correction method (6.5)-(6.6) is the definition of the Lagrange multiplier $\lambda^{n+1}$. The Lagrange multiplier $\lambda^{n+1}$ in (2.7) approximates $\nu\nabla u^{n+1} \cdot \mathbf{n}_f$ while the Lagrange multiplier $\lambda^{n+1}$ in (6.6) approximates $b\nu\nabla u^{n+1} \cdot \mathbf{s}$.
Note that if the interface $\Sigma$ is not perpendicular to two sides of the domain, $\lambda_0^{n+1}$ does not converge to $\nu\nabla u^{n+1} \cdot n_f$ at those endpoints. That is why we modify the original method such that the Lagrange multiplier approximates $b\nu\nabla u^{n+1} \cdot s$ (which is always 0).

We then show the numerical results of the modified method (6.6) for Example 6.1 in Table 4. We can clearly see the improvement of $e_\lambda$, $e_{1,\lambda}$ and $e_{du1}$ that match with our observation.

<table><caption>TABLE 4. Error and convergence rates for modified defect-correction method (6.6)</caption><thead><tr><th>Δt</th><th>e<sub>u1</sub></th><th>rates</th><th>e<sub>w1</sub></th><th>rates</th><th>e<sub>λ</sub></th><th>rates</th><th>e<sub>1,λ</sub></th><th>rates</th><th>e<sub>du1</sub></th><th>rates</th></tr></thead><tbody><tr><td>(1/2)<sup>2</sup></td><td>2.20e-02</td><td>-</td><td>4.55e-02</td><td>-</td><td>7.65e-01</td><td>-</td><td>7.90e-01</td><td>-</td><td>2.31e-01</td><td>-</td></tr><tr><td>(1/2)<sup>3</sup></td><td>5.29e-03</td><td>2.06</td><td>5.43e-02</td><td>-0.26</td><td>5.19e-01</td><td>0.56</td><td>1.60e-01</td><td>2.30</td><td>3.32e-02</td><td>2.80</td></tr><tr><td>(1/2)<sup>4</sup></td><td>1.29e-02</td><td>-1.29</td><td>1.16e-02</td><td>2.22</td><td>2.33e-01</td><td>1.16</td><td>1.31e-01</td><td>0.29</td><td>7.09e-02</td><td>-1.09</td></tr><tr><td>(1/2)<sup>5</sup></td><td>7.37e-03</td><td>0.81</td><td>4.18e-03</td><td>1.48</td><td>6.57e-02</td><td>1.83</td><td>2.92e-02</td><td>2.17</td><td>3.60e-02</td><td>0.98</td></tr><tr><td>(1/2)<sup>6</sup></td><td>2.05e-03</td><td>1.84</td><td>1.63e-03</td><td>1.36</td><td>1.89e-02</td><td>1.80</td><td>4.59e-03</td><td>2.67</td><td>9.49e-03</td><td>1.92</td></tr><tr><td>(1/2)<sup>7</sup></td><td>5.09e-04</td><td>2.01</td><td>4.42e-04</td><td>1.88</td><td>6.78e-03</td><td>1.48</td><td>8.66e-04</td><td>2.41</td><td>2.32e-03</td><td>2.03</td></tr><tr><td>(1/2)<sup>8</sup></td><td>1.24e-04</td><td>2.04</td><td>1.12e-04</td><td>1.98</td><td>2.85e-03</td><td>1.25</td><td>1.85e-04</td><td>2.23</td><td>5.62e-04</td><td>2.04</td></tr><tr><td>(1/2)<sup>9</sup></td><td>3.03e-05</td><td>2.03</td><td>2.80e-05</td><td>2.00</td><td>1.31e-03</td><td>1.13</td><td>4.25e-05</td><td>2.12</td><td>1.40e-04</td><td>2.01</td></tr></tbody></table>

## 7. FUTURE DIRECTIONS

There are several future directions. One, of course, is to prove Assumption 3.1 for problems with more general interface. This seems to be a challenging problem and new ideas are needed. The second is to extend the methodology to the wave-parabolic interface problem. Preliminary results are very promising. Once we accomplish this, we will be ready to tackle the linear FSI problem. Finally, we will consider the dynamic FSI problem.

## ACKNOWLEDGEMENT

This material is based upon work supported by the National Science Foundation under Grant No. DMS-1929284 while the last author was in residence at the Institute for Computational and Experimental Research in Mathematics in Providence, RI, during the "Numerical PDEs: Analysis, Algorithms, and Data Challenges" program. EB was supported by EPSRC grants EP/J002313/2 and EP/W007460/1. EB and MF were supported by the IMFIBIO Inria-UCL associated team.

## APPENDIX A. NOTATIONS

For readability, we collect various definitions in Section 3 here. We define $\mathcal{Y}$ as,

$$ \begin{aligned} \mathcal{Y} := & \frac{1}{\nu_s} \|\partial_t^3 \mathbf{w}\|_{L^2(0,T;L^2(\Omega_s))}^2 + \left(\frac{1}{\nu_f} + 1\right) \|\partial_t^3 \mathbf{u}\|_{L^2(0,T;L^2(\Omega_f))}^2 \\ & + \left(\frac{\nu_f}{\alpha^2} + \frac{\nu_f^2}{\alpha}\right) \|\partial_t^2 \mathbf{u}\|_{L^2(0,T;H^1(\Omega_f))}^2 + \frac{(\nu_f)^3}{\alpha^2} \|\partial_t^2 \mathbf{u}\|_{L^2(0,T;H^2(\Omega_f))}^2 \\ & + \frac{\alpha^2}{\nu_s} \|\partial_t^2 \mathbf{u}\|_{L^2(0,T;L^2(\Sigma))}^2 + \frac{1}{\nu_f} \|\partial_t^2 \mathbf{l}\|_{L^2(0,T;L^2(\Sigma))}^2 + \frac{\nu_f^2}{\alpha^2} \|\partial_t^2 \mathbf{u}\|_{L^\infty(0,T;H^1(\Omega_f))}^2, \end{aligned} $$
Note that if the interface $\Sigma$ is not perpendicular to two sides of the domain, $\lambda_0^{n+1}$ does not converge to $\nu\nabla u^{n+1} \cdot n_f$ at those endpoints. That is why we modify the original method such that the Lagrange multiplier approximates $b\nu\nabla u^{n+1} \cdot s$ (which is always 0).

We then show the numerical results of the modified method (6.6) for Example 6.1 in Table 4. We can clearly see the improvement of $e_\lambda$, $e_{1,\lambda}$ and $e_{du1}$ that match with our observation.

<table><caption>TABLE 4. Error and convergence rates for modified defect-correction method (6.6)</caption><thead><tr><th>Δt</th><th>e<sub>u1</sub></th><th>rates</th><th>e<sub>w1</sub></th><th>rates</th><th>e<sub>λ</sub></th><th>rates</th><th>e<sub>1,λ</sub></th><th>rates</th><th>e<sub>du1</sub></th><th>rates</th></tr></thead><tbody><tr><td>(1/2)<sup>2</sup></td><td>2.20e-02</td><td>-</td><td>4.55e-02</td><td>-</td><td>7.65e-01</td><td>-</td><td>7.90e-01</td><td>-</td><td>2.31e-01</td><td>-</td></tr><tr><td>(1/2)<sup>3</sup></td><td>5.29e-03</td><td>2.06</td><td>5.43e-02</td><td>-0.26</td><td>5.19e-01</td><td>0.56</td><td>1.60e-01</td><td>2.30</td><td>3.32e-02</td><td>2.80</td></tr><tr><td>(1/2)<sup>4</sup></td><td>1.29e-02</td><td>-1.29</td><td>1.16e-02</td><td>2.22</td><td>2.33e-01</td><td>1.16</td><td>1.31e-01</td><td>0.29</td><td>7.09e-02</td><td>-1.09</td></tr><tr><td>(1/2)<sup>5</sup></td><td>7.37e-03</td><td>0.81</td><td>4.18e-03</td><td>1.48</td><td>6.57e-02</td><td>1.83</td><td>2.92e-02</td><td>2.17</td><td>3.60e-02</td><td>0.98</td></tr><tr><td>(1/2)<sup>6</sup></td><td>2.05e-03</td><td>1.84</td><td>1.63e-03</td><td>1.36</td><td>1.89e-02</td><td>1.80</td><td>4.59e-03</td><td>2.67</td><td>9.49e-03</td><td>1.92</td></tr><tr><td>(1/2)<sup>7</sup></td><td>5.09e-04</td><td>2.01</td><td>4.42e-04</td><td>1.88</td><td>6.78e-03</td><td>1.48</td><td>8.66e-04</td><td>2.41</td><td>2.32e-03</td><td>2.03</td></tr><tr><td>(1/2)<sup>8</sup></td><td>1.24e-04</td><td>2.04</td><td>1.12e-04</td><td>1.98</td><td>2.85e-03</td><td>1.25</td><td>1.85e-04</td><td>2.23</td><td>5.62e-04</td><td>2.04</td></tr><tr><td>(1/2)<sup>9</sup></td><td>3.03e-05</td><td>2.03</td><td>2.80e-05</td><td>2.00</td><td>1.31e-03</td><td>1.13</td><td>4.25e-05</td><td>2.12</td><td>1.40e-04</td><td>2.01</td></tr></tbody></table>

## 7. FUTURE DIRECTIONS

There are several future directions. One, of course, is to prove Assumption 3.1 for problems with more general interface. This seems to be a challenging problem and new ideas are needed. The second is to extend the methodology to the wave-parabolic interface problem. Preliminary results are very promising. Once we accomplish this, we will be ready to tackle the linear FSI problem. Finally, we will consider the dynamic FSI problem.

## ACKNOWLEDGEMENT

This material is based upon work supported by the National Science Foundation under Grant No. DMS-1929284 while the last author was in residence at the Institute for Computational and Experimental Research in Mathematics in Providence, RI, during the "Numerical PDEs: Analysis, Algorithms, and Data Challenges" program. EB was supported by EPSRC grants EP/J002313/2 and EP/W007460/1. EB and MF were supported by the IMFIBIO Inria-UCL associated team.

## APPENDIX A. NOTATIONS

For readability, we collect various definitions in Section 3 here. We define $\mathcal{Y}$ as,

$$ \begin{aligned} \mathcal{Y} := & \frac{1}{\nu_s} \|\partial_t^3 \mathbf{w}\|_{L^2(0,T;L^2(\Omega_s))}^2 + \left(\frac{1}{\nu_f} + 1\right) \|\partial_t^3 \mathbf{u}\|_{L^2(0,T;L^2(\Omega_f))}^2 \\ & + \left(\frac{\nu_f}{\alpha^2} + \frac{\nu_f^2}{\alpha}\right) \|\partial_t^2 \mathbf{u}\|_{L^2(0,T;H^1(\Omega_f))}^2 + \frac{(\nu_f)^3}{\alpha^2} \|\partial_t^2 \mathbf{u}\|_{L^2(0,T;H^2(\Omega_f))}^2 \\ & + \frac{\alpha^2}{\nu_s} \|\partial_t^2 \mathbf{u}\|_{L^2(0,T;L^2(\Sigma))}^2 + \frac{1}{\nu_f} \|\partial_t^2 \mathbf{l}\|_{L^2(0,T;L^2(\Sigma))}^2 + \frac{\nu_f^2}{\alpha^2} \|\partial_t^2 \mathbf{u}\|_{L^\infty(0,T;H^1(\Omega_f))}^2, \end{aligned} $$
define $\mathfrak{Z}$ as,

$$ \mathfrak{Z} := \frac{1}{\nu_s} \|\partial_t^4 \mathbf{w}\|_{L^2(0,T;L^2(\Omega_s))}^2 + \left(\frac{1}{\nu_f} + 1\right) \|\partial_t^4 \mathbf{u}\|_{L^2(0,T;L^2(\Omega_f))}^2 \\ + \left(\frac{\nu_f}{\alpha^2} + \frac{\nu_f^2}{\alpha}\right) \|\partial_t^3 \mathbf{u}\|_{L^2(0,T;H^1(\Omega_f))}^2 + \frac{(\nu_f)^3}{\alpha^2} \|\partial_t^3 \mathbf{u}\|_{L^2(0,T;H^2(\Omega_f))}^2 \\ + \frac{\alpha^2}{\nu_s} \|\partial_t^3 \mathbf{u}\|_{L^2(0,T;L^2(\Sigma))}^2 + \frac{1}{\nu_f} \|\partial_t^3 \mathbf{l}\|_{L^2(0,T;L^2(\Sigma))}^2 + \frac{\nu_f^2}{\alpha^2} \|\partial_t^3 \mathbf{u}\|_{L^\infty(0,T;H^1(\Omega_f))}^2, $$

define Y as,

$$ Y := \frac{1}{\nu_s} \|\partial_x \partial_t^2 \mathbf{w}\|_{L^2(0,T;L^2(\Omega_s))}^2 + \left(\frac{1}{\nu_f} + 1\right) \|\partial_x \partial_t^2 \mathbf{u}\|_{L^2(0,T;L^2(\Omega_f))}^2 \\ + \left(\frac{\nu_f}{\alpha^2} + \frac{\nu_f^2}{\alpha}\right) \|\partial_x \partial_t \mathbf{u}\|_{L^2(0,T;H^1(\Omega_f))}^2 + \frac{(\nu_f)^3}{\alpha^2} \|\partial_x \partial_t \mathbf{u}\|_{L^2(0,T;H^2(\Omega_f))}^2 \\ + \frac{\alpha^2}{\nu_s} \|\partial_x \partial_t \mathbf{u}\|_{L^2(0,T;L^2(\Sigma))}^2 + \frac{1}{\nu_f} \|\partial_x \partial_t \mathbf{l}\|_{L^2(0,T;L^2(\Sigma))}^2 + \frac{\nu_f^2}{\alpha^2} \|\partial_x \partial_t \mathbf{u}\|_{L^\infty(0,T;H^1(\Omega_f))}^2, $$

and define $\mathfrak{Y}$ as,

$$ \mathfrak{Y} := \frac{1}{\nu_s} \|\partial_x \partial_t^3 \mathbf{w}\|_{L^2(0,T;L^2(\Omega_s))}^2 + \left(\frac{1}{\nu_f} + 1\right) \|\partial_x \partial_t^3 \mathbf{u}\|_{L^2(0,T;L^2(\Omega_f))}^2 \\ + \left(\frac{\nu_f}{\alpha^2} + \frac{\nu_f^2}{\alpha}\right) \|\partial_x \partial_t^2 \mathbf{u}\|_{L^2(0,T;H^1(\Omega_f))}^2 + \frac{(\nu_f)^3}{\alpha^2} \|\partial_x \partial_t^2 \mathbf{u}\|_{L^2(0,T;H^2(\Omega_f))}^2 \\ + \frac{\alpha^2}{\nu_s} \|\partial_x \partial_t^2 \mathbf{u}\|_{L^2(0,T;L^2(\Sigma))}^2 + \frac{1}{\nu_f} \|\partial_x \partial_t^2 \mathbf{u}\|_{L^2(0,T;L^2(\Sigma))}^2 + \frac{\nu_f^2}{\alpha^2} \|\partial_x \partial_t^2 \mathbf{u}\|_{L^\infty(0,T;H^1(\Omega_f))}^2. $$

## REFERENCES

[1] M. Aggul, J. M. Connors, D. Erkmen, and A. E. Labovsky. A defect-deferred correction method for fluid-fluid interaction. *SIAM Journal on Numerical Analysis*, 56(4):2484–2512, 2018.

[2] M. Alnaes, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M. E. Rognes, and G. N. Wells. The FEiNCS project version 1.5. *Archive of numerical software*, 3(100), 2015.

[3] S. Badia, F. Nobile, and C. Vergara. Fluid-structure partitioned procedures based on Robin transmission conditions. *Journal of Computational Physics*, 227(14):7027–7051, 2008.

[4] F. Ballarin and G. Rozza. https://multiphenics.github.io/, 2016.

[5] J. W. Banks, W. D. Henshaw, and D. W. Schwendeman. An analysis of a new stable partitioned algorithm for FSI problems. Part I: Incompressible flow and elastic solids. *J. Comput. Phys.*, 269:108–137, 2014.

[6] M. Beneš. Convergence and stability analysis of heterogeneous time step coupling schemes for parabolic problems. *Appl. Numer. Math.*, 121:198–222, 2017.

[7] M. Beneš, A. Nekvinda, and M. K. Yadav. Multi-time-step domain decomposition method with non-matching grids for parabolic problems. *Appl. Math. Comput.*, 267:571–582, 2015.

[8] K. Böhmer. Discrete Newton methods and iterated defect corrections. *Numer. Math.*, 37(2):167–192, 1981.

[9] K. Böhmer, P. W. Hemker, and H. J. Stetter. The defect correction approach. *Defect Correction Methods: Theory and Applications*, pages 1–32, 1984.

[10] M. Bucelli, M. Geraint Gabriel, A. Quarteroni, G. Gigante, and C. Vergara. A stable loosely-coupled scheme for cardiac electro-fluid-structure interaction. *Journal of Computational Physics*, 490:112326, 2023.

[11] M. Bukač, A. Seboldt, and C. Trenchea. Refactorization of Cauchy’s method: A second-order partitioned method for fluid-thick structure interaction problems. *Journal of Mathematical Fluid Mechanics*, 23(3):1–25, 2021.

[12] M. Bukač, S. Čanić, R. Glowinski, B. Muha, and A. Quaini. A modular, operator-splitting scheme for fluid-structure interaction problems with thick structures. *Internat. J. Numer. Methods Fluids*, 74(8):577–604, 2014.

[13] M. Bukač. An extension of explicit coupling for fluid-structure interaction problems. *Mathematics*, 9(15), 2021.

[14] E. Burman, R. Durst, M. Fernández, and J. Guzmán. Loosely coupled, non-iterative time-splitting scheme based on Robin-Robin coupling: Unified analysis for parabolic/parabolic and parabolic/hyperbolic problems. *Journal of Numerical Mathematics*, 31(1):59–77, 2023.

[15] E. Burman, R. Durst, M. Fernández, J. Guzmán, and S. Liu. Estimates of discrete time derivatives for the parabolic-parabolic Robin-Robin coupling methods. *Preprint*, 2024.
define $\mathfrak{Z}$ as,

$$ \mathfrak{Z} := \frac{1}{\nu_s} \|\partial_t^4 \mathbf{w}\|_{L^2(0,T;L^2(\Omega_s))}^2 + \left(\frac{1}{\nu_f} + 1\right) \|\partial_t^4 \mathbf{u}\|_{L^2(0,T;L^2(\Omega_f))}^2 \\ + \left(\frac{\nu_f}{\alpha^2} + \frac{\nu_f^2}{\alpha}\right) \|\partial_t^3 \mathbf{u}\|_{L^2(0,T;H^1(\Omega_f))}^2 + \frac{(\nu_f)^3}{\alpha^2} \|\partial_t^3 \mathbf{u}\|_{L^2(0,T;H^2(\Omega_f))}^2 \\ + \frac{\alpha^2}{\nu_s} \|\partial_t^3 \mathbf{u}\|_{L^2(0,T;L^2(\Sigma))}^2 + \frac{1}{\nu_f} \|\partial_t^3 \mathbf{l}\|_{L^2(0,T;L^2(\Sigma))}^2 + \frac{\nu_f^2}{\alpha^2} \|\partial_t^3 \mathbf{u}\|_{L^\infty(0,T;H^1(\Omega_f))}^2, $$

define $Y$ as,

$$ Y := \frac{1}{\nu_s} \|\partial_x \partial_t^2 \mathbf{w}\|_{L^2(0,T;L^2(\Omega_s))}^2 + \left(\frac{1}{\nu_f} + 1\right) \|\partial_x \partial_t^2 \mathbf{u}\|_{L^2(0,T;L^2(\Omega_f))}^2 \\ + \left(\frac{\nu_f}{\alpha^2} + \frac{\nu_f^2}{\alpha}\right) \|\partial_x \partial_t \mathbf{u}\|_{L^2(0,T;H^1(\Omega_f))}^2 + \frac{(\nu_f)^3}{\alpha^2} \|\partial_x \partial_t \mathbf{u}\|_{L^2(0,T;H^2(\Omega_f))}^2 \\ + \frac{\alpha^2}{\nu_s} \|\partial_x \partial_t \mathbf{u}\|_{L^2(0,T;L^2(\Sigma))}^2 + \frac{1}{\nu_f} \|\partial_x \partial_t \mathbf{l}\|_{L^2(0,T;L^2(\Sigma))}^2 + \frac{\nu_f^2}{\alpha^2} \|\partial_x \partial_t \mathbf{u}\|_{L^\infty(0,T;H^1(\Omega_f))}^2, $$

and define $\mathfrak{Y}$ as,

$$ \mathfrak{Y} := \frac{1}{\nu_s} \|\partial_x \partial_t^3 \mathbf{w}\|_{L^2(0,T;L^2(\Omega_s))}^2 + \left(\frac{1}{\nu_f} + 1\right) \|\partial_x \partial_t^3 \mathbf{u}\|_{L^2(0,T;L^2(\Omega_f))}^2 \\ + \left(\frac{\nu_f}{\alpha^2} + \frac{\nu_f^2}{\alpha}\right) \|\partial_x \partial_t^2 \mathbf{u}\|_{L^2(0,T;H^1(\Omega_f))}^2 + \frac{(\nu_f)^3}{\alpha^2} \|\partial_x \partial_t^2 \mathbf{u}\|_{L^2(0,T;H^2(\Omega_f))}^2 \\ + \frac{\alpha^2}{\nu_s} \|\partial_x \partial_t^2 \mathbf{u}\|_{L^2(0,T;L^2(\Sigma))}^2 + \frac{1}{\nu_f} \|\partial_x \partial_t^2 \mathbf{u}\|_{L^2(0,T;L^2(\Sigma))}^2 + \frac{\nu_f^2}{\alpha^2} \|\partial_x \partial_t^2 \mathbf{u}\|_{L^\infty(0,T;H^1(\Omega_f))}^2. $$

## REFERENCES

[1] M. Aggul, J. M. Connors, D. Erkmen, and A. E. Labovsky. A defect-deferred correction method for fluid-fluid interaction. *SIAM Journal on Numerical Analysis*, 56(4):2484–2512, 2018.

[2] M. Alnaes, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M. E. Rognes, and G. N. Wells. The FEiNCS project version 1.5. *Archive of numerical software*, 3(100), 2015.

[3] S. Badia, F. Nobile, and C. Vergara. Fluid-structure partitioned procedures based on Robin transmission conditions. *Journal of Computational Physics*, 227(14):7027–7051, 2008.

[4] F. Ballarin and G. Rozza. https://multiphenics.github.io/, 2016.

[5] J. W. Banks, W. D. Henshaw, and D. W. Schwendeman. An analysis of a new stable partitioned algorithm for FSI problems. Part I: Incompressible flow and elastic solids. *J. Comput. Phys.*, 269:108–137, 2014.

[6] M. Beneš. Convergence and stability analysis of heterogeneous time step coupling schemes for parabolic problems. *Appl. Numer. Math.*, 121:198–222, 2017.

[7] M. Beneš, A. Nekvinda, and M. K. Yadav. Multi-time-step domain decomposition method with non-matching grids for parabolic problems. *Appl. Math. Comput.*, 267:571–582, 2015.

[8] K. Böhmer. Discrete Newton methods and iterated defect corrections. *Numer. Math.*, 37(2):167–192, 1981.

[9] K. Böhmer, P. W. Hemker, and H. J. Stetter. The defect correction approach. *Defect Correction Methods: Theory and Applications*, pages 1–32, 1984.

[10] M. Bucelli, M. Geraint Gabriel, A. Quarteroni, G. Gigante, and C. Vergara. A stable loosely-coupled scheme for cardiac electro-fluid-structure interaction. *Journal of Computational Physics*, 490:112326, 2023.

[11] M. Bukač, A. Seboldt, and C. Trenchea. Refactorization of Cauchy’s method: A second-order partitioned method for fluid-thick structure interaction problems. *Journal of Mathematical Fluid Mechanics*, 23(3):1–25, 2021.

[12] M. Bukač, S. Čanić, R. Glowinski, B. Muha, and A. Quaini. A modular, operator-splitting scheme for fluid-structure interaction problems with thick structures. *Internat. J. Numer. Methods Fluids*, 74(8):577–604, 2014.

[13] M. Bukač. An extension of explicit coupling for fluid-structure interaction problems. *Mathematics*, 9(15), 2021.

[14] E. Burman, R. Durst, M. Fernández, and J. Guzmán. Loosely coupled, non-iterative time-splitting scheme based on Robin-Robin coupling: Unified analysis for parabolic/parabolic and parabolic/hyperbolic problems. *Journal of Numerical Mathematics*, 31(1):59–77, 2023.

[15] E. Burman, R. Durst, M. Fernández, J. Guzmán, and S. Liu. Estimates of discrete time derivatives for the parabolic-parabolic Robin-Robin coupling methods. *Preprint*, 2024.
[16] E. Burman, R. Durst, M. A. Fernández, and J. Guzmán. Fully discrete loosely coupled Robin-Robin scheme for incompressible fluid-structure interaction: stability and error analysis. *Numerische Mathematik*, 151(4):807-840, 2022.

[17] E. Burman, R. Durst, M. A. Fernández, J. Guzmán, and O. Ruz. Robin-Robin loose coupling for incompressible fluid-structure interaction: non-linear setting and nearly-optimal error analysis. 2023. https://inria.hal.science/hal-04258861.

[18] E. Burman, R. Durst, and J. Guzmán. Stability and error analysis of a splitting method using Robin-Robin coupling applied to a fluid-structure interaction problem. *Numerical Methods for Partial Differential Equations*, 38(5):1396-1406, 2022.

[19] E. Burman and M. A. Fernández. Stabilization of explicit coupling in fluid-structure interaction involving fluid incompressibility. *Computer Methods in Applied Mechanics and Engineering*, 198(5-8):766-784, 2009.

[20] E. Burman and M. A. Fernández. Explicit strategies for incompressible fluid-structure interaction problems: Nitsche type mortaring versus Robin-Robin coupling. *International Journal for Numerical Methods in Engineering*, 97(10):739-758, 2014.

[21] C. Canuto and A. Lo Giudice. A multi-timestep Robin-Robin domain decomposition method for time dependent advection-diffusion problems. *Appl. Math. Comput.*, 363:124596, 14, 2019.

[22] J. M. Connors and J. S. Howell. A fluid-fluid interaction method using decoupled subproblems and differing time steps. *Numer. Methods Partial Differential Equations*, 28(4):1283-1308, 2012.

[23] J. M. Connors, J. S. Howell, and W. J. Layton. Partitioned time stepping for a parabolic two domain problem. *SIAM Journal on Numerical Analysis*, 47(5):3526-3549, 2009.

[24] J. M. Connors, J. S. Howell, and W. J. Layton. Decoupled time stepping methods for fluid-fluid interaction. *SIAM J. Numer. Anal.*, 50(3):1297-1319, 2012.

[25] R. Durst. *Recent Advances in Splitting Methods Based on Robin-Robin Coupling Conditions*. PhD thesis, Brown University, 2022.

[26] M. Fernández and J. Mullaert. Convergence and error analysis for a class of splitting schemes in incompressible fluid-structure interaction. *IMA J. Numer. Anal.*, 36(4):1748-1782, 2016.

[27] G. Gigante and C. Vergara. On the choice of interface parameters in Robin-Robin loosely coupled schemes for fluid-structure interaction. *Fluids*, 6(6), 2021.

[28] G. Gigante and C. Vergara. On the stability of a loosely-coupled scheme based on a Robin interface condition for fluid-structure interaction. *Comput. Math. Appl.*, 96:109-119, 2021.

[29] W. Li, P. Huang, and Y. He. An unconditionally energy stable finite element scheme for a nonlinear fluid-fluid interaction model. *IMA J. Numer. Anal.*, 44(1):157-191, 2024.

[30] D. A. Serino, J. W. Banks, W. D. Henshaw, and D. W. Schwendeman. A stable added-mass partitioned (amp) algorithm for elastic solids and incompressible flow: Model problem analysis. *SIAM Journal on Scientific Computing*, 41(4):A2464-A2484, 2019.

[31] K. C. Sockwell, P. Bochev, K. Peterson, and P. Kuberry. Interface flux recovery framework for constructing partitioned heterogeneous time-integration methods. *Numer. Methods Partial Differential Equations*, 39(5):3572-3593, 2023.

[32] H. J. Stetter. The defect correction principle and discretization methods. *Numer. Math.*, 29(4):425-443, 1977/78.

[33] H. Zhang, Z. Liu, E. Constantinescu, and R. Jacob. Stability analysis of interface conditions for ocean-atmosphere coupling. *J. Sci. Comput.*, 84(3):Paper No. 44, 25, 2020.

[34] Y. Zhang, Y. Hou, and L. Shan. Error estimates of a decoupled algorithm for a fluid-fluid interaction problem. *J. Comput. Appl. Math.*, 333:266-291, 2018.

[35] Y. Zhang, L. Shan, and Y. Hou. New approach to prove the stability of a decoupled algorithm for a fluid-fluid interaction problem. *J. Comput. Appl. Math.*, 371:112695, 19, 2020.

<footnote n="1">¹DEPARTMENT OF MATHEMATICS, UNIVERSITY COLLEGE LONDON, LONDON, UK-WC1E 6BT, UNITED KINGDOM
*Email address:* e.burman@ucl.ac.uk</footnote>

<footnote n="2">²DEPARTMENT OF MATHEMATICS, UNIVERSITY OF PITTSBURGH, PITTSBURGH, PENNSYLVANIA 15261, USA
*Email address:* rebecca.durst@alumni.brown.edu</footnote>

<footnote n="3">³SORBONNE UNIVERSITÉ & CNRS, UMR 7598 LJLL, 75005 PARIS FRANCE - INRIA, 75012 PARIS, FRANCE
*Email address:* miguel.fernandez@inria.fr</footnote>

<footnote n="4">⁴BROWN UNIVERSITY
*Email address:* johnny_guzman@brown.edu</footnote>

<footnote n="5">⁵THE INSTITUTE FOR COMPUTATIONAL AND EXPERIMENTAL RESEARCH IN MATHEMATICS, BROWN UNIVERSITY, PROVIDENCE, RI 02912, USA
*Email address:* sijing.liu@brown.edu</footnote>
[16] E. Burman, R. Durst, M. A. Fernández, and J. Guzmán. Fully discrete loosely coupled Robin-Robin scheme for incompressible fluid-structure interaction: stability and error analysis. *Numerische Mathematik*, 151(4):807-840, 2022.

[17] E. Burman, R. Durst, M. A. Fernández, J. Guzmán, and O. Ruz. Robin-Robin loose coupling for incompressible fluid-structure interaction: non-linear setting and nearly-optimal error analysis. 2023. https://inria.hal.science/hal-04258861.

[18] E. Burman, R. Durst, and J. Guzmán. Stability and error analysis of a splitting method using Robin-Robin coupling applied to a fluid-structure interaction problem. *Numerical Methods for Partial Differential Equations*, 38(5):1396-1406, 2022.

[19] E. Burman and M. A. Fernández. Stabilization of explicit coupling in fluid-structure interaction involving fluid incompressibility. *Computer Methods in Applied Mechanics and Engineering*, 198(5-8):766-784, 2009.

[20] E. Burman and M. A. Fernández. Explicit strategies for incompressible fluid-structure interaction problems: Nitsche type mortaring versus Robin-Robin coupling. *International Journal for Numerical Methods in Engineering*, 97(10):739-758, 2014.

[21] C. Canuto and A. Lo Giudice. A multi-timestep Robin-Robin domain decomposition method for time dependent advection-diffusion problems. *Appl. Math. Comput.*, 363:124596, 14, 2019.

[22] J. M. Connors and J. S. Howell. A fluid-fluid interaction method using decoupled subproblems and differing time steps. *Numer. Methods Partial Differential Equations*, 28(4):1283-1308, 2012.

[23] J. M. Connors, J. S. Howell, and W. J. Layton. Partitioned time stepping for a parabolic two domain problem. *SIAM Journal on Numerical Analysis*, 47(5):3526-3549, 2009.

[24] J. M. Connors, J. S. Howell, and W. J. Layton. Decoupled time stepping methods for fluid-fluid interaction. *SIAM J. Numer. Anal.*, 50(3):1297-1319, 2012.

[25] R. Durst. *Recent Advances in Splitting Methods Based on Robin-Robin Coupling Conditions*. PhD thesis, Brown University, 2022.

[26] M. Fernández and J. Mullaert. Convergence and error analysis for a class of splitting schemes in incompressible fluid-structure interaction. *IMA J. Numer. Anal.*, 36(4):1748-1782, 2016.

[27] G. Gigante and C. Vergara. On the choice of interface parameters in Robin-Robin loosely coupled schemes for fluid-structure interaction. *Fluids*, 6(6), 2021.

[28] G. Gigante and C. Vergara. On the stability of a loosely-coupled scheme based on a Robin interface condition for fluid-structure interaction. *Comput. Math. Appl.*, 96:109-119, 2021.

[29] W. Li, P. Huang, and Y. He. An unconditionally energy stable finite element scheme for a nonlinear fluid-fluid interaction model. *IMA J. Numer. Anal.*, 44(1):157-191, 2024.

[30] D. A. Serino, J. W. Banks, W. D. Henshaw, and D. W. Schwendeman. A stable added-mass partitioned (amp) algorithm for elastic solids and incompressible flow: Model problem analysis. *SIAM Journal on Scientific Computing*, 41(4):A2464-A2484, 2019.

[31] K. C. Sockwell, P. Bochev, K. Peterson, and P. Kuberry. Interface flux recovery framework for constructing partitioned heterogeneous time-integration methods. *Numer. Methods Partial Differential Equations*, 39(5):3572-3593, 2023.

[32] H. J. Stetter. The defect correction principle and discretization methods. *Numer. Math.*, 29(4):425-443, 1977/78.

[33] H. Zhang, Z. Liu, E. Constantinescu, and R. Jacob. Stability analysis of interface conditions for ocean-atmosphere coupling. *J. Sci. Comput.*, 84(3):Paper No. 44, 25, 2020.

[34] Y. Zhang, Y. Hou, and L. Shan. Error estimates of a decoupled algorithm for a fluid-fluid interaction problem. *J. Comput. Appl. Math.*, 333:266-291, 2018.

[35] Y. Zhang, L. Shan, and Y. Hou. New approach to prove the stability of a decoupled algorithm for a fluid-fluid interaction problem. *J. Comput. Appl. Math.*, 371:112695, 19, 2020.

<footnote n="1">¹DEPARTMENT OF MATHEMATICS, UNIVERSITY COLLEGE LONDON, LONDON, UK-WC1E 6BT, UNITED KINGDOM
*Email address:* e.burman@ucl.ac.uk</footnote>

<footnote n="2">²DEPARTMENT OF MATHEMATICS, UNIVERSITY OF PITTSBURGH, PITTSBURGH, PENNSYLVANIA 15261, USA
*Email address:* rebecca.durst@alumni.brown.edu</footnote>

<footnote n="3">³SORBONNE UNIVERSITÉ & CNRS, UMR 7598 LJLL, 75005 PARIS FRANCE - INRIA, 75012 PARIS, FRANCE
*Email address:* miguel.fernandez@inria.fr</footnote>

<footnote n="4">⁴BROWN UNIVERSITY
*Email address:* johnny_guzman@brown.edu</footnote>

<footnote n="5">⁵THE INSTITUTE FOR COMPUTATIONAL AND EXPERIMENTAL RESEARCH IN MATHEMATICS, BROWN UNIVERSITY, PROVIDENCE, RI 02912, USA
*Email address:* sijing.liu@brown.edu</footnote>