# A space-time multiscale mortar mixed finite element method for parabolic equations

Manu Jayadharan, Michel Kern, Martin Vohralík, Ivan Yotov

► To cite this version:

Manu Jayadharan, Michel Kern, Martin Vohralík, Ivan Yotov. A space-time multiscale mortar mixed finite element method for parabolic equations. SIAM Journal on Numerical Analysis, 2023, 61 (2), pp.675 - 706. 10.1137/21M1447945. hal-03355088

HAL Id: hal-03355088

https://inria.hal.science/hal-03355088

Submitted on 27 Sep 2021

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 space-time multiscale mortar mixed finite element method for parabolic equations*

Manu Jayadharan†

Michel Kern‡§

Martin Vohralík‡§

Ivan Yotov†

September 21, 2021

## Abstract

We develop a space-time mortar mixed finite element method for parabolic problems. The domain is decomposed into a union of subdomains discretized with non-matching spatial grids and asynchronous time steps. The method is based on a space-time variational formulation that couples mixed finite elements in space with discontinuous Galerkin in time. Continuity of flux (mass conservation) across space-time interfaces is imposed via a coarse-scale space-time mortar variable that approximates the primary variable. Uniqueness, existence, and stability, as well as a priori error estimates for the spatial and temporal errors are established. A space-time non-overlapping domain decomposition method is developed that reduces the global problem to a space-time coarse-scale mortar interface problem. Each interface iteration involves solving in parallel space-time subdomain problems. The spectral properties of the interface operator and the convergence of the interface iteration are analyzed. Numerical experiments are provided that illustrate the theoretical results and the flexibility of the method for modeling problems with features that are localized in space and time.

## 1 Introduction

The multiscale mortar mixed finite element method of [3,5] allows for highly efficient and accurate discretization of elliptic problems. Let a spatial domain $\Omega$ be given, which is decomposed into subdomains $\Omega_i$. Then, on each $\Omega_i$, an individual mesh is set, and a standard mixed finite element scheme is considered. A stand-alone mortar variable approximating the primary variable is further introduced on an independent interface mesh, which is typically coarser but where one possibly employs polynomials of higher degree. It is used to couple the subdomain problems and to ensure (multiscale) weak continuity of the normal component of the mixed finite element flux variable over the interfaces between subdomains. Moreover, the mortar variable enables a very efficient parallelization in space via a non-overlapping domain decomposition algorithm based on reduction to an interface problem [3, 5, 28].

We introduce here a *space-time discretization* of a model parabolic equation, extending the above philosophy to time-dependent problems. Let a time interval $(0, T)$ be given. For each subdomain

<footnote>†Department of Mathematics, University of Pittsburgh, Pittsburgh, PA 15260, USA (manu.jayadharan@pitt.edu, yotov@math.pitt.edu).</footnote>

<footnote>‡Inria, 2 rue Simone Iff, 75589 Paris, France (michel.kern@inria.fr, martin.vohralik@inria.fr).</footnote>

<footnote>§CERMICS, Ecole des Ponts, 77455 Marne-la-Valle, France.</footnote>

<footnote>*This project has received partial funding from the US National Science Foundation grants DMS 1818775 and DMS 2111129 and the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No 647134 GATIPOR).</footnote>
$\Omega_i$, our approach considers an individual space mesh of $\Omega_i$ along with individual time stepping on $(0, T)$. On each *space-time subdomain* $\Omega_i \times (0, T)$, any standard mixed finite element scheme is combined with the discontinuous Galerkin (DG) time discretization [50]. Then a stand-alone mortar variable approximating the primary variable is introduced on an independent *space-time interface mesh*, which is typically coarse and where higher polynomial degrees may be used. It is used to couple the space-time subdomain problems and to ensure (multiscale) weak continuity of the normal component of the mixed finite element flux variable (and consequently mass conservation) over the space-time interfaces. This setting allows for high flexibility with individual discretizations of each space-time subdomain $\Omega_i \times (0, T)$, and in particular for *local time stepping*. Moreover, *space-time parallelization* can be achieved via reduction to a space-time interface problem requiring the solution of discrete problems on the individual space-time subdomains $\Omega_i \times (0, T)$, exchanging space-time boundary data through transmission conditions, in the spirit of space-time domain decomposition methods as in, e.g., [10, 21, 23-25, 32-34, 45, 51].

We mention some of the related previous works on local time stepping for parabolic problems in mixed formulations. The early work [19] studies finite difference methods on grids with local refinement in space and time. A similar approach is employed in [16] for transport equations. Two space-time domain decomposition methods are considered in [33] – a space-time Steklov-Poincaré operator and optimized Schwarz waveform relaxation (OSWR) [23] with Robin transmission conditions. Asynchronous time stepping is allowed, but the spatial grids are assumed matching. The focus is on the analysis of the iterative convergence of the OSWR method. A posteriori error estimates for these methods are developed in [2, 4] for nested time grids, with [4] also allowing for non-matching spatial grids through the use of mortar finite elements. The methods from [33] are extended to fracture modeling in [34]. Overlapping Schwarz domain decomposition with local grid refinement in space and time for two-phase flow in porous media is developed in [39]. Domain decomposition methods for mortar mixed finite element methods for parabolic problems with non-matching spatial grids and uniform time stepping are studied in [6, 22]. For parabolic problems in non-mixed form, domain decomposition methods with local time stepping have been studied, e.g., in [10, 15, 21, 24, 31, 32, 41, 45, 51]. Parallelism in time has also been explored, such as the Parareal algorithm [27, 42] and multigrid in time [20, 26]. Local time stepping techniques have been developed for multiphysics systems coupled through interface conditions, e.g., for the Stokes-Darcy system [35, 47]. Finally, we mention some earlier works on space-time methods for parabolic problems coupling mixed finite element discretizations in space with DG in time on a single domain. In [8], a method using continuous trial and discontinuous test functions in time is developed. A posteriori error estimation and space-time adaptivity for mixed finite element - DG methods is studied in [12, 40].

To the best of the authors' knowledge, the solvability, stability, and a priori error analysis for space-time domain decomposition methods with non-matching spatial grids and asynchronous time stepping have not been studied in the literature, which is the main goal of this paper. A key tool in the analysis is the construction of an interpolant in a space-time weakly continuous velocity space, which is used to prove a discrete divergence inf-sup condition on this space. Another key component is establishing a discrete space-time mortar inf-sup condition under a suitable assumption on the mortar space. In addition to performing complete analysis in the general case, we also consider conforming time discretizations. In this case, we provide stability and error bounds for the velocity divergence and improved error estimates. To the best of our knowledge, such result has not been established in the literature for space-time mixed finite element methods with a DG time discretization, even on a single domain. Finally, we develop a parallel non-overlapping domain decomposition algorithm for the solution of the resulting algebraic problem. In particular, we utilize a time-dependent Steklov-Poincaré operator approach to reduce the global problem to a space-time
interface problem. We show that the interface operator is positive definite and analyze its spectral properties. We employ an interface GMRES algorithm, which involves solving in parallel space-time subdomain problems at each iteration. The current iteration value of the mortar variable provides a Dirichlet boundary data on the space-time interfaces for the subdomain problems. We emphasize that, due to the discontinuous time discretization, the space-time subdomain problems are solved using classical time marching over the local time grid. We utilize the spectral bound of the interface operator to obtain an estimate for the number of interface GMRES iterations through field-of-values analysis.

This contribution is organized as follows. In Section 2, we describe the model problem and its domain decomposition weak formulation. Our space-time multiscale mortar discretization is introduced in Section 3, and we prove its existence, uniqueness, and stability with respect to data in Section 4. Section 5 derives a priori error estimates. Control of the velocity divergence and improved error estimates are established in Section 6. The reduction to a space-time interface problem and its analysis are presented in Section 7. We finally present numerical illustrations in Section 8 and close with conclusions in Section 9.

## 2 Setting

In this section we introduce the setting for our study.

### 2.1 Model problem

We consider a parabolic partial differential equation in a mixed form, modeling single phase flow in porous media. Let $\Omega \subset \mathbb{R}^d$, $d = 2, 3$, be a spatial polytopal domain with Lipschitz boundary and let $T > 0$ be the final time. The governing equations are

$$ \mathbf{u} = -K \nabla p, \quad \frac{\partial p}{\partial t} + \nabla \cdot \mathbf{u} = q \quad \text{in } \Omega \times (0, T], \qquad (2.1a) $$

where $p$ is the fluid pressure, $\mathbf{u}$ is the Darcy velocity, $q$ is a source term, and $K$ is a tensor representing the rock permeability divided by the fluid viscosity. We assume for simplicity homogeneous Dirichlet boundary condition

$$ p(x, t) = 0 \quad \text{on} \quad \partial\Omega \times (0, T] \qquad (2.1b) $$

and assign the initial pressure

$$ p(x, 0) = p_0(x) \quad \text{on} \quad \Omega. \qquad (2.1c) $$

We assume that $q \in L^2(0, T; L^2(\Omega))$, $p_0 \in H_0^1(\Omega)$, $\nabla \cdot K \nabla p_0 \in L^2(\Omega)$, and that $K$ is a spatially-dependent, uniformly bounded, symmetric, and positive definite tensor, i.e., for constants $0 < k_{\min} \le k_{\max} < \infty$,

$$ \text{for a.e. } x \in \Omega, \quad k_{\min} \zeta^T \zeta \le \zeta^T K(x) \zeta \le k_{\max} \zeta^T \zeta \quad \forall \zeta \in \mathbb{R}^d, \ d = 2, 3. \qquad (2.2) $$

Moreover, we suppose a scaling such that the diameter of $\Omega$ and the final time $T$ are of order one.

### 2.2 Space-time subdomains

Let $\Omega$ be a union of non-overlapping polytopal subdomains with Lipschitz boundary, $\bar{\Omega} = \cup \bar{\Omega}_i$. Let $\Gamma_i = \partial \Omega_i \setminus \partial \Omega$ be the interior boundary of $\Omega_i$, let $\Gamma_{ij} = \Gamma_i \cap \Gamma_j$ be the interface between two adjacent subdomains $\Omega_i$ and $\Omega_j$, and let $\Gamma = \cup \Gamma_{ij}$ be the union of all subdomain interfaces. We also introduce the space-time counterparts $\Omega^T = \Omega \times (0, T)$, $\Omega_i^T = \Omega_i \times (0, T)$, $\Gamma_i^T = \Gamma_i \times (0, T)$,
and $\Gamma_{ij}^T = \Gamma_{ij} \times (0, T)$. We will introduce space-time domain decomposition discretizations based on $\Omega_i^T$.

## 2.3 Basic notation

We will utilize the following notation. For a domain $\mathcal{O} \subset \mathbb{R}^d$, the $L^2(\mathcal{O})$ inner product and norm for scalar and vector-valued functions are denoted by $(\cdot, \cdot)_\mathcal{O}$ and $\| \cdot \|_\mathcal{O}$, respectively. The norms and seminorms of the Sobolev spaces $W^{k,p}(\mathcal{O})$, $k \in \mathbb{R}$, $p \ge 1$, are denoted by $\| \cdot \|_{k,p,\mathcal{O}}$ and $| \cdot |_{k,p,\mathcal{O}}$, respectively. The norms and seminorms of the Hilbert spaces $H^k(\mathcal{O})$ are denoted by $\| \cdot \|_{k,\mathcal{O}}$ and $| \cdot |_{k,\mathcal{O}}$, respectively. For a section of a subdomain boundary $S \subset \mathbb{R}^{d-1}$ we write $\langle \cdot, \cdot \rangle_S$ and $\| \cdot \|_S$ for the $L^2(S)$ inner product (or duality pairing) and norm, respectively. By $\mathbf{M}$ we denote the vectorial counterpart of a generic scalar space $M$.

The above notation is extended to space-time domains as follows. For $\mathcal{O}^T = \mathcal{O} \times (0, T)$ and $S^T = S \times (0, T)$, let $(\cdot, \cdot)_{\mathcal{O}^T} = \int_0^T (\cdot, \cdot)_{\mathcal{O}}$ and $\langle \cdot, \cdot \rangle_{S^T} = \int_0^T \langle \cdot, \cdot \rangle_S$. For space-time norms we use the standard Bochner notation. For example, given a spatial norm $\| \cdot \|_V$, we denote, for $p > 0$,

$$ \| \cdot \|_{L^p(0,T;V)} = \left( \int_0^T \| \cdot \|_V^p \right)^{\frac{1}{p}}, \quad \| \cdot \|_{L^\infty(0,T;V)} = \operatorname{ess\,sup} \| \cdot \|_V, $$

with the usual extension for $\| \cdot \|_{W^{k,p}(0,T;V)}$ and $\| \cdot \|_{H^k(0,T;V)}$. Let $\| \cdot \|_{S^T} = \| \cdot \|_{L^2(0,T;L^2(S))}$. We will use the space

$$ \mathbf{H}(\mathrm{div}; \mathcal{O}) = \{ \mathbf{v} \in \mathbf{L}^2(\mathcal{O}) : \nabla \cdot \mathbf{v} \in \mathbf{L}^2(\mathcal{O}) \}, $$

equipped with the norm

$$ \|\mathbf{v}\|_{\mathrm{div};\mathcal{O}} = (\|\mathbf{v}\|_{\mathcal{O}}^2 + \|\nabla \cdot \mathbf{v}\|_{\mathcal{O}}^2)^{\frac{1}{2}}. $$

Finally, throughout the paper, $C$ will denote a generic constant that is independent of the spatial and temporal discretization parameters.

## 2.4 Weak formulation

The weak formulation of problem (2.1) reads: find $(\mathbf{u}, p) : [0, T] \mapsto \mathbf{H}(\mathrm{div}; \mathcal{O}) \times L^2(\mathcal{O})$ such that $p(x, 0) = p_0$ and for a.e. $t \in (0, T)$,

$$ (K^{-1}\mathbf{u}, \mathbf{v})_{\mathcal{O}} - (p, \nabla \cdot \mathbf{v})_{\mathcal{O}} = 0 \quad \forall \mathbf{v} \in \mathbf{H}(\mathrm{div}; \mathcal{O}), \qquad (2.3a) $$

$$ (\partial_t p, w)_{\mathcal{O}} + (\nabla \cdot \mathbf{u}, w)_{\mathcal{O}} = (q, w)_{\mathcal{O}} \quad \forall w \in L^2(\mathcal{O}). \qquad (2.3b) $$

The following well-posedness result is rather standard and presented in, e.g., [33, Theorem 2.1].

**Theorem 2.1 (Well-posedness).** *Problem (2.3) has a unique solution $\mathbf{u} \in L^2(0, T; \mathbf{H}(\mathrm{div}; \mathcal{O})) \cap L^\infty(0, T; \mathbf{L}^2(\mathcal{O}))$, $p \in H^1(0, T; H_0^1(\mathcal{O}))$.*

We note that in particular the inclusion $p \in H^1(0, T; H_0^1(\mathcal{O}))$ follows from (2.3a), which implies that for a.e. $t \in (0, T)$, $\nabla p = -K^{-1}\mathbf{u}$ in the sense of distributions.

## 2.5 Domain decomposition weak formulation

We now give a domain decomposition weak formulation of (2.3). Introduce the subdomain velocity and pressure spaces

$$ \mathbf{V}_i = \mathbf{H}(\mathrm{div}; \Omega_i), \quad \mathbf{V} = \bigoplus \mathbf{V}_i, \quad W_i = L^2(\Omega_i), \quad W = \bigoplus W_i = L^2(\Omega), $$
endowed with the norms

$$ \|\mathbf{v}\|_{\mathbf{V}_i} = \|\mathbf{v}\|_{\text{div};\Omega_i}, \quad \|\mathbf{v}\|_{\mathbf{V}} = \left(\sum_i \|\mathbf{v}\|_{\mathbf{V}_i}^2\right)^{\frac{1}{2}}, \quad \|w\|_{W} = \|w\|_{\Omega}. $$

We also introduce the following spatial bilinear forms, which will prove useful below:

$$ a_i(\mathbf{u}, \mathbf{v}) = (K^{-1}\mathbf{u}, \mathbf{v})_{\Omega_i}, \quad a(\mathbf{u}, \mathbf{v}) = \sum_i a_i(\mathbf{u}, \mathbf{v}), \qquad (2.4a) $$

$$ b_i(\mathbf{v}, w) = -(\nabla \cdot \mathbf{v}, w)_{\Omega_i}, \quad b(\mathbf{v}, w) = \sum_i b_i(\mathbf{v}, w), \qquad (2.4b) $$

$$ b_{\Gamma}(\mathbf{v}, \mu) = \sum_{i} \langle \mathbf{v} \cdot \mathbf{n}_{i}, \mu \rangle_{\Gamma_{i}}. \qquad (2.4c) $$

In addition, for any spatial bilinear form $s(\cdot, \cdot)$, let $s^T(\cdot, \cdot) = \int_0^T s(\cdot, \cdot)$.

Now, since $p \in H^1(0, T; H_0^1(\Omega))$, we can consider the trace of the pressure $p$ on the interfaces, $\lambda = p|_{\Gamma}$. Thus, integrating in time, it is easy to see that the solution $(\mathbf{u}, p)$ of (2.3) satisfies

$$ a^T(\mathbf{u}, \mathbf{v}) + b^T(\mathbf{v}, p) + b_{\Gamma}^T(\mathbf{v}, \lambda) = 0 \quad \forall \mathbf{v} \in L^2(0, T; \mathbf{V}), \qquad (2.5a) $$

$$ (\partial_t p, w)_{\Omega^T} - b^T(\mathbf{u}, w) = (q, w)_{\Omega^T} \quad \forall w \in L^2(0, T; W). \qquad (2.5b) $$

# 3 Space-time mortar mixed finite element method

We consider a space-time discretization of (2.5), motivated by [33]. It employs a mortar finite element variable to approximate the pressure trace $\lambda$ from (2.5) and uses it as a Lagrange multiplier to impose weakly the continuity of flux across space-time interfaces.

## 3.1 Space-time grids and spaces

Let $\mathcal{T}_{h,i}$ be a shape-regular partition of the subdomain $\Omega_i$ into parallelepipeds or simplices in the sense of [13]. We stress that this allows for grids that do not match along the interfaces $\Gamma_{ij}$ between subdomains $\Omega_i$ and $\Omega_j$. Let $h_i = \max_{E \in \mathcal{T}_{h,i}} \text{diam } E$ and $h = \max_i h_i$. In some parts of the analysis, we will additionally require $\mathcal{T}_{h,i}$ to be quasi-uniform in that $h_i \le C \text{diam } E \ \forall E \in \mathcal{T}_{h,i}$, as well as that the mesh sizes in all subdomains be comparable in that $h \le C h_i \ \forall i$. Similarly, let $\mathcal{T}_{i}^{\Delta t} : 0 = t_i^0 < t_i^1 < \dots < t_i^{N_i} = T$ be a partition of the time interval $(0, T)$ corresponding to subdomain $\Omega_i$. This means that we consider different time discretizations on different subdomains. Let $\Delta t_i = \max_{1 \le k \le N_i} |t_i^k - t_i^{k-1}|$ and $\Delta t = \max_i \Delta t_i$. Though we admit non-uniform time stepping, we will sometimes require that $\mathcal{T}_{i}^{\Delta t}$ be quasi-uniform in that $\Delta t_i \le C |t_i^k - t_i^{k-1}| \ \forall (t_i^{k-1}, t_i^k) \in \mathcal{T}_{i}^{\Delta t}$, as well as that the time steps in all subdomains be comparable in that $\Delta t \le C \Delta t_i \ \forall i$. Composing $\mathcal{T}_{h,i}$ and $\mathcal{T}_{i}^{\Delta t}$ by tensor product results in a space-time partition

$$ \mathcal{T}_{h,i}^{\Delta t} = \mathcal{T}_{h,i} \times \mathcal{T}_{i}^{\Delta t} $$

of the space-time subdomain $\Omega_i^T$. An illustration is given in Figure 1, where yet a different, mortar space-time grid, is also shown in the middle.

For discretization in space, we consider any of the inf-sup stable mixed finite element spaces $\mathbf{V}_{h,i} \times W_{h,i} \subset \mathbf{V}_i \times W_i$ such as the Raviart-Thomas or the Brezzi-Douglas-Marini spaces, see,
<figure><img src="image_2.png" /><figcaption>Figure 1: Non-matching space-time subdomain and mortar grids in two spatial dimensions.</figcaption></figure>

e.g., [11]. For discretization in time, we will in turn utilize the discontinuous Galerkin (DG) method, cf. [50], which is based on a discontinuous piecewise polynomial approximation of the solution on the mesh $\mathcal{T}_i^{\Delta t}$. Denote by $W_i^{\Delta t}$ the subdomain time discretizations of the velocity and pressure. Composing the space and time discretizations

$$ \mathbf{V}_{h,i}^{\Delta t} = \mathbf{V}_{h,i} \times W_i^{\Delta t}, \quad W_{h,i}^{\Delta t} = W_{h,i} \times W_i^{\Delta t} $$

results in the space-time mixed finite element spaces $\mathbf{V}_{h,i}^{\Delta t} \times W_{h,i}^{\Delta t}$ in each space-time subdomain $\Omega_i^T$. We will also need the spatial variable only spaces

$$ \mathbf{V}_h = \bigoplus \mathbf{V}_{h,i}, \quad W_h = \bigoplus W_{h,i}. $$

Let $\mathcal{T}_{H,ij}$ be a shape-regular finite element partition of $\Gamma_{ij}$, where $H = \max_{i,j} \max_{e \in \mathcal{T}_{H,ij}} \text{diam } e$, see Figure 1, middle. The use of index $H$ indicates a possibly coarser interface grid compared to the subdomain grids, resulting in a multiscale approximation. Let $\mathcal{T}_{ij}^{\Delta T} : 0 = t_{ij}^0 < t_{ij}^1 < \dots < t_{ij}^{N_{ij}} = T$ be a partition of $(0, T)$ corresponding to $\Gamma_{ij}$, which may be different from (and again possibly coarser than) the time-partitions for the neighboring subdomains. Let $\Delta T = \max_{i,j} \max_{1 \le k \le N_{ij}} |t_{ij}^k - t_{ij}^{k-1}|$. Composing $\mathcal{T}_{H,ij}$ and $\mathcal{T}_{ij}^{\Delta T}$ by tensor product gives a space-time partition

$$ \mathcal{T}_{H,ij}^{\Delta T} = \mathcal{T}_{H,ij} \times \mathcal{T}_{ij}^{\Delta T} $$

of the space-time interface $\Gamma_{ij}^T$. Finally, let

$$ \Lambda_{H,ij}^{\Delta T} = \Lambda_{H,ij} \times \Lambda_{ij}^{\Delta T} $$

be a space-time mortar finite element space on $\mathcal{T}_{H,ij}^{\Delta T}$ consisting of continuous or discontinuous piecewise polynomials in space and in time. We will also need the spatial variable only space

$$ \Lambda_H = \bigoplus \Lambda_{H,ij}. $$

Finally, the global space-time finite element spaces are defined as

$$ \mathbf{V}_h^{\Delta t} = \bigoplus \mathbf{V}_{h,i}^{\Delta t}, \quad W_h^{\Delta t} = \bigoplus W_{h,i}^{\Delta t}, \quad \Lambda_H^{\Delta T} = \bigoplus \Lambda_{H,ij}^{\Delta T}. \qquad (3.1) $$

In particular, the Lagrange multiplier will be sought for in the mortar space $\Lambda_H^{\Delta T}$. For the purpose of the analysis, we also define the space of velocities with space-time weakly continuous normal components

$$ \mathbf{V}_{h,0}^{\Delta t} = \{ \mathbf{v} \in \mathbf{V}_h^{\Delta t} : b_{\Gamma}^{T}(\mathbf{v}, \mu) = 0 \quad \forall \mu \in \Lambda_H^{\Delta T} \}. \qquad (3.2) $$

The discrete velocity and pressure spaces inherit the norms $\| \cdot \|_{\mathbf{V}}$ and $\| \cdot \|_{W}$, respectively. The mortar space is equipped with the spatial norm $\|\mu\|_{\Lambda_H} = \|\mu\|_{L^2(\Gamma)}$.
## 3.2 Space-time multiscale mortar mixed finite element method

For the DG time discretization, we introduce the notation for $\varphi(x, \cdot)$, $\phi(x, \cdot) \in W_i^{\Delta t}$, $x \in \Omega_i$, see [50],

$$ \int_0^T \tilde{\partial}_t \varphi \, d\varphi = \sum_{k=1}^{N_i} \int_{t_i^{k-1}}^{t_i^k} \partial_t \varphi \, d\varphi + \sum_{k=1}^{N_i} [\varphi]_{k-1} \phi_{k-1}^+ \quad (3.3) $$

where $[\varphi]_k = \varphi_k^+ - \varphi_k^-$, with $\varphi_k^+ = \lim_{t \to t_i^{k,+}} \varphi$ and $\varphi_k^- = \lim_{t \to t_i^{k,-}} \varphi$.

**Remark 3.1 (Initial value).** In what follows, we will tacitly assume that a function $\varphi(x, \cdot) \in W_i^{\Delta t}$ has an associated initial value $\varphi_0^-$, which will be defined if it is explicitly used.

The space-time multiscale mortar mixed finite element method for approximating (2.5) is: find $\mathbf{u}_h^{\Delta t} \in \mathbf{V}_h^{\Delta t}$, $p_h^{\Delta t} \in W_h^{\Delta t}$, and $\lambda_H^{\Delta T} \in \Lambda_H^{\Delta T}$ such that

$$ a^T(\mathbf{u}_h^{\Delta t}, \mathbf{v}) + b^T(\mathbf{v}, p_h^{\Delta t}) + b_{\Gamma}^T(\mathbf{v}, \lambda_H^{\Delta T}) = 0 \quad \forall \mathbf{v} \in \mathbf{V}_h^{\Delta t}, \quad (3.4a) $$

$$ (\tilde{\partial}_t p_h^{\Delta t}, w)_{\Omega^T} - b^T(\mathbf{u}_h^{\Delta t}, w) = (q, w)_{\Omega^T} \quad \forall w \in W_h^{\Delta t}, \quad (3.4b) $$

$$ b_{\Gamma}^{T}(\mathbf{u}_{h}^{\Delta t}, \mu) = 0 \quad \forall \mu \in \Lambda_{H}^{\Delta T}, \quad (3.4c) $$

where the obvious notation $(\tilde{\partial}_t p_h^{\Delta t}, w)_{\Omega^T} = \sum_i (\tilde{\partial}_t p_h^{\Delta t}, w)_{\Omega_i^T}$ has been used. We note that $(\tilde{\partial}_t p_h^{\Delta t}, w)_{\Omega^T}$ involves the term $((p_h^{\Delta t})_0^+ - (p_h^{\Delta t})_0^-, w_0^+)_\Omega$, see the last term in (3.3) for $k=1$. Here, $(p_h^{\Delta t})_0^+$ is computed by the method, while $(p_h^{\Delta t})_0^-$ is determined by the initial condition. We discuss the construction of initial data in Section 4.4.

The above method provides a highly general and flexible framework, allowing for different spatial and temporal discretizations in different subdomains. We note that according to (3.4c), continuity of the flux is imposed weakly on the space-time interfaces $\Gamma_{ij}^T$, requiring that the jump in flux is orthogonal to the space-time mortar space $\Lambda_{H,ij}^{\Delta T}$. This formulation results in a correct notion of mass conservation across interfaces for time-dependent domain decomposition problems with non-matching grids in both space and time. In the case of discontinuous mortars, (3.4c) implies that the total flux across any space-time interface cell $e \times (t_{ij}^{k-1}, t_{ij}^k)$, $e \in \mathcal{T}_{H,ij}$, is continuous.

# 4 Well-posedness analysis

In this section we analyze the existence, uniqueness, and stability of the solution to (3.4).

## 4.1 Space-time interpolants

We will make use of several space-time interpolants. Let $\mathcal{P}_{h,i}$ be the $L^2$-orthogonal projection onto $W_{h,i}$ and let $\mathcal{P}_i^{\Delta t}$ be the $L^2$-orthogonal projection onto $W_i^{\Delta t}$. We then define the $L^2$-orthogonal projection in space and time on subdomain $\Omega_i$ by

$$ \mathcal{P}_{h,i}^{\Delta t} = \mathcal{P}_{h,i} \times \mathcal{P}_i^{\Delta t} : L^2(0, T; L^2(\Omega_i)) \to W_{h,i}^{\Delta t} $$

and globally by

$$ \mathcal{P}_h^{\Delta t} : L^2(0, T; L^2(\Omega)) \to W_h^{\Delta t}, \quad \mathcal{P}_h^{\Delta t}|_{\Omega_i} = \mathcal{P}_{h,i}^{\Delta t}. $$

Setting $\mathcal{P}_h|_{\Omega_i} = \mathcal{P}_{h,i}$ and $\mathcal{P}_i^{\Delta t}|_{\Omega_i} = \mathcal{P}_i^{\Delta t}$, we will also write $\mathcal{P}_h^{\Delta t} = \mathcal{P}_h \times \mathcal{P}^{\Delta t}$. Since $\nabla \cdot \mathbf{V}_{h,i} = W_{h,i}$, we have, for all $\varphi \in L^2(0, T; L^2(\Omega_i))$,

$$ (\mathcal{P}_h^{\Delta t} \varphi - \varphi, \nabla \cdot \mathbf{v})_{\Omega_i^T} = 0 \quad \forall \mathbf{v} \in \mathbf{V}_{h,i}^{\Delta t}. \quad (4.1) $$
For $\epsilon > 0$, denote $\mathbf{H}^\epsilon(\text{div}; \Omega_i) := \mathbf{H}^\epsilon(\Omega_i) \cap \mathbf{H}(\text{div}; \Omega_i)$. Let $\Pi_{h,i} : \mathbf{H}^\epsilon(\text{div}; \Omega_i) \to \mathbf{V}_{h,i}$ be the canonical mixed interpolant [11] and let

$$ \Pi_{h,i}^{\Delta t} = \Pi_{h,i} \times \mathcal{P}_{i}^{\Delta t} : L^2(0, T; \mathbf{H}^{\epsilon}(\text{div}; \Omega_i)) \to \mathbf{V}_{h,i}^{\Delta t}. $$

In particular, this space-time interpolant satisfies, for all $\psi \in L^2(0, T; \mathbf{H}^\epsilon(\text{div}; \Omega_i))$,

$$ (\nabla \cdot (\Pi_{h,i}^{\Delta t} \psi - \psi), w)_{\Omega_i^T} = 0 \quad \forall w \in W_{h,i}^{\Delta t}, \qquad (4.2a) $$

$$ \langle (\Pi_{h,i}^{\Delta t} \psi - \psi) \cdot \mathbf{n}_i, \mathbf{v} \cdot \mathbf{n}_i \rangle_{\partial \Omega_i^T} = 0 \quad \forall \mathbf{v} \in \mathbf{V}_{h,i}^{\Delta t}, \qquad (4.2b) $$

$$ \| \Pi_{h,i}^{\Delta t} \psi \|_{L^2(0,T; \mathbf{V}_i)} \le C (\| \psi \|_{L^2(0,T; \mathbf{H}^\epsilon(\Omega_i))} + \| \nabla \cdot \psi \|_{L^2(0,T; L^2(\Omega_i))}). \qquad (4.2c) $$

Let $Q_{h,i} : L^2(\partial \Omega_i) \to \mathbf{V}_{h,i} \cdot \mathbf{n}_i$ be the $L^2$-orthogonal projection and let

$$ Q_{h,i}^{\Delta t} = Q_{h,i} \times \mathcal{P}_{i}^{\Delta t} : L^2(0, T; L^2(\partial \Omega_i)) \to \mathbf{V}_{h,i}^{\Delta t} \cdot \mathbf{n}_i. \qquad (4.3) $$

Finally, let $\mathcal{P}_{H,\Gamma_{ij}} : L^2(\Gamma_{ij}) \to \Lambda_{H,ij}$ and $\mathcal{P}_{ij}^{\Delta T} : L^2(0,T) \to \Lambda_{ij}^{\Delta T}$ be the $L^2$-orthogonal projections and let

$$ \mathcal{P}_{H,\Gamma_{ij}}^{\Delta T} = \mathcal{P}_{H,\Gamma_{ij}} \times \mathcal{P}_{ij}^{\Delta T} : L^2(0,T; L^2(\Gamma_{ij})) \to \Lambda_{H,ij}^{\Delta T}, \quad \mathcal{P}_{H,\Gamma|\Gamma_{ij}}^{\Delta T} = \mathcal{P}_{H,\Gamma_{ij}}^{\Delta T} \qquad (4.4) $$

be the mortar space-time $L^2$-orthogonal projection.

## 4.2 Assumptions on the mortar grids

We make the following assumptions on the mortar grids, which are needed to guarantee that the method (3.4) is well posed: there exists a positive constant $C$ independent of the spatial mesh sizes $h$ and $H$ (as well as of the temporal mesh sizes $\Delta t$ and $\Delta T$) such that

$$ \forall \mu \in \Lambda_H, \forall i, j, \quad \|\mu\|_{\Gamma_{ij}} \le C(\|Q_{h,i}\mu\|_{\Gamma_{ij}} + \|Q_{h,j}\mu\|_{\Gamma_{ij}}), \qquad (4.5a) $$

$$ \forall i, j, \quad \Lambda_{ij}^{\Delta T} \subset W_i^{\Delta t} \cap W_j^{\Delta t}. \qquad (4.5b) $$

The spatial mortar assumption (4.5a) is the same as the assumption made in [3,5]. Note that it is in particular satisfied with $C = \frac{1}{2}$ when $\mathcal{T}_{H,ij}$ is a coarsening of both $\mathcal{T}_{h,i}$ and $\mathcal{T}_{h,j}$ on the interface $\Gamma_{ij}$ and the space $\Lambda_{H,ij}$ consists of discontinuous piecewise polynomials contained in $\mathbf{V}_{h,i} \cdot \mathbf{n}_i$ and $\mathbf{V}_{h,j} \cdot \mathbf{n}_j$ on $\Gamma_{ij}$. In general, it requires that the mortar space $\Lambda_H$ is sufficiently coarse, so that it is controlled by the normal traces of the neighboring subdomain velocity spaces.

The temporal mortar assumption (4.5b) similarly provides control of the mortar time discretization by the subdomain time discretizations. It requires that each subdomain time discretization be a refinement of the mortar time discretization. We also note that (4.5a) and (4.5b) imply

$$ \forall \mu \in \Lambda_H^{\Delta T}, \forall i, j, \quad \|\mu\|_{L^2(0,T;L^2(\Gamma_{ij}))} \le C(\|Q_{h,i}^{\Delta t}\mu\|_{L^2(0,T;L^2(\Gamma_{ij}))} + \|Q_{h,j}^{\Delta t}\mu\|_{L^2(0,T;L^2(\Gamma_{ij}))}) \qquad (4.6) $$

for a constant $C$ independent of $h, H, \Delta t$, and $\Delta T$.

## 4.3 Discrete inf-sup conditions

Recall the form $b^T(\cdot, \cdot)$ from (2.4b). Under the above assumptions on the mortar grids, the weakly continuous velocity space $\mathbf{V}_{h,0}^{\Delta t}$ of (3.2) satisfies the following inf-sup condition.

**Lemma 4.1** (Discrete divergence inf-sup condition on $\mathbf{V}_{h,0}^{\Delta t}$). Let (4.5) hold. Then there exists a constant $\beta > 0$, independent of $h, H, \Delta t$, and $\Delta T$, such that

$$ \forall w \in W_h^{\Delta t}, \quad \sup_{0 \neq \mathbf{v} \in \mathbf{V}_{h,0}^{\Delta t}} \frac{b^T(\mathbf{v}, w)}{\|\mathbf{v}\|_{L^2(0,T;\mathbf{V})}} \ge \beta \|w\|_{L^2(0,T;L^2(\Omega))}. \qquad (4.7) $$
Proof. Let $V_{h,0} = \{\mathbf{v} \in V_h : b_{\Gamma}(\mathbf{v}, \mu) = 0 \quad \forall \mu \in \Lambda_H\}$. It is shown in [3, 5] that if (4.5a) holds, then there is an interpolant $\Pi_{h,0}: \mathbf{H}^{\frac{1}{2}+\epsilon}(\text{div}; \Omega) \to V_{h,0}$ such that, for all $\psi \in \mathbf{H}^{\frac{1}{2}+\epsilon}(\text{div}; \Omega)$,

$$ \sum_i (\nabla \cdot (\Pi_{h,0}\psi - \psi), w)_{\Omega_i} = 0 \quad \forall w \in W_h, \quad (4.8a) $$

$$ \|\Pi_{h,0}\psi\|_{V} \le C(\|\psi\|_{\mathbf{H}^{\frac{1}{2}+\epsilon}(\Omega)} + \|\nabla \cdot \psi\|_{L^2(\Omega)}), \quad (4.8b) $$

for a constant C independent of h and H. Define

$$ \Pi_{h,0}^{\Delta t} = \Pi_{h,0} \times \mathcal{P}^{\Delta t}. $$

We claim that $\Pi_{h,0}^{\Delta t}: L^2(0, T; \mathbf{H}^{\frac{1}{2}+\epsilon}(\text{div}; \Omega)) \to V_{h,0}^{\Delta t}$. To see this, note first that, for all functions $\psi \in L^2(0, T; \mathbf{H}^{\frac{1}{2}+\epsilon}(\text{div}; \Omega))$, clearly $\Pi_{h,0}^{\Delta t}\psi \in V_h^{\Delta t}$. Thus (4.5b) implies

$$ \begin{aligned} b_{\Gamma}^{T}(\Pi_{h,0}^{\Delta t}\psi, \mu) &= \sum_{i} \int_{0}^{T} \langle \Pi_{h,0}^{\Delta t}\psi \cdot \mathbf{n}_{i}, \mu \rangle_{\Gamma_{i}} = \sum_{i} \int_{0}^{T} \langle \Pi_{h,0}\psi \cdot \mathbf{n}_{i}, \mu \rangle_{\Gamma_{i}} \\ &= \int_{0}^{T} b_{\Gamma}(\Pi_{h,0}\psi, \mu) = 0 \quad \forall \mu \in \Lambda_{H}^{\Delta T}, \end{aligned} $$

i.e., indeed $\Pi_{h,0}^{\Delta t}\psi \in V_{h,0}^{\Delta t}$ by virtue of (3.2). Moreover, (4.8a) and (4.8b) imply

$$ \sum_i (\nabla \cdot (\Pi_{h,0}^{\Delta t}\psi - \psi), w)_{\Omega_i^T} = 0 \quad \forall w \in W_h^{\Delta t}, \quad (4.9a) $$

$$ \|\Pi_{h,0}^{\Delta t}\psi\|_{L^2(0,T;V)} \le C(\|\psi\|_{L^2(0,T;\mathbf{H}^{\frac{1}{2}+\epsilon}(\Omega))} + \|\nabla \cdot \psi\|_{L^2(0,T;L^2(\Omega))}). \quad (4.9b) $$

The inf-sup condition (4.7) then follows from the classical continuous inf-sup condition for $b^T(\cdot, \cdot)$, the existence of the interpolant $\Pi_{h,0}^{\Delta t}$, and Fortin's lemma [11]. $\square$

To control the mortar variable, we need the following mortar inf-sup condition.

**Lemma 4.2** (Discrete mortar inf-sup condition on $V_h^{\Delta t}$). Let (4.6) hold. Then there exists a constant $\beta_{\Gamma} > 0$, independent of $h$, $H$, $\Delta t$, and $\Delta T$, such that

$$ \forall \mu \in \Lambda_{H}^{\Delta T}, \quad \sup_{0 \neq \mathbf{v} \in V_{h}^{\Delta t}} \frac{b_{\Gamma}^{T}(\mathbf{v}, \mu)}{\|\mathbf{v}\|_{L^2(0,T;\mathbf{v})}} \ge \beta_{\Gamma} \|\mu\|_{L^2(0,T;L^2(\Gamma))}. \quad (4.10) $$

Proof. Let $\mu \in \Lambda_{H}^{\Delta T}$ be given. In the following we assume that $\mu$ is extended by zero on $\partial\Omega$. We consider a set of auxiliary subdomain problems. Let $\varphi_i(x, t)$ be the solution for a.e. $t \in (0, T)$ of the problem

$$ \nabla \cdot \nabla \varphi_i(\cdot, t) = \overline{(Q_{h,i}^{\Delta t} \mu)(\cdot, t)} \quad \text{in } \Omega_i, \quad (4.11a) $$

$$ \nabla \varphi_i(\cdot, t) \cdot \mathbf{n}_i = (Q_{h,i}^{\Delta t} \mu)(\cdot, t) \quad \text{on } \partial\Omega_i, \quad (4.11b) $$

where $\overline{Q_{h,i}^{\Delta t}} \mu$ denotes the mean value of $Q_{h,i}^{\Delta t} \mu$ on $\partial\Omega_i$. Let $\psi_i = \nabla \varphi_i$. Elliptic regularity [30, 43] implies that for a.e. $t \in (0, T)$,

$$ \|\psi_i\|_{\frac{1}{2}, \Omega_i} + \|\nabla \cdot \psi_i\|_{\Omega_i} \le C \|Q_{h,i}^{\Delta t} \mu\|_{\partial \Omega_i}. \quad (4.12) $$
Let $\mathbf{v}_i = \mathbf{\Pi}_{h,i}^{\Delta t} \psi_i \in \mathbf{V}_{h,i}^{\Delta t}$. Note that (4.2b) together with (4.3) and (4.11b) imply that $\mathbf{v}_i \cdot \mathbf{n}_i = \mathcal{Q}_{h,i}^{\Delta t} \mu$ on $\partial\Omega_i$. Thus, using definition (2.4c) of $b_{\Gamma}^T$, the fact that $\mu$ is extended by zero on $\partial\Omega_i \setminus \Gamma_i$, and definition (4.3) of the projection $\mathcal{Q}_{h,i}^{\Delta t}$, we have

$$ 
\begin{aligned}
b_{\Gamma}^{T}(\mathbf{v}, \mu) &= \sum_{i} \langle \mathbf{\Pi}_{h,i}^{\Delta t} \psi_i \cdot \mathbf{n}_i, \mu \rangle_{\Gamma_i^T} = \sum_{i} \langle \mathbf{\Pi}_{h,i}^{\Delta t} \psi_i \cdot \mathbf{n}_i, \mu \rangle_{\partial\Omega_i^T} = \sum_{i} \langle \mathbf{\Pi}_{h,i}^{\Delta t} \psi_i \cdot \mathbf{n}_i, \mathcal{Q}_{h,i}^{\Delta t} \mu \rangle_{\partial\Omega_i^T} \\
&= \sum_{i} \| \mathcal{Q}_{h,i}^{\Delta t} \mu \|_{L^2(0,T;L^2(\partial\Omega_i))}^2 \geq C \sum_{i} \|\mu\|_{L^2(0,T;L^2(\Gamma_i))}^2,
\end{aligned} 
\quad (4.13) $$

where we used (4.6) in the inequality. On the other hand, (4.2c) with $\epsilon = \frac{1}{2}$ and (4.12), along with the stability of $L^2$-orthogonal projection $\mathcal{Q}_{h,i}^{\Delta t}$, imply

$$ \|\mathbf{v}_i\|_{L^2(0,T;\mathbf{V}_i)} \leq C \|\mu\|_{L^2(0,T;L^2(\Gamma_i))}. \quad (4.14) $$

The assertion of the lemma follows from combining (4.13) and (4.14). $\square$

## 4.4 Initial data

We next discuss the construction of discrete initial data for all variables. The data need to be compatible in the sense that they satisfy the equations without time derivatives in the method, (3.4a) and (3.4c). Recall that we are given initial pressure datum $p(0) = p_0 \in H_0^1(\Omega)$ with $\nabla \cdot K\nabla p_0 \in L^2(\Omega)$. Let us define $\mathbf{u}_0 = -K\nabla p_0$ and $\lambda_0 = p_0|_{\Gamma}$. Then the solution to (2.5) satisfies $\mathbf{u}(0) = \mathbf{u}_0$ and $\lambda(0) = \lambda_0$. Moreover, we have

$$ 
\begin{aligned}
& a(\mathbf{u}_0, \mathbf{v}) + b(\mathbf{v}, p_0) + b_{\Gamma}(\mathbf{v}, \lambda_0) = 0 \quad \forall \mathbf{v} \in \mathbf{V}_h, \\
& b_{\Gamma}(\mathbf{u}_0, \mu) = 0 \quad \forall \mu \in \Lambda_H.
\end{aligned} 
$$

Next, define the discrete initial data $(\mathbf{u}_{h,0}, p_{h,0}, \lambda_{H,0}) \in \mathbf{V}_h \times W_h \times \Lambda_H$ as the elliptic projection of $(\mathbf{u}_0, p_0, \lambda_0)$, i.e., the unique solution to the problem

$$ a(\mathbf{u}_{h,0}, \mathbf{v}) + b(\mathbf{v}, p_{h,0}) + b_{\Gamma}(\mathbf{v}, \lambda_{H,0}) = a(\mathbf{u}_0, \mathbf{v}) + b(\mathbf{v}, p_0) + b_{\Gamma}(\mathbf{v}, \lambda_0) = 0 \quad \forall \mathbf{v} \in \mathbf{V}_h, \quad (4.15a) $$

$$ b(\mathbf{u}_{h,0}, w) = b(\mathbf{u}_0, w) = -(\nabla \cdot K \nabla p_0, w) \quad \forall w \in W_h, \quad (4.15b) $$

$$ b_{\Gamma}(\mathbf{u}_{h,0}, \mu) = b_{\Gamma}(\mathbf{u}_0, \mu) = 0 \quad \forall \mu \in \Lambda_H. \quad (4.15c) $$

The well-posedness of (4.15) is shown in [3, 5] under the spatial mortar assumption (4.5a). In particular, it follows from the analysis in [3, 5] that

$$ \|\mathbf{u}_{h,0}\|_{\mathbf{V}} + \|p_{h,0}\|_{W} + \|\lambda_{H,0}\|_{\Lambda_H} \leq C \|\nabla \cdot K \nabla p_0\|_{\Omega}, \quad (4.16) $$

$$ 
\begin{aligned}
& \|\mathbf{u}_0 - \mathbf{u}_{h,0}\|_{\mathbf{V}} + \|p_0 - p_{h,0}\|_{W} + \|\lambda_0 - \lambda_{H,0}\|_{\Lambda_H} \\
& \leq C (\|\mathbf{u}_0 - \mathbf{\Pi}_{h,0}\mathbf{u}_0\|_{\mathbf{V}} + \|p_0 - P_h p_0\|_{W} + \|\lambda_0 - P_{H,\Gamma}\lambda_0\|_{\Lambda_H}).
\end{aligned} 
\quad (4.17) $$

We now set

$$ (p_h^{\Delta t})_0^- = p_{h,0}, \quad (\mathbf{u}_h^{\Delta t})_0^- = \mathbf{u}_{h,0}, \quad (\lambda_H^{\Delta T})_0^- = \lambda_{h,0}. \quad (4.18) $$

As we noted earlier, $(p_h^{\Delta t})_0^-$ provides initial condition for the method (3.4), cf. (3.4b). The data $(\mathbf{u}_h^{\Delta t})_0^-$ and $(\lambda_H^{\Delta T})_0^-$ are not needed in the method, but it will be utilized in the analysis of $\nabla \cdot \mathbf{u}_h^{\Delta t}$ in Section 6.
## 4.5 Existence, uniqueness, and stability with respect to data

In the analysis we will utilize the following auxiliary result.

**Lemma 4.3 (Summation in time).** For all $\Omega_i$ and for any $\phi(x, \cdot) \in W_i^{\Delta t}$, $x \in \Omega_i$, there holds

$$ \int_0^T \tilde{\partial}_t \varphi \varphi = \frac{1}{2} ((\varphi_{N_i}^-)^2 - (\varphi_0^-)^2) + \frac{1}{2} \sum_{k=1}^{N_i} ([\varphi]_{k-1})^2. \quad (4.19) $$

*Proof.* Using the definition (3.3) of $\tilde{\partial}_t w$, we have

$$ \begin{aligned} \int_0^T \tilde{\partial}_t \varphi \varphi &= \sum_{k=1}^{N_i} \int_{t_{k-1}}^{t_k} \frac{1}{2} \frac{\partial}{\partial t} \varphi^2 + \sum_{k=1}^{N_i} [\varphi]_{k-1} \varphi_{k-1}^+ \\ &= \frac{1}{2} \sum_{k=1}^{N_i} ((\varphi_k^-)^2 - (\varphi_{k-1}^+)^2 + (\varphi_{k-1}^+)^2 - (\varphi_{k-1}^-)^2 + (\varphi_{k-1}^- - \varphi_{k-1}^+)^2) \\ &= \frac{1}{2} ((\varphi_{N_i}^-)^2 - (\varphi_0^-)^2) + \frac{1}{2} \sum_{k=1}^{N_i} (\varphi_{k-1}^+ - \varphi_{k-1}^-)^2. \end{aligned} $$

□

To simplify the presentation, we introduce the notation

$$ \|\varphi\|_{\text{DG}}^2 = \sum_i (\|\varphi_{N_i}^-\|_{\Omega_i}^2 + \sum_{k=1}^{N_i} \|[\varphi]_{k-1}\|_{\Omega_i}^2). \quad (4.20) $$

**Theorem 4.4 (Existence and uniqueness of the discrete solution, stability with respect to data).** Assume that conditions (4.5) hold. Then the space-time mortar method (3.4) has a unique solution. Moreover, for some constant $C > 0$ independent of $h$, $H$, $\Delta t$, and $\Delta T$,

$$ \|p_h^{\Delta t}\|_{\text{DG}} + \|u_h^{\Delta t}\|_{\Omega^T} + \|p_h^{\Delta t}\|_{\Omega^T} + \|\lambda_H^{\Delta T}\|_{\Gamma^T} \le C(\|q\|_{\Omega^T} + \|\nabla \cdot K \nabla p_0\|_{\Omega}). \quad (4.21) $$

*Proof.* We begin with establishing the stability bound (4.21). Taking $v = u_h^{\Delta t}$, $w = p_h^{\Delta t}$, and $\mu = \lambda_H^{\Delta T}$ in (3.4) and combining the equations, we obtain, using (4.19) and Young's inequality,

$$ \frac{1}{2} \sum_i \left( \| (p_h^{\Delta t})_{N_i}^- \|_{\Omega_i}^2 + \sum_{k=1}^{N_i} \| [p_h^{\Delta t}]_{k-1} \|_{\Omega_i}^2 \right) + \| K^{-\frac{1}{2}} u_h^{\Delta t} \|_{\Omega^T}^2 \le \frac{\epsilon}{2} \| p_h^{\Delta t} \|_{\Omega^T}^2 + \frac{1}{2\epsilon} \| q \|_{\Omega^T}^2 + \frac{1}{2} \| p_{h,0} \|_{\Omega}^2. $$

The inf-sup condition for the weakly continuous velocity (4.7) and (3.4a) imply

$$ \|p_h^{\Delta t}\|_{\Omega^T} \le C \|K^{-\frac{1}{2}} u_h^{\Delta t}\|_{\Omega^T}. $$

Furthermore, the mortar inf-sup condition (4.10) and (3.4a) imply

$$ \|\lambda_H^{\Delta T}\|_{\Gamma^T} \le C(\|K^{-\frac{1}{2}} u_h^{\Delta T}\|_{\Omega^T} + \|p_h^{\Delta T}\|_{\Omega^T}). $$

Combining the above three inequalities, taking $\epsilon$ sufficiently small, and using (2.2) and (4.16), we obtain (4.21). The existence and uniqueness of a solution follows from (4.21) by taking $q = 0$ and $p_0 = 0$. □

**Remark 4.5 (Control of divergence).** Control on $\|\nabla \cdot u_h^{\Delta t}\|_{L^2(0,T;L^2(\Omega_i))}$ can be established under the assumption of matching time steps between subdomains and choosing the mortar finite element space in time to match the subdomains. We present this result later in Section 6, along with improved error estimates.
# 5 A priori error analysis

In this section we derive a priori error estimates for the solution of the space-time mortar MFE method (3.4).

## 5.1 Approximation properties of the space-time interpolants

Assume that the spaces $V_h^{\Delta t}$ and $W_h^{\Delta t}$ from (3.1) contain on each space-time element polynomials in $P_k$ and $P_l$, respectively, in space and $P_q$ in time, where $P_r$ denotes the space of polynomials of degree up to $r$. Let $\Lambda_H^{\Delta T}$ contain on each space-time mortar element polynomials in $P_m$ in space and $P_s$ in time. We have the following approximation properties for the space-time interpolants $\mathcal{P}_h^{\Delta t}$ and $\mathcal{P}_{H,\Gamma}^{\Delta T}$ of Section 4.1 and $\Pi_{h,0}^{\Delta t}$ of the proof of Lemma 4.1:

$$ \|\psi - \Pi_{h,0}^{\Delta t} \psi\|_{\Omega^T} \le C \sum_i \|\psi\|_{H^{r_q}(0,T;H^{r_k}(\Omega_i))} (h^{r_k} + \Delta t^{r_q}) + C \|\psi\|_{H^{r_q}(0,T;\tilde{H}^{r_k+\frac{1}{2}}(\Omega))} (h^{\tilde{r}_k} H^{\frac{1}{2}} + \Delta t^{r_q}), \\ 0 < r_k \le k+1, \quad 0 < \tilde{r}_k \le k+1, \quad 0 \le r_q \le q+1, \qquad (5.1a) $$

$$ \|\varphi - \mathcal{P}_h^{\Delta t} \varphi\|_{\Omega_i^T} \le C \|\varphi\|_{H^{r_q}(0,T;H^{r_l}(\Omega_i))} (h^{r_l} + \Delta t^{r_q}), \quad 0 \le r_l \le l+1, \quad 0 \le r_q \le q+1, \qquad (5.1b) $$

$$ \|\varphi - \mathcal{P}_{H,\Gamma}^{\Delta T} \varphi\|_{\Gamma_{ij}^T} \le C \|\varphi\|_{H^{r_s}(0,T;H^{r_m}(\Gamma_{ij}))} (H^{r_m} + \Delta T^{r_s}), \quad 0 \le r_m \le m+1, \quad 0 \le r_s \le s+1. \qquad (5.1c) $$

Bound (5.1a) follows from the approximation properties of $\Pi_{h,0}$ obtained in [3, 5]. Bounds (5.1b) and (5.1c) are standard approximation properties of the $L^2$ projection [13].

In the analysis we will also use the following approximation property, which follows from the stability of the $L^2$ projection in $L^\infty$ [14]:

$$ \|\varphi - \mathcal{P}_h^{\Delta t} \varphi\|_{L^\infty(0,T;L^2(\Omega_i))} \le C \|\varphi\|_{W^{r_q,\infty}(0,T;H^{r_l}(\Omega_i))} (h^{r_l} + \Delta t^{r_q}), \quad 0 \le r_l \le l+1, \quad 0 \le r_q \le q+1. \qquad (5.2) $$

We also recall the well-known discrete trace (inverse) inequality for a quasi-uniform mesh $\mathcal{T}_{h,i}$: for all $\mathbf{v} \in \mathbf{V}_{h,i}$, $\|\mathbf{v} \cdot \mathbf{n}_i\|_{\Gamma_i} \le C h_i^{-\frac{1}{2}} \|\mathbf{v}\|_{\Omega_i}$. This implies

$$ \forall \mathbf{v} \in \mathbf{V}_{h,i}^{\Delta t}, \quad \|\mathbf{v} \cdot \mathbf{n}_i\|_{\Gamma_i^T} \le C h_i^{-\frac{1}{2}} \|\mathbf{v}\|_{\Omega_i^T}. \qquad (5.3) $$

## 5.2 A priori error estimate

We proceed with the error estimate for the space-time mortar MFE method (3.4).

**Theorem 5.1 (A priori error estimate).** Assume that conditions (4.5) hold and that the solution to (2.5) is sufficiently smooth. Let the space and time meshes $\mathcal{T}_{h,i}$ and $\mathcal{T}_i^{\Delta t}$ be quasi-uniform, as well as $h \le Ch_i$ and $\Delta t \le C\Delta t_i \forall i$. Then there exists a constant $C > 0$ independent of the mesh sizes $h$, $H$, $\Delta t$, and $\Delta T$, such that the solution to the space-time mortar MFE method (3.4) satisfies

$$ \|p - p_h^{\Delta t}\|_{\mathrm{DG}} + \|\mathbf{u} - \mathbf{u}_h^{\Delta t}\|_{\Omega^T} + \|p - p_h^{\Delta t}\|_{\Omega^T} + \|\lambda - \lambda_H^{\Delta T}\|_{\Gamma^T} \\ \le C \left( \sum_i \|\mathbf{u}\|_{H^{r_q}(0,T;H^{r_k}(\Omega_i))} (h^{r_k} + \Delta t^{r_q}) + \|\mathbf{u}\|_{H^{r_q}(0,T;\tilde{H}^{r_k+\frac{1}{2}}(\Omega))} (h^{\tilde{r}_k} H^{\frac{1}{2}} + \Delta t^{r_q}) \right. \\ \left. + \sum_i \|p\|_{W^{r_q,\infty}(0,T;H^{r_l}(\Omega_i))} \Delta t^{-\frac{1}{2}} (h^{r_l} + \Delta t^{r_q}) + \sum_{i,j} \|\lambda\|_{H^{r_s}(0,T;H^{r_m}(\Gamma_{ij}))} h^{-\frac{1}{2}} (H^{r_m} + \Delta T^{r_s}) \right), \\ 0 < r_k \le k+1, \quad 0 < \tilde{r}_k \le k+1, \quad 0 \le r_q \le q+1, \quad 0 \le r_l \le l+1, \quad 0 \le r_m \le m+1, \quad 0 \le r_s \le s+1. \qquad (5.4) $$
Proof. For the purpose of the analysis, we consider the following equivalent formulation of (3.4) in the space of weakly continuous velocities $V_{h,0}^{\Delta t}$ given by (3.2): find $u_{h,0}^{\Delta t} \in V_{h,0}^{\Delta t}$ and $p_h^{\Delta t} \in W_h^{\Delta t}$ such that $(p_h^{\Delta t})^{-} = p_{h,0}$ and

$$a^T(\mathbf{u}_h^{\Delta t}, \mathbf{v}) + b^T(\mathbf{v}, p_h^{\Delta t}) = 0 \quad \forall \mathbf{v} \in V_{h,0}^{\Delta t}, \qquad (5.5a)$$

$$(\tilde{\partial}_t p_h^{\Delta t}, w)_{\Omega^T} - b^T(\mathbf{u}_h^{\Delta t}, w) = (q, w)_{\Omega^T} \quad \forall w \in W_h^{\Delta t}. \qquad (5.5b)$$

The fact that $\mathcal{P}_{H,\Gamma}^{\Delta T}$ defined in (4.4) maps to $\Lambda_H^{\Delta T}$ and definition (3.2) imply that $b_{\Gamma}^T(\mathbf{v}, \mathcal{P}_{H,\Gamma}^{\Delta T}\lambda) = 0$ for all $\mathbf{v} \in V_{h,0}^{\Delta t}$, where $\lambda = p|_{\Gamma}$ is the pressure trace from (2.5). Then, subtracting (5.5a)-(5.5b) from (2.5a)-(2.5b), we obtain the error equations

$$(\partial_t p - \tilde{\partial}_t p_h^{\Delta t}, w)_{\Omega^T} - b^T(\Pi_{h,0}^{\Delta t}\mathbf{u} - \mathbf{u}_h^{\Delta t}, w) = 0 \quad \forall w \in W_h^{\Delta t}, \qquad (5.6a)$$

$$(\partial_t p - \tilde{\partial}_t p_h^{\Delta t}, w)_{\Omega^T} - b^T(\Pi_{h,0}^{\Delta t}\mathbf{u} - \mathbf{u}_h^{\Delta t}, w) = 0 \quad \forall w \in W_h^{\Delta t}, \qquad (5.6b)$$

where we have also used (4.1) and (4.9a) to incorporate the interpolants $\mathcal{P}_h^{\Delta t}$ and $\Pi_{h,0}^{\Delta t}$. We take $\mathbf{v} = \Pi_{h,0}^{\Delta t}\mathbf{u} - \mathbf{u}_h^{\Delta t}$ and $w = \mathcal{P}_h^{\Delta t}p - p_h^{\Delta t}$ and sum the two equations, resulting in

$$\begin{aligned} & a^T(\Pi_{h,0}^{\Delta t}\mathbf{u} - \mathbf{u}_h^{\Delta t}, \Pi_{h,0}^{\Delta t}\mathbf{u} - \mathbf{u}_h^{\Delta t}) + (\partial_t p - \tilde{\partial}_t p_h^{\Delta t}, \mathcal{P}_h^{\Delta t}p - p_h^{\Delta t})_{\Omega^T} \\ & = a^T(\Pi_{h,0}^{\Delta t}\mathbf{u} - \mathbf{u}, \Pi_{h,0}^{\Delta t}\mathbf{u} - \mathbf{u}_h^{\Delta t}) - b_{\Gamma}^T(\Pi_{h,0}^{\Delta t}\mathbf{u} - \mathbf{u}_h^{\Delta t}, \lambda - \mathcal{P}_{H,\Gamma}^{\Delta T}\lambda). \end{aligned} \qquad (5.7)$$

For the second term on the left of (5.7), restricted to a subdomain, we write

$$\begin{aligned} & \int_0^T (\partial_t p - \tilde{\partial}_t p_h^{\Delta t}, \mathcal{P}_h^{\Delta t} p - p_h^{\Delta t})_{\Omega_i} = \int_0^T (\tilde{\partial}_t(p - p_h^{\Delta t}), \mathcal{P}_h^{\Delta t} p - p_h^{\Delta t})_{\Omega_i} \\ & = \int_0^T (\tilde{\partial}_t(p - p_h^{\Delta t}), p - p_h^{\Delta t})_{\Omega_i} + \int_0^T (\tilde{\partial}_t(p - p_h^{\Delta t}), \mathcal{P}_h^{\Delta t} p - p)_{\Omega_i} =: I_1^i + I_2^i. \end{aligned} \qquad (5.8)$$

Using (4.19) and notation (4.20), for the first term, we have

$$\sum_i I_1^i = \frac{1}{2} \|p - p_h^{\Delta t}\|_{\text{DG}}^2 - \frac{1}{2} \|p_0 - p_{h,0}\|_{\Omega}^2. \qquad (5.9)$$

Using (3.3), the second term is

$$I_2^i = \sum_{k=1}^{N_i} \int_{t_i^{k-1}}^{t_i^k} (\partial_t(p - p_h^{\Delta t}), \mathcal{P}_h^{\Delta t}p - p)_{\Omega_i} + \sum_{k=1}^{N_i} ([p - p_h^{\Delta t}]_{k-1}, (\mathcal{P}_h^{\Delta t}p - p)_{k-1}^{+})_{\Omega_i}. \qquad (5.10)$$

For the first term on the right above, using the orthogonality property of $\mathcal{P}_h^{\Delta t}$, we develop

$$\begin{aligned} & \sum_{k=1}^{N_i} \int_{t_i^{k-1}}^{t_i^k} (\partial_t(p - p_h^{\Delta t}), \mathcal{P}_h^{\Delta t}p - p)_{\Omega_i} = \sum_{k=1}^{N_i} \int_{t_i^{k-1}}^{t_i^k} (\partial_t(p - \mathcal{P}_h^{\Delta t}p), \mathcal{P}_h^{\Delta t}p - p)_{\Omega_i} \\ & = -\frac{1}{2} \sum_{k=1}^{N_i} \int_{t_i^{k-1}}^{t_i^k} \partial_t \|\mathcal{P}_h^{\Delta t}p - p\|_{\Omega_i}^2 = -\frac{1}{2} \sum_{k=1}^{N_i} \|\mathcal{P}_h^{\Delta t}p - p\|_{\Omega_i}^2 \Big|_{t_i^{k-1}}^{t_i^k}. \end{aligned} \qquad (5.11)$$
Combining (5.7)–(5.11), and using (2.2) and the Cauchy–Schwarz and Young inequalities, we obtain,

$$
\begin{align*}
& \| \mathbf{\Pi}_{h,0}^{\Delta t} \mathbf{u} - \mathbf{u}_h^{\Delta t} \|_{\Omega^T}^2 + \| p - p_h^{\Delta t} \|_{\text{DG}}^2 \\
& \le C \left( \| \mathbf{\Pi}_{h,0}^{\Delta t} \mathbf{u} - \mathbf{u} \|_{\Omega^T} \| \mathbf{\Pi}_{h,0}^{\Delta t} \mathbf{u} - \mathbf{u}_h^{\Delta t} \|_{\Omega^T} + \sum_i \| (\mathbf{\Pi}_{h,0}^{\Delta t} \mathbf{u} - \mathbf{u}_h^{\Delta t}) \cdot \mathbf{n}_i \|_{\Gamma_i^T} \| \lambda - \mathcal{P}_{H,\Gamma}^{\Delta T} \lambda \|_{\Gamma_i^T} \right. \\
& \qquad \left. + \sum_i \| [p - p_h^{\Delta t}]_{k-1} \|_{\Omega_i} \| (\mathcal{P}_h^{\Delta t} p - p)_{k-1}^+ \|_{\Omega_i} + \sum_i \sum_{k=1}^{N_i} \| \mathcal{P}_h^{\Delta t} p - p \|_{\Omega_i}^2 \Big|_{t^{k-1}}^{t^k} + \| p_0 - p_{h,0} \|_{\Omega}^2 \right) \\
& \le \epsilon \left( \| \mathbf{\Pi}_{h,0}^{\Delta t} \mathbf{u} - \mathbf{u}_h^{\Delta t} \|_{\Omega^T}^2 + \| p - p_h^{\Delta t} \|_{\text{DG}}^2 \right) \\
& \qquad + C_{\epsilon} \left( \| \mathbf{\Pi}_{h,0}^{\Delta t} \mathbf{u} - \mathbf{u} \|_{\Omega^T}^2 + h^{-1} \| \lambda - \mathcal{P}_{H,\Gamma}^{\Delta T} \lambda \|_{\Gamma^T}^2 + \sum_i \sum_{k=1}^{N_i} \| (\mathcal{P}_h^{\Delta t} p - p)_{k-1}^+ \|_{\Omega_i}^2 \right) \\
& \qquad + C \left( \sum_i \sum_{k=1}^{N_i} \| \mathcal{P}_h^{\Delta t} p - p \|_{\Omega_i}^2 \Big|_{t^{k-1}}^{t^k} + \| p_0 - p_{h,0} \|_{\Omega}^2 \right),
\end{align*}
$$

where we used the trace (inverse) inequality (5.3) for a quasi-uniform mesh $\mathcal{T}_{h,i}$ and $h \le Ch_i$ in the last estimate. Taking $\epsilon$ sufficiently small gives

$$
\begin{equation}
\begin{split}
\| \mathbf{\Pi}_{h,0}^{\Delta t} \mathbf{u} - \mathbf{u}_h^{\Delta t} \|_{\Omega^T} + \| p - p_h^{\Delta t} \|_{\text{DG}} &\le C \left( \| \mathbf{\Pi}_{h,0}^{\Delta t} \mathbf{u} - \mathbf{u} \|_{\Omega^T} + h^{-\frac{1}{2}} \| \lambda - \mathcal{P}_{H,\Gamma}^{\Delta T} \lambda \|_{\Gamma^T} \right. \\
&\qquad \left. + \Delta t^{-\frac{1}{2}} \| \mathcal{P}_h^{\Delta t} p - p \|_{L^{\infty}(0,T;L^2(\Omega))} + \| p_0 - p_{h,0} \|_{\Omega} \right),
\end{split}
\tag{5.12}
\end{equation}
$$

where we used that $N_i \le \frac{CT}{\Delta t_i}$ for a quasi-uniform time mesh $\mathcal{T}_i^{\Delta t}$ and $\Delta t \le C\Delta t_i$ to obtain the factor $\Delta t^{-\frac{1}{2}}$.

Next, the inf-sup condition for the weakly continuous velocity (4.7) and (5.6a) imply, using (5.3) and $h \le Ch_i$,

$$
\| \mathcal{P}_h^{\Delta t} p - p_h^{\Delta t} \|_{\Omega^T} \le C \left( \| \mathbf{u} - \mathbf{u}_h^{\Delta t} \|_{\Omega^T} + h^{-\frac{1}{2}} \| \lambda - \mathcal{P}_{H,\Gamma}^{\Delta T} \lambda \|_{\Gamma^T} \right). \quad (5.13)
$$

Finally, to obtain a bound on $\lambda_H^{\Delta T}$, we subtract (3.4a) from (2.5a), to obtain the error equation

$$
a^T(\mathbf{u} - \mathbf{u}_h^{\Delta t}, \mathbf{v}) + b^T(\mathbf{v}, \mathbf{p} - \mathbf{p}_h^{\Delta t}) + b_{\Gamma}^T(\mathbf{v}, \mathcal{P}_{H,\Gamma}^{\Delta T}\lambda - \lambda_H^{\Delta T}) = b_{\Gamma}^T(\mathbf{v}, \mathcal{P}_{H,\Gamma}^{\Delta T}\lambda - \lambda) \quad \forall \mathbf{v} \in \mathbf{V}_h^{\Delta t}. \quad (5.14)
$$

The mortar inf-sup condition (4.10) and (5.14) imply, using (5.3) and $h \le Ch_i$,

$$
\| \mathcal{P}_{H,\Gamma}^{\Delta T} \lambda - \lambda_H^{\Delta T} \|_{\Gamma^T} \le C \left( \| \mathbf{u} - \mathbf{u}_h^{\Delta T} \|_{\Omega^T} + \| p - p_h^{\Delta T} \|_{\Omega^T} + h^{-\frac{1}{2}} \| \lambda - \mathcal{P}_{H,\Gamma}^{\Delta T} \lambda \|_{\Gamma^T} \right). \quad (5.15)
$$

The assertion of the theorem follows from combining (5.12), (5.13), and (5.15) and using the triangle inequality, (4.17), and the approximation bounds (5.1)–(5.2). $\square$

**Remark 5.2** (The factors $\Delta t^{-\frac{1}{2}}$ and $h^{-\frac{1}{2}}$ and appropriate choice of the polynomial degrees $m$ and $s$). The term $h^{-\frac{1}{2}}(H^{rm} + \Delta T^{rs})$ in the error bound appears due to the use of the discrete trace (inverse) inequality (5.3) to control the consistency error $b_{\Gamma}^{T}(\mathbf{\Pi}_{h,0}^{\Delta t} \mathbf{u} - \mathbf{u}_{h}^{\Delta t}, \lambda - \mathcal{P}_{H,\Gamma}^{\Delta T} \lambda)$. This term can be made comparable to the other error terms in (5.4) by choosing $m$ and $s$ sufficiently large, assuming that the solution is sufficiently smooth. Alternatively, this term can be improved if a bound on $\|\nabla \cdot (\mathbf{u} - \mathbf{u}_{h}^{\Delta t})\|_{\Omega_i^T}$ is available, with the use of the normal trace inequality for $\mathbf{H}(\text{div}; \Omega_i)$ functions. In this case the factor $\Delta t^{-\frac{1}{2}}$ in the term $\Delta t^{-\frac{1}{2}}(h^{rl} + \Delta t^{rq})$ can also be avoided by using a suitable time-interpolant for the pressure. We present this argument in the next section in the special case of matching time steps between subdomains; then, additionally, the assumptions on quasi-uniform space and time meshes $\mathcal{T}_{h,i}$ and $\mathcal{T}_i^{\Delta t}$ can be avoided.
# 6 Control of the velocity divergence and improved error estimates

In this section we establish stability and error estimates for the velocity divergence, along with an improved error bound for the rest of the variables, as noted in Remark 5.2 above. For this section only, we make the following assumption on the temporal discretization:

$$ \forall i, j, \quad W_i^{\Delta t} = \Lambda_{ij}^{\Delta T} = W_j^{\Delta t}. \qquad (6.1) $$

In particular, we assume that all subdomains and mortar interfaces have the same time discretization, which we denote by $W^{\Delta t}$. Let $t^k$, $k = 0, \dots, N$, be the discrete times and let $q$ be the polynomial degree in $W^{\Delta t}$. In this section, consequently, $\Delta T = \Delta t$ and $s = q$.

We will utilize the Radau reconstruction operator $\mathcal{I}$ [18, 44], which satisfies, for any $\varphi(x, \cdot) \in W^{\Delta t}$, $\mathcal{I}\varphi(x, \cdot) \in H^1(0, T)$, $\mathcal{I}\varphi(x, \cdot)|_{(t^{k-1}, t^k)} \in P_{q+1}$, such that

$$ \int_{t^{k-1}}^{t^k} \partial_t \mathcal{I} \varphi \, d\varphi = \int_{t^{k-1}}^{t^k} \partial_t \varphi \, d\varphi + [\varphi]_{k-1} \phi_{k-1}^{+} \quad \forall \varphi(x, \cdot) \in W^{\Delta t}. $$

Recalling (3.3), this implies that

$$ \int_{0}^{T} \partial_t \mathcal{I} \varphi \, d\varphi = \int_{0}^{T} \tilde{\partial}_t \varphi \, d\varphi \quad \forall \varphi(x, \cdot) \in W^{\Delta t}. \qquad (6.2) $$

Then the second equation (3.4b) of the space-time mortar mixed method can be rewritten as

$$ (\partial_t \mathcal{I} p_h^{\Delta t}, w)_{\Omega^T} - b^T(\mathbf{u}_h^{\Delta t}, w) = (q, w)_{\Omega^T} \quad \forall w \in W_h^{\Delta t}. \qquad (6.3) $$

For notational convenience, for $\mathbf{v} \in \mathbf{V}$, let henceforth $\nabla_h \cdot \mathbf{v} \in L^2(\Omega)$ be such that $\forall i$, $(\nabla_h \cdot \mathbf{v})|_{\Omega_i} = \nabla \cdot (\mathbf{v}|_{\Omega_i})$.

## 6.1 Stability bound

We proceed with the stability bound for $\|\nabla_h \cdot \mathbf{u}_h^{\Delta t}\|_{L^2(0,T;L^2(\Omega))}$; under assumption (6.1), this complements the bound (4.21) of Theorem 4.4.

**Theorem 6.1 (Control of divergence).** Assume that condition (6.1) holds. Then, for the solution of the space-time mortar method (3.4), there exists a constant $C > 0$ independent of $h$, $H$, $\Delta t$, and $\Delta T$ such that

$$ \|\partial_t \mathcal{I} p_h^{\Delta t}\|_{\Omega^T} + \|\nabla_h \cdot \mathbf{u}_h^{\Delta t}\|_{\Omega^T} + \|\mathbf{u}_h^{\Delta t}\|_{\text{DG}} \le C(\|q\|_{\Omega^T} + \|\nabla \cdot K \nabla p_0\|_{\Omega}). $$

*Proof.* Since $\partial_t \mathcal{I} p_h^{\Delta t} \in W_h^{\Delta t}$ and $\nabla_h \cdot \mathbf{u}_h^{\Delta t} \in W_h^{\Delta t}$, (6.3) implies that $\partial_t \mathcal{I} p_h^{\Delta t} + \nabla_h \cdot \mathbf{u}_h^{\Delta t} = \mathcal{P}_h^{\Delta t} q$. Therefore,

$$ \|\partial_t \mathcal{I} p_h^{\Delta t}\|_{\Omega^T}^2 + \|\nabla_h \cdot \mathbf{u}_h^{\Delta t}\|_{\Omega^T}^2 + 2(\partial_t \mathcal{I} p_h^{\Delta t}, \nabla_h \cdot \mathbf{u}_h^{\Delta t})_{\Omega^T} = \|\mathcal{P}_h^{\Delta t} q\|_{\Omega^T}^2. \qquad (6.4) $$

To control the third term on the left, we note that (3.4a) implies that for each $k = 1, \dots, N$ and every $t \in [t^{k-1}, t^k]$, it holds that

$$ a(\mathbf{u}_h^{\Delta t}, \mathbf{v}) + b(\mathbf{v}, p_h^{\Delta t}) + b_{\Gamma}(\mathbf{v}, \lambda_H^{\Delta T}) = 0 \quad \forall \mathbf{v} \in \mathbf{V}_h. $$

Therefore, using that the initial data satisfy the above equation, see (4.18) and (4.15a), we have

$$ a^T(\partial_t \mathcal{I} \mathbf{u}_h^{\Delta t}, \mathbf{v}) + b^T(\mathbf{v}, \partial_t \mathcal{I} p_h^{\Delta t}) + b_{\Gamma}^T(\mathbf{v}, \partial_t \mathcal{I} \lambda_H^{\Delta T}) = 0 \quad \forall \mathbf{v} \in \mathbf{V}_h^{\Delta t}. \qquad (6.5) $$

Taking $\mathbf{v} = \mathbf{u}_h^{\Delta t}$ in (6.5) and $\mu = \partial_t \mathcal{I} \lambda_H^{\Delta T}$ in (3.4c) and combining the equations, we obtain

$$ -b^T(\mathbf{u}_h^{\Delta t}, \partial_t \mathcal{I} p_h^{\Delta t}) = a^T(\partial_t \mathcal{I} \mathbf{u}_h^{\Delta t}, \mathbf{u}_h^{\Delta t}) = \frac{1}{2}\|K^{-\frac{1}{2}}\mathbf{u}_h^{\Delta t}\|_{\text{DG}}^2 - \frac{1}{2}\|K^{-\frac{1}{2}}(\mathbf{u}_h^{\Delta t})_0\|_{\Omega}^2, \qquad (6.6) $$

using (6.2), (4.19), and (4.20) for the second equality. The assertion of the lemma follows by combining (6.4) and (6.6) and using (2.2), (4.18), and (4.16). $\square$
## 6.2 Improved a priori error estimate

In this section we utilize the control on $||\nabla_h \cdot \mathbf{u}_h^{\Delta t}||_{\Omega^T}$ to obtain error estimates that avoid the factors $h^{-\frac{1}{2}}$ and $\Delta t^{-\frac{1}{2}}$ that appear in the error estimate (5.4), together with the quasi-uniformity assumption on the space and time meshes $\mathcal{T}_{h,i}$ and $\mathcal{T}_i^{\Delta t}$. To this end, we will use an alternative time interpolant. Let $\tilde{\mathcal{P}}^{\Delta t}: H^1(0,T) \to W^{\Delta t}$ be such that, for any $\varphi \in H^1(0,T)$,

$$ \forall k = 1, \dots, N, \quad \int_{t^{k-1}}^{t^k} (\tilde{\mathcal{P}}^{\Delta t} \varphi - \varphi) w = 0 \quad \forall w \in P_{q-1}, \quad (\tilde{\mathcal{P}}^{\Delta t} \varphi)_k^- = \varphi(t^k). \quad (6.7) $$

Let us further set $(\tilde{\mathcal{P}}^{\Delta t} \varphi)_0^- = \varphi(0)$ and define the space-time interpolant

$$ \tilde{\mathcal{P}}_h^{\Delta t} = \mathcal{P}_h \times \tilde{\mathcal{P}}^{\Delta t}. \quad (6.8) $$

The following properties of $\tilde{\mathcal{P}}^{\Delta t}$ and $\tilde{\mathcal{P}}_h^{\Delta t}$ will be useful in the analysis.

**Lemma 6.2** (Time derivative orthogonality). For all $\varphi \in H^1(0,T)$ and $w \in W^{\Delta t}$,

$$ \int_0^T \partial_t \varphi w = \int_0^T \tilde{\partial}_t \tilde{\mathcal{P}}^{\Delta t} \varphi w. \quad (6.9) $$

Furthermore, for all $\varphi \in H^1(0,T)$ and $w \in W_h^{\Delta t}$,

$$ \int_0^T (\partial_t \varphi, w)_\Omega = \int_0^T (\tilde{\partial}_t \tilde{\mathcal{P}}_h^{\Delta t} \varphi, w)_\Omega. \quad (6.10) $$

*Proof.* Using (6.7), we write

$$ \begin{aligned} \int_0^T \partial_t \varphi w &= \sum_{k=1}^N \int_{t^{k-1}}^{t^k} \partial_t \varphi w = \sum_{k=1}^N \left( - \int_{t^{k-1}}^{t^k} \varphi \partial_t w + \varphi(t^k) w_k^- - \varphi(t^{k-1}) w_{k-1}^+ \right) \\ &= \sum_{k=1}^N \left( - \int_{t^{k-1}}^{t^k} \tilde{\mathcal{P}}^{\Delta t} \varphi \partial_t w + (\tilde{\mathcal{P}}^{\Delta t} \varphi)_k^- w_k^- - (\tilde{\mathcal{P}}^{\Delta t} \varphi)_{k-1}^- w_{k-1}^+ \right) \\ &= \sum_{k=1}^N \left( \int_{t^{k-1}}^{t^k} \partial_t \tilde{\mathcal{P}}^{\Delta t} \varphi w + [\tilde{\mathcal{P}}^{\Delta t} \varphi]_{k-1} w_{k-1}^+ \right) = \int_0^T \tilde{\partial}_t \tilde{\mathcal{P}}^{\Delta t} \varphi w, \end{aligned} $$

where we used the definition (3.3) in the last equality. This establishes (6.9). The identity (6.10) follows from (6.9), using that $\int_0^T (\partial_t \varphi, w)_\Omega = \int_0^T (\partial_t \mathcal{P}_h \varphi, w)_\Omega$. $\square$

Consider the space-time interpolants $\tilde{\mathbf{\Pi}}_{h,0}^{\Delta t} = \mathbf{\Pi}_{h,0} \times \tilde{\mathcal{P}}^{\Delta t}$ in $\mathbf{V}_{h,0}^{\Delta t}$ and $\tilde{\mathcal{P}}_h^{\Delta t}$ from (6.8) in $W_h^{\Delta t}$. Similarly to (5.1a) and (5.1b), they satisfy

$$ \|\psi - \tilde{\mathbf{\Pi}}_{h,0}^{\Delta t} \psi\|_{\Omega^T} \le C \sum_i \|\psi\|_{H^{r_q}(0,T;\mathbf{H}^{r_k}(\Omega_i))} (h^{r_k} + \Delta t^{r_q}) + C \|\psi\|_{H^{r_q}(0,T;\mathbf{H}^{r_k+\frac{1}{2}}(\Omega))} (h^{\tilde{r}_k} H^{\frac{1}{2}} + \Delta t^{r_q}), \\ 0 < r_k \le k+1, \quad 0 < \tilde{r}_k \le k+1, \quad 1 \le r_q \le q+1, \quad (6.11) $$

$$ \|\nabla \cdot (\psi - \tilde{\mathbf{\Pi}}_{h,0}^{\Delta t} \psi)\|_{\Omega_i^T} \le C \|\nabla \cdot \psi\|_{H^{r_q}(0,T;H^{r_l}(\Omega_i))} (h^{r_l} + \Delta t^{r_q}), \\ 0 \le r_l \le l+1, \quad 1 \le r_q \le q+1, \quad (6.12) $$

$$ \|\varphi - \tilde{\mathcal{P}}_h^{\Delta t} \varphi\|_{\Omega_i^T} \le C \|\varphi\|_{H^{r_q}(0,T;H^{r_l}(\Omega_i))} (h^{r_l} + \Delta t^{r_q}), \quad 0 \le r_l \le l+1, \quad 1 \le r_q \le q+1, \quad (6.13) $$
with similar bounds in $\|\cdot\|_{L^{\infty}(0,T;L^2(\Omega_i))}$, cf. (5.2). The use of $\tilde{\mathcal{P}}_h^{\Delta t}$ will allow us to avoid the term $\Delta t^{-\frac{1}{2}}$ in the error estimate.

We will also utilize the Scott-Zhang interpolant [48] $\tilde{\mathcal{P}}_{H,\Gamma}: H^1(\Gamma) \to \Lambda_{H,\Gamma} \cap C(\Gamma)$, which can be defined to preserve the trace on $\partial\Gamma$. Thus the function $\varphi - \tilde{\mathcal{P}}_{H,\Gamma}\varphi$ can be extended continuously by zero to $\partial\Omega \setminus \Gamma$ in the $H^{\frac{1}{2}}$-norm. Let us denote this extension by $E(\varphi - \tilde{\mathcal{P}}_{H,\Gamma}\varphi)$. Let $\tilde{\mathcal{P}}_{H,\Gamma}^{\Delta T} = \tilde{\mathcal{P}}_{H,\Gamma} \times \tilde{\mathcal{P}}^{\Delta t}$ be the space-time mortar interpolant. It has the approximation property [48]

$$ \begin{gathered} \|\varphi - \tilde{\mathcal{P}}_{H,\Gamma}^{\Delta T}\|_{t,\Gamma_{ij}^T} \le C \|\varphi\|_{H^{rs}(0,T;H^{rm}(\Gamma_{ij}))} (H^{rm-t} + \Delta T^{rs}), \tag{6.14} \\ 1 \le r_m \le m+1, \quad 1 \le r_s \le s+1, \quad 0 \le t \le 1. \end{gathered} $$

The use of $\tilde{\mathcal{P}}_{H,\Gamma}^{\Delta T}$, along with the well-known normal trace inequality, for all $\mathbf{v} \in \mathbf{H}(\text{div};\Omega_i)$, $\|\mathbf{v} \cdot \mathbf{n}_i\|_{-\frac{1}{2},\partial\Omega_i} \le C\|\mathbf{v}\|_{\text{div};\Omega_i}$, which implies

$$ \forall \mathbf{v} \in L^2(0, T; \mathbf{V}_i), \quad \|\mathbf{v} \cdot \mathbf{n}_i\|_{L^2(0,T;H^{-\frac{1}{2}}(\partial\Omega_i))} \le C \|\mathbf{v}\|_{L^2(0,T;\mathbf{V}_i)}, \quad (6.15) $$

will allow us to avoid the $h^{-\frac{1}{2}}$ term in the error estimate.

We are ready to prove the following error estimate that improves the bound (5.4) of Theorem 5.1 under assumption (6.1).

**Theorem 6.3** (Improved a priori error estimate). Assume that conditions (4.5a) and (6.1) hold and that the solution to (2.5) is sufficiently smooth. Then there exists a constant $C > 0$ independent of the mesh sizes $h$, $H$, $\Delta t$, and $\Delta T$, such that the solution to the space-time mortar MFE method (3.4) satisfies

$$ \begin{aligned} & \|\mathbf{u}(t^N) - (\mathbf{u}_h^{\Delta t})^-\|_{\Omega} + \|\mathbf{u} - \mathbf{u}_h^{\Delta t}\|_{\Omega^T} + \|\nabla_h \cdot (\mathbf{u} - \mathbf{u}_h^{\Delta t})\|_{\Omega^T} \\ & \qquad + \|p(t^N) - (p_h^{\Delta t})^-\|_{\Omega} + \|p - p_h^{\Delta t}\|_{\Omega^T} + \|\lambda - \lambda_H^{\Delta T}\|_{\Gamma^T} \\ & \le C \left( \sum_i \|\mathbf{u}\|_{W^{rq,\infty}(0,T;H^{rk}(\Omega_i))} (h^{rk} + \Delta t^{rq}) + \|\mathbf{u}\|_{W^{rq,\infty}(0,T;H^{\tilde{r}k+\frac{1}{2}}(\Omega))} (h^{\tilde{r}k} H^{\frac{1}{2}} + \Delta t^{rq}) \right. \\ & \qquad \left. + \sum_i \|p\|_{W^{rq,\infty}(0,T;H^{r_l}(\Omega_i))} (h^{r_l} + \Delta t^{rq}) + \sum_{i,j} \|\lambda\|_{H^{rs}(0,T;H^{rm}(\Gamma_{ij}))} (H^{rm-\frac{1}{2}} + \Delta T^{rs}) \right. \\ & \qquad \left. + \sum_i \|\mathbf{u}\|_{H^1(0,T;H^{rk}(\Omega_i))} h^{rk} + \|\mathbf{u}\|_{H^1(0,T;H^{\tilde{r}k+\frac{1}{2}}(\Omega))} h^{\tilde{r}k} H^{\frac{1}{2}} + \sum_{i,j} \|\lambda\|_{H^1(0,T;H^{rm}(\Gamma_{ij}))} H^{rm-\frac{1}{2}} \right), \\ & 0 < r_k \le k+1, \quad 0 < \tilde{r}_k \le k+1, \quad 1 \le r_q \le q+1, \quad 0 \le r_l \le l+1, \quad \frac{1}{2} \le r_m \le m+1, \quad 1 \le r_s \le s+1. \end{aligned} \quad (6.16) $$

*Proof.* Subtracting (5.5a)-(5.5b) from (2.5a)-(2.5b), we obtain the error equations, $\forall \mathbf{v} \in \mathbf{V}_{h,0}^{\Delta t}$, $w \in W_h^{\Delta t}$,

$$ a^T(\mathbf{u} - \mathbf{u}_h^{\Delta t}, \mathbf{v}) + b^T(\mathbf{v}, \tilde{\mathcal{P}}_h^{\Delta t}p - p_h^{\Delta t}) + b^T(\mathbf{v}, p - \tilde{\mathcal{P}}_h^{\Delta t}p) + b_\Gamma^T(\mathbf{v}, \lambda - \tilde{\mathcal{P}}_{H,\Gamma}^{\Delta T}\lambda) = 0, \quad (6.17a) $$

$$ (\partial_t p - \partial_t p_h^{\Delta t}, w)_{\Omega^T} - b^T(\tilde{\Pi}_{h,0}^{\Delta t}\mathbf{u} - \mathbf{u}_h^{\Delta t}, w) - b^T(\mathbf{u} - \tilde{\Pi}_{h,0}^{\Delta t}\mathbf{u}, w) = 0. \quad (6.17b) $$

We note the extra third terms in (6.17a) and (6.17b) compared to (5.6a) and (5.6b), which appear since $\tilde{\mathcal{P}}^{\Delta t}$ has orthogonality only for $P_{q-1}$. We take $\mathbf{v} = \tilde{\Pi}_{h,0}^{\Delta t}\mathbf{u} - \mathbf{u}_h^{\Delta t}$ and $w = \tilde{\mathcal{P}}_h^{\Delta t}p - p_h^{\Delta t}$ and sum the two equations, obtaining

$$ \begin{aligned} & a^T(\tilde{\Pi}_{h,0}^{\Delta t}\mathbf{u} - \mathbf{u}_h^{\Delta t}, \tilde{\Pi}_{h,0}^{\Delta t}\mathbf{u} - \mathbf{u}_h^{\Delta t}) + (\partial_t p - \partial_t p_h^{\Delta t}, \tilde{\mathcal{P}}_h^{\Delta t}p - p_h^{\Delta t})_{\Omega^T} \\ & = a^T(\tilde{\Pi}_{h,0}^{\Delta t}\mathbf{u} - \mathbf{u}, \tilde{\Pi}_{h,0}^{\Delta t}\mathbf{u} - \mathbf{u}_h^{\Delta t}) - b^T(\tilde{\Pi}_{h,0}^{\Delta t}\mathbf{u} - \mathbf{u}_h^{\Delta t}, p - \tilde{\mathcal{P}}_h^{\Delta t}p) \quad (6.18) \\ & \qquad - b_\Gamma^T(\tilde{\Pi}_{h,0}^{\Delta t}\mathbf{u} - \mathbf{u}_h^{\Delta t}, \lambda - \tilde{\mathcal{P}}_{H,\Gamma}^{\Delta T}\lambda) + b^T(\mathbf{u} - \tilde{\Pi}_{h,0}^{\Delta t}\mathbf{u}, \tilde{\mathcal{P}}_h^{\Delta t}p - p_h^{\Delta t}). \end{aligned} $$
For the second term on the left, using (6.10), (4.19), and (4.20), we write

$$ (\partial_t p - \tilde{\partial}_t p_h^{\Delta t}, \tilde{\mathcal{P}}_h^{\Delta t} p - p_h^{\Delta t})_{\Omega^T} = (\tilde{\partial}_t (\tilde{\mathcal{P}}_h^{\Delta t} p - p_h^{\Delta t}), \tilde{\mathcal{P}}_h^{\Delta t} p - p_h^{\Delta t})_{\Omega^T} \\ = \frac{1}{2} \|\tilde{\mathcal{P}}_h^{\Delta t} p - p_h^{\Delta t}\|_{\text{DG}}^2 - \frac{1}{2} \|\mathcal{P}_h p_0 - p_{h,0}\|_{\Omega}^2. \quad (6.19) $$

Combining (6.18) and (6.19), and using the Cauchy-Schwarz and Young inequalities, we obtain,

$$ \begin{aligned} & \|\tilde{\mathbf{\Pi}}_{h,0}^{\Delta t} \mathbf{u} - \mathbf{u}_h^{\Delta t}\|_{\Omega^T}^2 + \|\tilde{\mathcal{P}}_h^{\Delta t} p - p_h^{\Delta t}\|_{\text{DG}}^2 \\ & \le C \left( \|\tilde{\mathbf{\Pi}}_{h,0}^{\Delta t} \mathbf{u} - \mathbf{u}\|_{\Omega^T} \|\tilde{\mathbf{\Pi}}_{h,0}^{\Delta t} \mathbf{u} - \mathbf{u}_h^{\Delta t}\|_{\Omega^T} + \|\nabla_h \cdot (\tilde{\mathbf{\Pi}}_{h,0}^{\Delta t} \mathbf{u} - \mathbf{u}_h^{\Delta t})\|_{\Omega^T} \|p - \tilde{\mathcal{P}}_h^{\Delta t} p\|_{\Omega^T} + \|\mathcal{P}_h p_0 - p_{h,0}\|_{\Omega}^2 \right. \\ & \qquad \left. + \|\nabla_h \cdot (\mathbf{u} - \tilde{\mathbf{\Pi}}_{h,0}^{\Delta t} \mathbf{u})\|_{\Omega^T} \|\tilde{\mathcal{P}}_h^{\Delta t} p - p_h^{\Delta t}\|_{\Omega^T} \right. \\ & \qquad \left. + \sum_i \|(\tilde{\mathbf{\Pi}}_{h,0}^{\Delta t} \mathbf{u} - \mathbf{u}_h^{\Delta t}) \cdot \mathbf{n}_i\|_{L^2(0,T;H^{-\frac{1}{2}}(\partial\Omega_i))} \|E(\lambda - \tilde{\mathcal{P}}_{H,\Gamma}^{\Delta T} \lambda)\|_{L^2(0,T;H^{\frac{1}{2}}(\partial\Omega_i))} \right) \\ & \le \epsilon_1 \left( \|\tilde{\mathbf{\Pi}}_{h,0}^{\Delta t} \mathbf{u} - \mathbf{u}_h^{\Delta t}\|_{\Omega^T}^2 + \|\nabla_h \cdot (\tilde{\mathbf{\Pi}}_{h,0}^{\Delta t} \mathbf{u} - \mathbf{u}_h^{\Delta t})\|_{\Omega^T}^2 + \|\tilde{\mathcal{P}}_h^{\Delta t} p - p_h^{\Delta t}\|_{\Omega^T}^2 \right) \\ & \qquad + C \|\mathcal{P}_h p_0 - p_{h,0}\|_{\Omega}^2 \\ & \qquad + C_{\epsilon_1} \left( \|\tilde{\mathbf{\Pi}}_{h,0}^{\Delta t} \mathbf{u} - \mathbf{u}\|_{\Omega^T}^2 + \|\nabla_h \cdot (\mathbf{u} - \tilde{\mathbf{\Pi}}_{h,0}^{\Delta t} \mathbf{u})\|_{\Omega^T}^2 + \|p - \tilde{\mathcal{P}}_h^{\Delta t} p\|_{\Omega^T}^2 + \|\lambda - \tilde{\mathcal{P}}_{H,\Gamma}^{\Delta T} \lambda\|_{L^2(0,T;H^{\frac{1}{2}}(\Gamma))}^2 \right), \end{aligned} \quad (6.20) $$

where we used (6.15) in the last inequality. To complete the estimate, we bound the pressure, mortar, and divergence errors. The inf-sup condition for the weakly continuous velocity (4.7) and (6.17a) imply, using (6.15),

$$ \|\tilde{\mathcal{P}}_h^{\Delta t} p - p_h^{\Delta t}\|_{\Omega^T} \le C \left( \|\tilde{\mathbf{\Pi}}_{h,0}^{\Delta t} \mathbf{u} - \mathbf{u}_h^{\Delta t}\|_{\Omega^T} + \|\mathbf{u} - \tilde{\mathbf{\Pi}}_{h,0}^{\Delta t} \mathbf{u}\|_{\Omega^T} \right. \\ \left. + \|p - \tilde{\mathcal{P}}_h^{\Delta t} p\|_{\Omega^T} + \|\lambda - \tilde{\mathcal{P}}_{H,\Gamma}^{\Delta T} \lambda\|_{L^2(0,T;H^{\frac{1}{2}}(\Gamma))} \right). \quad (6.21) $$

To bound the mortar error using the mortar inf-sup condition (4.10), we subtract (3.4a) from (2.5a), to obtain the error equation

$$ a^T(\mathbf{u} - \mathbf{u}_h^{\Delta t}, \mathbf{v}) + b^T(\mathbf{v}, p - p_h^{\Delta t}) + b_{\Gamma}^T(\mathbf{v}, \tilde{\mathcal{P}}_{H,\Gamma}^{\Delta T} \lambda - \lambda_H^{\Delta T}) + b_{\Gamma}^T(\mathbf{v}, \lambda - \tilde{\mathcal{P}}_{H,\Gamma}^{\Delta T} \lambda) = 0 \quad \forall \mathbf{v} \in \mathbf{V}_h^{\Delta t}. \quad (6.22) $$

The mortar inf-sup condition (4.10) and (6.22) imply, using (6.15),

$$ \|\tilde{\mathcal{P}}_{H,\Gamma}^{\Delta T} \lambda - \lambda_H^{\Delta T}\|_{\Gamma^T} \le C \left( \|\tilde{\mathbf{\Pi}}_{h,0}^{\Delta t} \mathbf{u} - \mathbf{u}_h^{\Delta t}\|_{\Omega^T} + \|\mathbf{u} - \tilde{\mathbf{\Pi}}_{h,0}^{\Delta t} \mathbf{u}\|_{\Omega^T} \right. \\ \left. + \|\tilde{\mathcal{P}}_h^{\Delta t} p - p_h^{\Delta t}\|_{\Omega^T} + \|p - \tilde{\mathcal{P}}_h^{\Delta t} p\|_{\Omega^T} + \|\lambda - \tilde{\mathcal{P}}_{H,\Gamma}^{\Delta T} \lambda\|_{L^2(0,T;H^{\frac{1}{2}}(\Gamma))} \right). \quad (6.23) $$

Next, to bound on the divergence error, using (6.2) and (6.10), we rewrite (6.17b) as

$$ (\partial_t (\mathcal{I} \tilde{\mathcal{P}}_h^{\Delta t} p - \mathcal{I} p_h^{\Delta t}), w)_{\Omega^T} - b^T (\tilde{\mathbf{\Pi}}_{h,0}^{\Delta t} \mathbf{u} - \mathbf{u}_h^{\Delta t}, w) = b^T (\mathbf{u} - \tilde{\mathbf{\Pi}}_{h,0}^{\Delta t} \mathbf{u}, w), $$

concluding that $\partial_t (\mathcal{I} \tilde{\mathcal{P}}_h^{\Delta t} p - \mathcal{I} p_h^{\Delta t}) + \nabla_h \cdot (\tilde{\mathbf{\Pi}}_{h,0}^{\Delta t} \mathbf{u} - \mathbf{u}_h^{\Delta t}) = - \mathcal{P}_h^{\Delta t} (\nabla_h \cdot (\mathbf{u} - \tilde{\mathbf{\Pi}}_{h,0}^{\Delta t} \mathbf{u}))$. Therefore,

$$ \|\partial_t (\mathcal{I} \tilde{\mathcal{P}}_h^{\Delta t} p - \mathcal{I} p_h^{\Delta t})\|_{\Omega^T}^2 + \|\nabla_h \cdot (\tilde{\mathbf{\Pi}}_{h,0}^{\Delta t} \mathbf{u} - \mathbf{u}_h^{\Delta t})\|_{\Omega^T}^2 \\ + 2(\partial_t (\mathcal{I} \tilde{\mathcal{P}}_h^{\Delta t} p - \mathcal{I} p_h^{\Delta t}), \nabla_h \cdot (\tilde{\mathbf{\Pi}}_{h,0}^{\Delta t} \mathbf{u} - \mathbf{u}_h^{\Delta t}))_{\Omega^T} = \|\mathcal{P}_h^{\Delta t} (\nabla_h \cdot (\mathbf{u} - \tilde{\mathbf{\Pi}}_{h,0}^{\Delta t} \mathbf{u}))\|_{\Omega^T}^2. \quad (6.24) $$
To control the third term on the left, we note that, since the subdomain and mortar time discretizations are the same, the error equation (6.17a) holds for every $t \in [t^{k-1}, t^k]$, $k = 1, \dots, N$, with a test function $\mathbf{v} \in \mathbf{V}_{h,0}$. Therefore, similarly to (6.5), it implies that

$$a^T(\partial_t(\mathbf{u} - \mathcal{I}\mathbf{u}_h^{\Delta t}), \mathbf{v}) + b^T(\mathbf{v}, \partial_t(p - \mathcal{I}p_h^{\Delta t})) + b_{\Gamma}^T(\mathbf{v}, \partial_t(\lambda - \mathcal{I}\lambda_{H}^{\Delta T})) = 0 \quad \forall \mathbf{v} \in \mathbf{V}_h^{\Delta t}. \quad (6.25)$$

The first term is manipulated as

$$\begin{aligned} a^T(\partial_t(\mathbf{u} - \mathcal{I}\mathbf{u}_h^{\Delta t}), \mathbf{v}) &= a^T(\partial_t(\mathbf{u} - \mathbf{\Pi}_{h,0}\mathbf{u}), \mathbf{v}) + a^T(\partial_t(\mathbf{\Pi}_{h,0}\mathbf{u} - \mathcal{I}\mathbf{u}_h^{\Delta t}), \mathbf{v}) \\ &= a^T(\partial_t(\mathbf{u} - \mathbf{\Pi}_{h,0}\mathbf{u}), \mathbf{v}) + a^T(\tilde{\partial}_t(\tilde{\mathbf{\Pi}}_{h,0}^{\Delta t}\mathbf{u} - \mathbf{u}_h^{\Delta t}), \mathbf{v}), \end{aligned} \quad (6.26)$$

using (6.9) and (6.2) in the last equality. For the second term in (6.25) we write

$$b^T(\mathbf{v}, \partial_t(p - \mathcal{I}p_h^{\Delta t})) = b^T(\mathbf{v}, \partial_t(\mathcal{I}\tilde{\mathcal{P}}_h^{\Delta t}p - \mathcal{I}p_h^{\Delta t})), \quad (6.27)$$

using (6.10), (6.2), and the fact that $\nabla_h \cdot \mathbf{v} \in W_h^{\Delta t}$. The third term in (6.25) is manipulated as

$$\begin{aligned} b_{\Gamma}^{T}(\mathbf{v}, \partial_t(\lambda - \mathcal{I}\lambda_{H}^{\Delta T})) &= b_{\Gamma}^{T}(\mathbf{v}, \partial_t(\lambda - \tilde{\mathcal{P}}_{H,\Gamma}\lambda)) + b_{\Gamma}^{T}(\mathbf{v}, \partial_t(\tilde{\mathcal{P}}_{H,\Gamma}\lambda - \mathcal{I}\lambda_{H}^{\Delta T})) \\ &= b_{\Gamma}^{T}(\mathbf{v}, \partial_t(\lambda - \tilde{\mathcal{P}}_{H,\Gamma}\lambda)) + b_{\Gamma}^{T}(\mathbf{v}, \partial_t(\mathcal{I}\tilde{\mathcal{P}}_{H,\Gamma}^{\Delta t}\lambda - \mathcal{I}\lambda_{H}^{\Delta T})), \end{aligned} \quad (6.28)$$

using (6.9) and (6.2) in the last equality. Now, combining (6.25)–(6.28), taking $\mathbf{v} = \tilde{\mathbf{\Pi}}_{h,0}^{\Delta t}\mathbf{u} - \mathbf{u}_h^{\Delta t} \in \mathbf{V}_{h,0}^{\Delta t}$, and using that $\partial_t(\mathcal{I}\tilde{\mathcal{P}}_{H,\Gamma}^{\Delta t}\lambda - \mathcal{I}\lambda_{H}^{\Delta T}) \in \Lambda_{H}^{\Delta T}$ to deduce that the last term in (6.28) is zero, we obtain

$$\begin{aligned} (\partial_t(\mathcal{I}\tilde{\mathcal{P}}_h^{\Delta t}p - \mathcal{I}p_h^{\Delta t}), \nabla_h \cdot (\tilde{\mathbf{\Pi}}_{h,0}^{\Delta t}\mathbf{u} - \mathbf{u}_h^{\Delta t}))_{\Omega^T} &= a^T(\tilde{\partial}_t(\tilde{\mathbf{\Pi}}_{h,0}^{\Delta t}\mathbf{u} - \mathbf{u}_h^{\Delta t}), \tilde{\mathbf{\Pi}}_{h,0}^{\Delta t}\mathbf{u} - \mathbf{u}_h^{\Delta t}) \\ &\quad + a^T(\partial_t(\mathbf{u} - \mathbf{\Pi}_{h,0}\mathbf{u}), \tilde{\mathbf{\Pi}}_{h,0}^{\Delta t}\mathbf{u} - \mathbf{u}_h^{\Delta t}) + b_{\Gamma}^T(\tilde{\mathbf{\Pi}}_{h,0}^{\Delta t}\mathbf{u} - \mathbf{u}_h^{\Delta t}, \partial_t(\lambda - \tilde{\mathcal{P}}_{H,\Gamma}\lambda)). \end{aligned} \quad (6.29)$$

Combining (6.24) and (6.29) and using (4.19) and the Cauchy–Schwarz and Young inequalities, we arrive at

$$\begin{aligned} & \|\partial_t(\mathcal{I}\tilde{\mathcal{P}}_h^{\Delta t}p - \mathcal{I}p_h^{\Delta t})\|_{\Omega^T}^2 + \|\nabla_h \cdot (\tilde{\mathbf{\Pi}}_{h,0}^{\Delta t}\mathbf{u} - \mathbf{u}_h^{\Delta t})\|_{\Omega^T}^2 + \|\tilde{\mathbf{\Pi}}_{h,0}^{\Delta t}\mathbf{u} - \mathbf{u}_h^{\Delta t}\|_{\text{DG}}^2 \\ & \le \epsilon_2 \left( \|\tilde{\mathbf{\Pi}}_{h,0}^{\Delta t}\mathbf{u} - \mathbf{u}_h^{\Delta t}\|_{\Omega^T}^2 + \|\nabla_h \cdot (\tilde{\mathbf{\Pi}}_{h,0}^{\Delta t}\mathbf{u} - \mathbf{u}_h^{\Delta t})\|_{\Omega^T}^2 \right) + \|\nabla_h \cdot (\mathbf{u} - \tilde{\mathbf{\Pi}}_{h,0}^{\Delta t}\mathbf{u})\|_{\Omega^T}^2 \\ & \quad + C_{\epsilon_2} \left( \|\partial_t(\mathbf{u} - \mathbf{\Pi}_{h,0}\mathbf{u})\|_{\Omega^T}^2 + \|\partial_t(\lambda - \tilde{\mathcal{P}}_{H,\Gamma}\lambda)\|_{L^2(0,T;H^{\frac{1}{2}}(\Gamma))}^2 \right) + C \|\mathbf{\Pi}_{h,0}\mathbf{u}_0 - \mathbf{u}_{h,0}\|_{\Omega}^2. \end{aligned} \quad (6.30)$$

Combining (6.20), (6.21), (6.23), and (6.30), taking $\epsilon_1$ and $\epsilon_2$ small enough, we obtain

$$\begin{aligned} & \|\tilde{\mathbf{\Pi}}_{h,0}^{\Delta t}\mathbf{u} - \mathbf{u}_h^{\Delta t}\|_{\Omega^T}^2 + \|\nabla_h \cdot (\tilde{\mathbf{\Pi}}_{h,0}^{\Delta t}\mathbf{u} - \mathbf{u}_h^{\Delta t})\|_{\Omega^T}^2 + \|\tilde{\mathbf{\Pi}}_{h,0}^{\Delta t}\mathbf{u} - \mathbf{u}_h^{\Delta t}\|_{\text{DG}}^2 \\ & \quad + \|\tilde{\mathcal{P}}_h^{\Delta t}p - p_h^{\Delta t}\|_{\Omega^T}^2 + \|\tilde{\mathcal{P}}_h^{\Delta t}p - p_h^{\Delta t}\|_{\text{DG}}^2 + \|\tilde{\mathcal{P}}_{H,\Gamma}^{\Delta T}\lambda - \lambda_H^{\Delta T}\|_{\Gamma^T}^2 \\ & \le C \left( \|\mathbf{u} - \tilde{\mathbf{\Pi}}_{h,0}^{\Delta t}\mathbf{u}\|_{L^2(0,T;\mathbf{v})}^2 + \|p - \tilde{\mathcal{P}}_h^{\Delta t}p\|_{\Omega^T}^2 + \|\lambda - \tilde{\mathcal{P}}_{H,\Gamma}^{\Delta T}\lambda\|_{L^2(0,T;H^{\frac{1}{2}}(\Gamma))}^2 + \|\partial_t(\mathbf{u} - \mathbf{\Pi}_{h,0}\mathbf{u})\|_{\Omega^T}^2 \right. \\ & \quad \left. + \|\partial_t(\lambda - \tilde{\mathcal{P}}_{H,\Gamma}\lambda)\|_{L^2(0,T;H^{\frac{1}{2}}(\Gamma))}^2 + \|\mathbf{\Pi}_{h,0}\mathbf{u}_0 - \mathbf{u}_{h,0}\|_{\Omega}^2 + \|\mathcal{P}_h p_0 - p_{h,0}\|_{\Omega}^2 \right). \end{aligned}$$

Finally, using the triangle inequality, the approximation properties (6.11)–(6.14), and the initial error bound (4.17), we arrive at (6.16). We remark that in the final bound, we have kept only the term at time $t^N$ from the norm $\| \cdot \|_{\text{DG}}$, since the approximation error in the full $\| \cdot \|_{\text{DG}}$ norm involves a $\Delta t^{-\frac{1}{2}}$ factor. $\square$

**Remark 6.4** (Significance of the improved error estimate). The error estimate (6.16) avoids the factors $h^{-\frac{1}{2}}$ and $\Delta t^{-\frac{1}{2}}$, which appeared in the earlier bound (5.4), and thus provides optimal order of convergence for all variables. Moreover, it provides a bound on the velocity divergence error. To the best of the authors’ knowledge, such result has not been established in the literature for space-time mixed finite element methods with a DG time discretization, even on a single domain.
# 7 Reduction to an interface problem

In this section we combine the time-dependent Steklov-Poincaré operator approach from [33] with the mortar domain decomposition algorithm from [3, 5] to reduce the global problem (3.4) to a space-time interface problem.

## 7.1 Decomposition of the solution

Consider a decomposition of the solution to (3.4) in the form

$$ \mathbf{u}_h^{\Delta t} = \mathbf{u}_h^{\Delta t,*}(\lambda_H^{\Delta T}) + \bar{\mathbf{u}}_h^{\Delta t}, \quad p_h^{\Delta t} = p_h^{\Delta t,*}(\lambda_H^{\Delta T}) + \bar{p}_h^{\Delta t}. \qquad (7.1) $$

Here, $\bar{\mathbf{u}}_h^{\Delta t} \in \mathbf{V}_h^{\Delta t}$, $\bar{p}_h^{\Delta t} \in W_h^{\Delta t}$ are such that for each $\Omega_i^T$, $(\bar{\mathbf{u}}_h^{\Delta t}|_{\Omega_i^T} \in \mathbf{V}_{h,i}^{\Delta t}, \bar{p}_h^{\Delta t}|_{\Omega_i^T} \in W_{h,i}^{\Delta t})$ is the solution to the space-time subdomain problem in $\Omega_i^T$ with zero Dirichlet data on the space-time interfaces and the prescribed source term, initial data, and boundary data on the external boundary:

$$ a_i^T(\bar{\mathbf{u}}_h^{\Delta t}, \mathbf{v}) + b_i^T(\mathbf{v}, \bar{p}_h^{\Delta t}) = 0 \quad \forall \mathbf{v} \in \mathbf{V}_{h,i}^{\Delta t} \qquad (7.2a) $$

$$ (\tilde{\partial}_t \bar{p}_h^{\Delta t}, w)_{\Omega_i^T} - b_i^T(\bar{\mathbf{u}}_h^{\Delta t}, w) = (q, w)_{\Omega_i^T} \quad \forall w \in W_{h,i}^{\Delta t}. \qquad (7.2b) $$

Furthermore, for a given $\mu \in \Lambda_H^{\Delta T}$, $\mathbf{u}_h^{\Delta t,*}(\mu) \in \mathbf{V}_h^{\Delta t}$, $p_h^{\Delta t,*}(\mu) \in W_h^{\Delta t}$ are such that for each $\Omega_i^T$, $(\mathbf{u}_h^{\Delta t,*}(\mu)|_{\Omega_i^T} \in \mathbf{V}_{h,i}^{\Delta t}, p_h^{\Delta t,*}(\mu)|_{\Omega_i^T} \in W_{h,i}^{\Delta t})$ is the solution to the space-time subdomain problem in $\Omega_i^T$ with Dirichlet data $\mu$ on the space-time interfaces and zero source term, initial data, and boundary data on the external boundary:

$$ a_i^T(\mathbf{u}_h^{\Delta t,*}(\mu), \mathbf{v}) + b_i^T(\mathbf{v}, p_h^{\Delta t,*}(\mu)) = -\langle \mathbf{v} \cdot \mathbf{n}_i, \mu \rangle_{\Gamma_i^T} \quad \forall \mathbf{v} \in \mathbf{V}_{h,i}^{\Delta t}, \qquad (7.3a) $$

$$ (\tilde{\partial}_t p_h^{\Delta t,*}(\mu), w)_{\Omega_i^T} - b_i^T(\mathbf{u}_h^{\Delta t,*}(\mu), w) = 0 \quad \forall w \in W_{h,i}^{\Delta t}. \qquad (7.3b) $$

Note that both (7.2) and (7.3) are posed in the individual space-time subdomains $\Omega_i^T$ and can thus be solved in parallel (on the entire space-time subdomains $\Omega_i^T$, without any synchronization on time steps). It is easy to check that (3.4) is equivalent to solving the space-time interface problem: find $\lambda_H^{\Delta T} \in \Lambda_H^{\Delta T}$ such that

$$ -b_{\Gamma}^{T}(\mathbf{u}_{h}^{\Delta t,*}(\lambda_{H}^{\Delta T}), \mu) = b_{\Gamma}^{T}(\bar{\mathbf{u}}_{h}^{\Delta t}, \mu) \quad \forall \mu \in \Lambda_{H}^{\Delta T}, \qquad (7.4) $$

and obtaining $\mathbf{u}_h^{\Delta t}$ and $p_h^{\Delta t}$ from (7.1)–(7.3).

## 7.2 Space-time Steklov-Poincaré operator

The above problem can be written in an operator form: find $\lambda_H^{\Delta T} \in \Lambda_H^{\Delta T}$ such that

$$ S \lambda_H^{\Delta T} = g, \qquad (7.5) $$

where $S: \Lambda_H^{\Delta T} \to \Lambda_H^{\Delta T}$ is the space-time Steklov-Poincaré operator defined as

$$ \langle S\lambda, \mu \rangle_{\Gamma^T} = \sum_i \langle S_i\lambda, \mu \rangle_{\Gamma_i^T}, \quad \langle S_i\lambda, \mu \rangle_{\Gamma_i^T} = -\langle \mathbf{u}_h^{\Delta t,*}(\lambda) \cdot \mathbf{n}_i, \mu \rangle_{\Gamma_i^T} \quad \forall \lambda, \mu \in \Lambda_H^{\Delta T}, \qquad (7.6) $$

and $g \in \Lambda_H^{\Delta T}$ is defined as $\langle g, \mu \rangle_{\Gamma^T} = b_{\Gamma}(\bar{\mathbf{u}}_h^{\Delta t}, \mu) \forall \mu \in \Lambda_H^{\Delta T}$.
**Lemma 7.1** (Space-time Steklov-Poincaré operator). Assume that conditions (4.5) hold. Then the operator S defined in (7.6) is positive definite.

*Proof.* For $\lambda, \mu \in \Lambda_H^{\Delta T}$, consider (7.3a) with data $\mu$ and test function $v = u_h^{\Delta t,*}(\lambda)$. This implies, using (7.6),

$$ \begin{aligned} \langle S\lambda, \mu \rangle_{\Gamma^T} &= a^T(\mathbf{u}_h^{\Delta t,*}(\mu), \mathbf{u}_h^{\Delta t,*}(\lambda)) + b^T(\mathbf{u}_h^{\Delta t,*}(\lambda), p_h^{\Delta t,*}(\mu)) \\ &= a^T(\mathbf{u}_h^{\Delta t,*}(\mu), \mathbf{u}_h^{\Delta t,*}(\lambda)) + (\tilde{\partial}_t p_h^{\Delta t,*}(\lambda), p_h^{\Delta t,*}(\mu))_{\Omega^T}, \end{aligned} \quad (7.7) $$

where we have used (7.3b) with data $\lambda$ and test function $p_h^{\Delta t,*}(\mu)$ in the second equality. Lemma 4.3 together with $p_h^{\Delta t,*}(\mu)(x, 0) = 0$ (recall that zero initial data is supposed in (7.3)) imply that

$$ \langle S\mu, \mu \rangle_{\Gamma^T} \geq a^T(\mathbf{u}_h^{\Delta t,*}(\mu), \mathbf{u}_h^{\Delta t,*}(\mu)) \geq 0 \quad \forall \mu \in \Lambda_H^{\Delta T}, \quad (7.8) $$

hence S is positive semi-definite. Assume that $\langle S\mu, \mu \rangle_{\Gamma^T} = 0$. Then $\mathbf{u}_h^{\Delta t,*}(\mu) = 0$. The inf-sup condition for the weakly continuous velocity (4.7) and (7.3a) imply $p_h^{\Delta t,*}(\mu) = 0$. Then the mortar inf-sup condition (4.10) and (7.3a) imply $\mu = 0$, thus S is positive definite. $\square$

Due to Lemma 7.1, GMRES can be employed to solve the interface problem (7.5). On each GMRES iteration, the dominant computational cost is the evaluation of the action of S, which requires solving space-time problems with prescribed Dirichlet interface data in each individual space-time subdomain $\Omega_i \times (0, T)$. The following result can be used to provide a bound on the number of GMRES iterations.

**Theorem 7.2** (Spectral bound). Assume that conditions (4.5) hold. Let $T_{h,i}$ be quasi-uniform and $h \le Ch_i \forall i$. Then there exist positive constants $C_0$ and $C_1$ independent of the mesh sizes h, H, $\Delta t$, and $\Delta T$, such that

$$ \forall \mu \in \Lambda_H^{\Delta T}, \quad C_0 \|\mu\|_{\Gamma^T}^2 \le \langle S\mu, \mu \rangle_{\Gamma^T} \le C_1 h^{-1} \|\mu\|_{\Gamma^T}^2. \quad (7.9) $$

*Proof.* Using (7.6), the Cauchy-Schwarz inequality, (5.3), and $h \le Ch_i$, we obtain

$$ \langle S_i \mu, \mu \rangle_{\Gamma_i^T} \le \|\mathbf{u}_h^{\Delta t,*}(\mu) \cdot \mathbf{n}_i\|_{\Gamma_i^T} \|\mu\|_{\Gamma_i^T} \le C h^{-\frac{1}{2}} \|\mathbf{u}_h^{\Delta t,*}(\mu)\|_{\Omega_i^T} \|\mu\|_{\Gamma_i^T} \le C h^{-\frac{1}{2}} \langle S_i \mu, \mu \rangle_{\Gamma_i^T}^{\frac{1}{2}} \|\mu\|_{\Gamma_i^T}, $$

where we used (7.8), which is also valid on each $\Omega_i^T$, in the last inequality. This implies the upper bound in (7.9).

To prove the lower bound in (7.9), we consider the set of auxiliary subdomain problems (4.11) with data $\mu$. Let $v_i = \Pi_{h,i}^{\Delta t} \psi_i$ and recall that $v_i \cdot n_i = Q_{h,i}^{\Delta t} \mu$. Using (4.6) and (7.3a), we have

$$ \begin{aligned} \|\mu\|_{\Gamma^T}^2 &\le C \sum_i \langle Q_{h,i}^{\Delta t} \mu, Q_{h,i}^{\Delta t} \mu \rangle_{\Gamma_i^T} = C \sum_i \langle Q_{h,i}^{\Delta t} \mu, \mu \rangle_{\Gamma_i^T} = C \sum_i \langle v_i \cdot n_i, \mu \rangle_{\Gamma_i^T}, \\ &= -C \sum_i \left( a_i^T (\mathbf{u}_h^{\Delta t,*}(\mu), v_i) + b_i^T (v_i, p_h^{\Delta t,*}(\mu)) \right) \\ &\le C \sum_i \left( \|\mathbf{u}_h^{\Delta t,*}(\mu)\|_{\Omega_i^T}^2 + \|p_h^{\Delta t,*}(\mu)\|_{\Omega_i^T}^2 \right)^{\frac{1}{2}} \|v_i\|_{L^2(0,T; v_i)} \\ &\le C \left\{ \sum_i \|\mathbf{u}_h^{\Delta t,*}(\mu)\|_{\Omega_i^T}^2 \right\}^{\frac{1}{2}} \left\{ \sum_i \|\mu\|_{\Gamma_i^T}^2 \right\}^{\frac{1}{2}} \\ &\le C \langle S\mu, \mu \rangle_{\Gamma^T}^{\frac{1}{2}} \|\mu\|_{\Gamma^T}. \end{aligned} $$
In the next to last inequality above, we used the Cauchy–Schwarz inequality together with the inf-sup condition (4.7) and (7.3a) to bound $||p_h^{\Delta t,*}(\mu)||_{\Omega^T}$ and the elliptic regularity (4.14) to bound $||v_i||_{L^2(0,T;V_i)}$. In the last inequality we used (7.8). This concludes the proof. □

## 7.3 GMRES convergence through the field-of-values estimates

Theorem 7.2 leads to convergence estimates for solving the interface problem (7.5) with GMRES. In [17, Theorem 3.3], a bound is shown for the k-th residual $r_k$ of the generalized conjugate residual method for solving a system with a positive definite matrix $S \in \mathbb{R}^{n \times n}$, which also applies to GMRES. It can be stated in terms of angle $\beta \in [0, \pi/2)$, see [9]:

$$ ||r_k|| \le \sin^k(\beta) ||r_0||, \quad \text{where } \cos(\beta) = \frac{\lambda_{\min}((S + S^T)/2)}{||S||}, \qquad (7.10) $$

where $||\cdot||$ denotes the Euclidean vector norm and the induced matrix norm. The quantities in (7.10) can be interpreted in terms of the field-of-values of $S$, defined as

$$ W(S) = \{ \zeta^T S \zeta : \zeta \in \mathbb{C}^n, ||\zeta|| = 1 \}. $$

It is known (see [36, Chapter 15]) that $W(S)$ is a compact and convex set in the complex plane that contains (but is usually much larger than) the eigenvalues of $S$. Because $S$ is positive definite, $0 \notin W(S)$, and because $S$ is real, the smallest eigenvalue of the symmetric part of $S$ is actually the distance from $0$ to $W(S)$, so that the angle $\beta$ can be improved to, see [9] or [29, Theorem 2.2.2],

$$ \cos(\beta) = \frac{\mathrm{dist}(0, W(S))}{||S||}. $$

The above bound, together with inequalities (7.9) obtained in Theorem 7.2, imply that the reduction in the k-th GMRES residual for solving the interface problem (7.5) is bounded by

$$ ||r_k|| \le \left( \sqrt{1 - (C_0/C_1)^2 h^2} \right)^k ||r_0||. \qquad (7.11) $$

A similar inequality, allowing for an explicit preconditioning matrix, has been obtained in [49].

## 8 Numerical results

In this section, we present several numerical results obtained with the space-time mortar method developed in Section 3.2, illustrating the convergence rates and other theoretical results obtained in the previous sections. The method is implemented using the deal.II finite element package [7].

In all the examples, we consider two-dimensional spatial domains and take the mixed finite element spaces $V_{h,i} \times W_{h,i}$ on the spatial subdomain $\Omega_i$ to be the lowest order Raviart–Thomas pair $RT_0 \times DGQ_0$ (i.e., $k = l = 0$) on quadrilateral meshes [11], where $DGQ_r$ denotes the space of discontinuous piecewise polynomials of degree up to $r$ in each variable. Combining this with the lowest-order DG (backward Euler, $q = 0$) for time discretization on the mesh $T_i^{\Delta t}$ gives us a space-time mixed finite element space $V_{h,i}^{\Delta t} \times W_{h,i}^{\Delta t}$ in $\Omega_i^T$ as detailed in Section 3.1. We test two different choices for the mortar finite element space $\Lambda_{H,ij}^{\Delta T}$ on the space-time interface mesh $T_{H,ij}^{\Delta T}$, with $\Delta T$ suitably chosen as a function of $\Delta t$, depending on the mortar space polynomial degree. These are discontinuous bilinear $DGQ_1$ ($m = s = 1$) and biquadratic $DGQ_2$ ($m = s = 2$) mortars.

For solving the interface problem identified in Section 7.2, we have implemented the GMRES algorithm without preconditioner. Developing a preconditioner for the iterative solver, which could significantly reduce the number of iterations, and its theoretical analysis is a subject of future research.
<table><caption>Table 1: Example 1, mesh size and number of degrees of freedom</caption><thead><tr><th rowspan="2">Ref.<br>No.</th><th colspan="3">Ω<sub>1</sub><sup>T</sup></th><th colspan="3">Ω<sub>2</sub><sup>T</sup></th><th colspan="3">Ω<sub>3</sub><sup>T</sup></th><th colspan="3">Ω<sub>4</sub><sup>T</sup></th><th colspan="3">Γ<sup>T</sup>(m = 1)</th><th colspan="3">Γ<sup>T</sup>(m = 2)</th></tr><tr><th>n<sub>1</sub></th><th>N<sub>1</sub></th><th>#DoF</th><th>n<sub>2</sub></th><th>N<sub>2</sub></th><th>#DoF</th><th>n<sub>3</sub></th><th>N<sub>3</sub></th><th>#DoF</th><th>n<sub>4</sub></th><th>N<sub>4</sub></th><th>#DoF</th><th>n<sub>Γ</sub></th><th>N<sub>Γ</sub></th><th>#DoF</th><th>n<sub>Γ</sub></th><th>N<sub>Γ</sub></th><th>#DoF</th></tr></thead><tbody><tr><td>0</td><td>3</td><td>3</td><td>33</td><td>2</td><td>2</td><td>16</td><td>4</td><td>4</td><td>56</td><td>3</td><td>3</td><td>33</td><td>1</td><td>1</td><td>16</td><td>1</td><td>1</td><td>36</td></tr><tr><td>1</td><td>6</td><td>6</td><td>120</td><td>4</td><td>4</td><td>56</td><td>8</td><td>8</td><td>208</td><td>6</td><td>6</td><td>120</td><td>2</td><td>2</td><td>64</td><td></td><td></td><td></td></tr><tr><td>2</td><td>12</td><td>12</td><td>456</td><td>8</td><td>8</td><td>208</td><td>16</td><td>16</td><td>800</td><td>12</td><td>12</td><td>456</td><td>4</td><td>4</td><td>256</td><td>2</td><td>2</td><td>144</td></tr><tr><td>3</td><td>24</td><td>24</td><td>1776</td><td>16</td><td>16</td><td>800</td><td>32</td><td>32</td><td>3136</td><td>24</td><td>24</td><td>1776</td><td>8</td><td>8</td><td>1024</td><td></td><td></td><td></td></tr><tr><td>4</td><td>48</td><td>48</td><td>7008</td><td>32</td><td>32</td><td>3136</td><td>64</td><td>64</td><td>12416</td><td>48</td><td>48</td><td>7008</td><td>16</td><td>16</td><td>4096</td><td>4</td><td>4</td><td>576</td></tr></tbody></table>

<table><caption>Table 2: Example 1, convergence with bilinear mortars</caption><thead><tr><th rowspan="2">Ref.</th><th rowspan="2"># GMRES</th><th colspan="2">‖u - u<sub>h</sub><sup>Δt</sup>‖<sub>L<sup>2</sup>(0,T;L<sup>2</sup>(Ω))</sub></th><th colspan="2">‖p - p<sub>h</sub><sup>Δt</sup>‖<sub>DG</sub></th><th colspan="2">‖p - p<sub>h</sub><sup>Δt</sup>‖<sub>L<sup>2</sup>(0,T;W)</sub></th><th colspan="2">‖λ - λ<sub>H</sub><sup>ΔT</sup>‖<sub>L<sup>2</sup>(0,T;Λ<sub>H</sub>)</sub></th></tr><tr><th>Rate</th><th>Rate</th><th>Rate</th><th>Rate</th><th>Rate</th><th>Rate</th><th>Rate</th><th>Rate</th></tr></thead><tbody><tr><td>0</td><td>11</td><td>6.50e-01</td><td>0.84</td><td>1.21e+00</td><td>0.75</td><td>7.91e-01</td><td>0.73</td><td>7.98e-01</td><td>0.64</td></tr><tr><td>1</td><td>23</td><td>-1.06</td><td>3.63e-01</td><td>7.21e-01</td><td>1.18</td><td>4.76e-01</td><td>0.95</td><td>5.11e-01</td><td>1.13</td></tr><tr><td>2</td><td>39</td><td>-0.76</td><td>1.74e-01</td><td>3.19e-01</td><td>1.18</td><td>2.46e-01</td><td>0.95</td><td>2.34e-01</td><td>1.13</td></tr><tr><td>3</td><td>59</td><td>-0.60</td><td>8.63e-02</td><td>1.02</td><td>1.46e-01</td><td>1.13</td><td>1.25e-01</td><td>0.98</td><td>1.20e-01</td><td>0.96</td></tr><tr><td>4</td><td>86</td><td>-0.54</td><td>4.29e-02</td><td>1.01</td><td>6.93e-02</td><td>1.08</td><td>6.25e-02</td><td>1.00</td><td>6.11e-02</td><td>0.97</td></tr></tbody></table>

## 8.1 Example 1: convergence test

In this example, we solve the parabolic problem (2.1) with a known solution to verify the accuracy of the space-time mortar method. We also discuss the correspondence of the number of interface GMRES iterations to the theoretical estimate. We further compare the accuracy and computational cost with $DGQ_1$ and $DGQ_2$ mortar spaces. We use the known pressure function $p(x, y, t) = \sin(8t) \sin(11x) \cos(11y - \frac{\pi}{4})$ along with permeability $K = I_{2 \times 2}$ to determine the right-hand side $q$ in (2.1) and impose Dirichlet boundary condition and initial condition on the space-time domain $\Omega^T = (0,1)^2 \times (0,0.5)$.

We partition the space domain $\Omega$ into four identical squares $\Omega_i$, $i = 1, 2, 3, 4$, with $\Omega_1 = (0,0.5) \times (0,0.5)$, $\Omega_2 = (0.5,1) \times (0,0.5)$, $\Omega_3 = (0,0.5) \times (0.5,1)$, and $\Omega_4 = (0.5,1) \times (0.5,1)$. The space-time domain $\Omega^T$ is correspondingly partitioned into four space-time subdomains $\Omega_i^T$, $i = 1, 2, 3, 4$. We start with an initial grid for each $\Omega_i^T$ and $\Gamma_{ij}^T$ and refine it successively 4 times to test the convergence rate of the solutions with respect to the actual known solution. The subdomains $\Omega_i^T$ maintain a checkerboard non-matching mesh structure throughout the refinement cycles. In particular, let $n_i$ be number of elements in either the $x$ or $y$-directions and let $N_i$ be the number of elements in the $t$-direction in subdomain $\Omega_i^T$. The initial grids are chosen as $n_1 = N_1 = 3$, $n_2 = N_2 = 2$, $n_3 = N_3 = 4$, and $n_4 = N_4 = 3$, see Table 1. Note that $h = \Delta t$. In the case of bilinear mortars, we maintain $H = 2h$ and $\Delta T = 2\Delta t$, halving the mesh sizes on each refinement cycle. For biquadratic mortars, we start with $H = 2h$ and $\Delta T = 2\Delta t$, but refine the mortar mesh only every other time to maintain $H = \sqrt{h}$ and $\Delta T = \sqrt{\Delta t}$. The coarser mesh on $\Gamma^T$ in the biquadratic

<table><caption>Table 3: Example 1, convergence with biquadratic mortars</caption><thead><tr><th rowspan="2">Ref.</th><th rowspan="2"># GMRES</th><th colspan="2">‖u - u<sub>h</sub><sup>Δt</sup>‖<sub>L<sup>2</sup>(0,T;L<sup>2</sup>(Ω))</sub></th><th colspan="2">‖p - p<sub>h</sub><sup>Δt</sup>‖<sub>DG</sub></th><th colspan="2">‖p - p<sub>h</sub><sup>Δt</sup>‖<sub>L<sup>2</sup>(0,T;W)</sub></th><th colspan="2">‖λ - λ<sub>H</sub><sup>ΔT</sup>‖<sub>L<sup>2</sup>(0,T;Λ<sub>H</sub>)</sub></th></tr><tr><th>Rate</th><th>Rate</th><th>Rate</th><th>Rate</th><th>Rate</th><th>Rate</th><th>Rate</th><th>Rate</th></tr></thead><tbody><tr><td>0</td><td>18</td><td>6.81e-01</td><td>1.00</td><td>1.35e+00</td><td>0.97</td><td>8.39e-01</td><td>0.87</td><td>2.13e+00</td><td>1.46</td></tr><tr><td>2</td><td>34</td><td>-0.46</td><td>1.70e-01</td><td>3.51e-01</td><td>0.97</td><td>2.51e-01</td><td>0.87</td><td>2.82e-01</td><td>1.46</td></tr><tr><td>4</td><td>57</td><td>-0.37</td><td>4.48e-02</td><td>0.96</td><td>8.59e-02</td><td>1.02</td><td>6.59e-02</td><td>0.96</td><td>9.20e-02</td><td>0.81</td></tr></tbody></table>
<figure><img src="image_3.png" /><figcaption>(a) On the whole space-time domain $\Omega^T$</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(b) On $\Omega_1^T \cup \Omega_4^T$</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(c) On $\Omega_2^T \cup \Omega_3^T$</figcaption></figure>

Figure 2: Example 1, pressure computed using bilinear mortars shown on the space-time grid at refinement 3.

<figure><img src="image_6.png" /><figcaption>(a) On $\Omega_1^T \cup \Omega_4^T$</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(b) On $\Omega_2^T \cup \Omega_3^T$</figcaption></figure>

Figure 3: Example 1, *x*-component of velocity computed using bilinear mortars shown on the space-time grid at refinement 3.
<figure><img src="image_8.png" /><figcaption>Figure 4: Example 1, y-component of velocity computed using bilinear mortars shown on the space-time grid at refinement 3.</figcaption></figure>

mortar case $DGQ_2$ is compensated by the higher degree of the space. In particular, the last term on the right hand side in the bound (5.4) from Theorem 5.1 gives $O(h^{-\frac{1}{2}}H^{m+1})$. With $m=1$ and $H = O(h)$, this results in $O(h^{\frac{3}{2}})$, while with $m=2$ and $H = O(h^{\frac{1}{2}})$, it gives $O(h)$. Similarly, the last term in the bound (6.16) from Theorem 6.3 gives $O(H^{m+\frac{1}{2}})$. With $m=1$ and $H = O(h)$, this results in $O(h^{\frac{3}{2}})$, while with $m=2$ and $H = O(h^{\frac{1}{2}})$, it gives $O(h^{\frac{5}{4}})$. In all cases, the order is no smaller than the optimal order $O(h)$ with respect to the $RT_0 \times DGQ_0$ finite element spaces. More details on the mesh refinement are given in Table 1. There we also report the number of spatial degrees of freedom of the spaces $RT_0 \times DGQ_0$ on $\Omega_i$, as well as the number of space-time degrees of freedom of the mortar space $DGQ_1$ or $DGQ_2$ on $\Gamma^T$.

In Tables 2 and 3 we report the relative errors with respect to the norm of the true solution, as well as the convergence rates as powers of the subdomain discretization parameters $h$ and $\Delta t$. We observe optimal first order convergence of the method using both bilinear and biquadratic mortars. We note that the $\Delta t^{-\frac{1}{2}}$ loss in convergence rate in the bound from Theorem 5.1 is not observed in the numerical results. In Tables 2 and 3 we also report the growth rate for the number of GMRES iterations in the case of bilinear and biquadratic mortars, respectively. We recall that Theorem 7.2 bounds the spectral ratio of the interface operator $S$ by $O(h^{-1})$. Thus, up to deviation from a normal matrix, the growth rate is expected to be $O(h^{-0.5})$ [37,38]. This is close to what we observe in Tables 2 and 3.

We further compare the performance of bilinear and biquadratic mortars. As expected, the accuracy of the two cases at the same refinement level is comparable, which is evident from Tables 2 and 3. On the other hand, since the biquadratic mortar space $DGQ_2$ has far fewer degrees of freedom compared to bilinear mortar space $DGQ_1$, the former results in a smaller number of GMRES iterations. Thus, choosing higher mortar degrees $m, s$ with a coarser mortar mesh results in a computationally more efficient method compared to using smaller $m, s$ and a finer mortar mesh.

Finally, the computed solution is presented in three-dimensional space-time plots in Figures 2–4. Note that the $z$-axis corresponds to time direction $t$. The plots clearly show the continuity of pressure and velocity, which is imposed in a weak sense, across the space-time interfaces.
## 8.2 Example 2: problem with a boundary layer

In this example, we demonstrate the advantages of applying the multiscale mortar space-time domain decomposition method to a problem where the solution variables, pressure and velocity, vary on different scales across the space-time domain. For this, we use the known solution, $p(x, y, t) = 1000xyt e^{-10(x^2+y^2+\frac{1}{4}t^2)}$ along with permeability $K = I_{2 \times 2}$ to determine the right-hand side $q$ in (2.1) and impose Dirichlet boundary condition and initial condition on the space-time domain $\Omega^T = (0, 1)^2 \times (0, 0.5)$. By construction, $p(x, y, t)$ varies rapidly in both space and time along the lower-left corner of $\Omega^T$, forming a sharp boundary layer. The pressure decays exponentially away from this corner and is close to zero over a large part of the domain. This calls for an efficient multiscale method that would take advantage of the multiscale nature of the problem and provide better resolution around the lower-left corner compared to the rest of $\Omega^T$.

We partition $\Omega^T$ into $4 \times 4$ identical space-time subdomains $\Omega_i^T$. From the knowledge about the variation of the true pressure, we design a multiscale space-time grid on $\Omega^T$, where the refinement of the grid on each $\Omega_i^T$ is proportional to the amount of pressure variation. The finest mesh on $\Omega_i^T$ has $h = 1/128$ and $\Delta t = 1/64$, and the coarsest mesh on $\Omega_i^T$ has $h = 1/8$ and $\Delta t = 1/8$, see Figures 5–6 for the mesh refinement. The coarser meshes on the majority of the space-time subdomains reduce the computational complexity of the subdomain solves associated with them. We use a linear mortar ($m = s = 1$) on the subdomain interfaces. The mortar mesh sizes in space are chosen as follows. For vertical interfaces (fixed $x$) between subdomains on the bottom row, the one along the boundary layer, we set $H = 1/32$. For the next row of subdomains we set $H = 1/16$, and for the other two rows, $H = 1/8$. Similarly, for the horizontal interfaces (fixed $y$) between subdomains on the left column we set $H = 1/32$, for the second column, $H = 1/16$, and for the other two columns, $H = 1/8$. We choose $\Delta T = 1/8$ on all interfaces. These choices guarantee that the mortar assumption (4.6) is satisfied and that the dimension of the interface problem is reduced, while at the same time provide suitable resolution to enforce weakly flux continuity across the space-time interfaces. For comparison, we solve the problem using a uniformly fine and matching subdomain mesh with $h = H = 1/128$ and $\Delta t = \Delta T = 1/64$. A comparison of the number of GMRES iterations and the relative errors from the multiscale and the fine-scale methods is given in Table 4. It shows that the multiscale and the fine-scale solutions attain comparable accuracy. We observe slightly smaller relative errors in the fine-scale case because of the matching grids and higher resolution throughout the space-time domain. The slightly higher errors for the multiscale method are compensated by the cheaper subdomain solves and smaller interface problem that converges in fewer iterations compared to the fine-scale method.

<table><caption>Table 4: Example 2, errors and GMRES iterations for the multiscale and fine-scale methods</caption><thead><tr><th>Method</th><th># GMRES</th><th>||u - u<sub>h</sub><sup>Δt</sup>||<sub>L<sup>2</sup>(0,T;L<sup>2</sup>(Ω))</sub></th><th>||p - p<sub>h</sub><sup>Δt</sup>||<sub>DG</sub></th><th>||p - p<sub>h</sub><sup>Δt</sup>||<sub>L<sup>2</sup>(0,T;W)</sub></th><th>||λ - λ<sub>H</sub><sup>ΔT</sup>||<sub>L<sup>2</sup>(0,T;λ<sub>H</sub>)</sub></th></tr></thead><tbody><tr><td>multiscale</td><td>102</td><td>5.657e-02</td><td>8.425e-02</td><td>6.319e-02</td><td>5.796e-02</td></tr><tr><td>fine-scale</td><td>140</td><td>1.524e-02</td><td>2.234e-02</td><td>2.154e-02</td><td>3.016e-02</td></tr></tbody></table>

The computed multiscale solution is presented in Figures 5–8. The plots show that the multiscale method provides good resolution where it matters – in the regions with high solution variation. Moreover, we observe very good enforcement of continuity of both pressure and velocity across various space and time interfaces. Side-to-side comparisons of the multiscale and fine-scale solutions are given on the right sides in Figure 6 and Figure 8. They show excellent agreement between the
<figure><img src="image_9.png" /><figcaption>Figure 5: Example 2, pressure from the multiscale method, cut along the plane $x = 0.25$.</figcaption></figure>

two solutions and once again confirm that the less expensive multiscale method provides comparable accuracy to the more expensive fine-scale method.

## 9 Conclusions

We presented a space-time domain decomposition mixed finite element method for parabolic problems that allows for non-matching spatial grids and local time stepping via space-time mortar finite elements. Well-posedness and a priori error estimates were established. A parallel non-overlapping domain decomposition algorithm was developed, which reduces the algebraic problem to a coarse-scale interface problem for the mortar variable. The theoretical results and the flexibility of the method were illustrated in a series of numerical experiments. Future work may include the development and analysis of a Neumann–Neumann preconditioner for the interface problem in Section 7 as in [33], using techniques from [46], as well as deriving a posteriori error estimates, possibly building upon the ideas from [1, 2, 18].

## References

[1] S. Ali Hassan, C. Japhet, M. Kern, and M. Vohralík. A posteriori stopping criteria for optimized Schwarz domain decomposition algorithms in mixed formulations. *Comput. Methods Appl. Math.*, 18(3):495–519, 2018.

[2] S. Ali Hassan, C. Japhet, and M. Vohralík. A posteriori stopping criteria for space-time domain decomposition for the heat equation in mixed formulations. *Electron. Trans. Numer. Anal.*, 49:151–181, 2018.

[3] T. Arbogast, L. C. Cowsar, M. F. Wheeler, and I. Yotov. Mixed finite element methods on nonmatching multiblock grids. *SIAM J. Numer. Anal.*, 37(4):1295–1315, 2000.

[4] T. Arbogast, D. Estep, B. Sheehan, and S. Tavener. A posteriori error estimates for mixed finite element and finite volume methods for parabolic problems coupled through a boundary. *SIAM/ASA J. Uncertain. Quantif.*, 3(1):169–198, 2015.
<figure><img src="image_10.png" /><figcaption>Figure 6: Example 2, left: pressure from the multiscale method, cut along the plane $t = 0.35$; right: pressure from the multiscale (top) and fine-scale (bottom) methods on the whole domain.</figcaption></figure>

[5] T. Arbogast, G. Pencheva, M. F. Wheeler, and I. Yotov. A multiscale mortar mixed finite element method. *Multiscale Model. Simul.*, 6(1):319–346, 2007.

[6] M. Arshad, E.-J. Park, and D. Shin. Multiscale mortar mixed domain decomposition approximations of nonlinear parabolic equations. *Comput. Math. Appl.*, 97:375–385, 2021.

[7] W. Bangerth, R. Hartmann, and G. Kanschat. deal.II—a general-purpose object-oriented finite element library. *ACM Trans. Math. Software*, 33(4):Art. 24, 27, 2007.

[8] M. Bause, F. A. Radu, and U. Köcher. Error analysis for discretizations of parabolic problems using continuous finite elements in time and mixed finite elements in space. *Numer. Math.*, 137(4):773–818, 2017.

[9] B. Beckermann, S. A. Goreinov, and E. E. Tyrtyshnikov. Some remarks on the Elman estimate for GMRES. *SIAM J. Matrix Anal. Appl.*, 27(3):772–778, 2005.

[10] 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.

[11] F. Brezzi and M. Fortin. *Mixed and hybrid finite element methods*. Springer-Verlag, New York, 1991.
<figure><img src="image_11.png" /><figcaption>Figure 7: Example 2, velocity magnitude from the multiscale method, cut along the plane $x = 0.25$.</figcaption></figure>

[12] J. M. Cascón, L. Ferragut, and M. I. Asensio. Space-time adaptive algorithm for the mixed parabolic problem. *Numer. Math.*, 103(3):367–392, 2006.

[13] P. G. Ciarlet. *The finite element method for elliptic problems*. North-Holland Publishing Co., Amsterdam-New York-Oxford, 1978. Studies in Mathematics and its Applications, Vol. 4.

[14] M. Crouzeix and V. Thomée. The stability in $L_p$ and $W_p^1$ of the $L_2$-projection onto finite element function spaces. *Math. Comp.*, 48(178):521–532, 1987.

[15] C. N. Dawson, Q. Du, and T. F. Dupont. A finite difference domain decomposition algorithm for numerical solution of the heat equation. *Math. Comp.*, 57(195):63–71, 1991.

[16] L. Delpopolo Carciopolo, M. Cusini, L. Formaggia, and H. Hajibeygi. Adaptive multilevel space-time-stepping scheme for transport in heterogeneous porous media (ADM-LTS). *J. Comput. Phys. X*, 6:100052, 21, 2020.

[17] S. C. Eisenstat, H. C. Elman, and M. H. Schultz. Variational iterative methods for nonsymmetric systems of linear equations. *SIAM J. Numer. Anal.*, 20(2):345–357, 1983.

[18] A. Ern, I. Smears, and M. Vohralík. Guaranteed, locally space-time efficient, and polynomial-degree robust a posteriori error estimates for high-order discretizations of parabolic problems. *SIAM J. Numer. Anal.*, 55(6):2811–2834, 2017.

[19] R. E. Ewing, R. D. Lazarov, and P. S. Vassilevski. Finite difference schemes on grids with local refinement in time and space for parabolic problems. I. Derivation, stability, and error analysis. *Computing*, 45(3):193–215, 1990.

[20] R. D. Falgout, S. Friedhoff, T. V. Kolev, S. P. MacLachlan, and J. B. Schroder. Parallel time integration with multigrid. *SIAM J. Sci. Comput.*, 36(6):C635–C661, 2014.

[21] V. Faucher and A. Combescure. A time and space mortar method for coupling linear modal subdomains and non-linear subdomains in explicit structural dynamics. *Comput. Methods Appl. Mech. Engrg.*, 192(5):509–533, 2003.

[22] S. Gaiffe, R. Glowinski, and R. Masson. Domain decomposition and splitting methods for mortar mixed finite element approximations to parabolic equations. *Numer. Math.*, 93(1):53–75, 2002.
<figure><img src="image_12.png" /><figcaption>Figure 8: Example 2, left: velocity magnitude from the multiscale method, cut along the plane $t = 0.35$; right: velocity magnitude from the multiscale (top) and fine-scale (bottom) methods on the whole domain.</figcaption></figure>

[23] M. J. Gander. 50 years of time parallel time integration. In *Multiple shooting and time domain decomposition methods. MuS-TDD, Heidelberg, Germany, May 6-8, 2013*, pages 69-113. Cham: Springer, 2015.

[24] M. J. Gander and L. Halpern. Optimized Schwarz waveform relaxation methods for advection-reaction diffusion problems. *SIAM J. Numer. Anal.*, 45(2):666-697, 2007.

[25] M. J. Gander, F. Kwok, and B. C. Mandal. Dirichlet-Neumann and Neumann-Neumann waveform relaxation algorithms for parabolic problems. *Electron. Trans. Numer. Anal.*, 45:424-456, 2016.

[26] M. J. Gander and M. Neumüller. Analysis of a new space-time parallel multigrid algorithm for parabolic problems. *SIAM J. Sci. Comput.*, 38(4):A2173-A2208, 2016.

[27] M. J. Gander and S. Vandewalle. Analysis of the parareal time-parallel time-integration method. *SIAM J. Sci. Comput.*, 29(2):556-578, 2007.

[28] B. Ganis and I. Yotov. Implementation of a mortar mixed finite element method using a multiscale flux basis. *Comput. Methods Appl. Mech. Engrg.*, 198(49):3989-3998, 2009.
[29] A. Greenbaum. *Iterative methods for solving linear systems*, volume 17 of *Frontiers in Applied Mathematics*. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1997.

[30] P. Grisvard. *Elliptic problems in nonsmooth domains*, volume 69 of *Classics in Applied Mathematics*. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2011.

[31] C. Hager, P. Hauret, P. Le Tallec, and B. I. Wohlmuth. Solving dynamic contact problems with local refinement in space and time. *Comput. Methods Appl. Mech. Engrg.*, 201/204:25-41, 2012.

[32] L. Halpern, C. Japhet, and J. Szeftel. Optimized Schwarz waveform relaxation and discontinuous Galerkin time stepping for heterogeneous problems. *SIAM J. Numer. Anal.*, 50(5):2588-2611, 2012.

[33] T.-T.-P. Hoang, J. Jaffré, C. Japhet, M. Kern, and J. E. Roberts. Space-time domain decomposition methods for diffusion problems in mixed formulations. *SIAM J. Numer. Anal.*, 51(6):3532-3559, 2013.

[34] T.-T.-P. Hoang, C. Japhet, M. Kern, and J. E. Roberts. Space-time domain decomposition for reduced fracture models in mixed formulation. *SIAM J. Numer. Anal.*, 54(1):288-316, 2016.

[35] T.-T.-P. Hoang and H. Lee. A global-in-time domain decomposition method for the coupled nonlinear Stokes and Darcy flows. *J. Sci. Comp.*, 87(1):1-22, 2021.

[36] R. A. Horn and C. R. Johnson. *Topics in matrix analysis*. Cambridge University Press, Cambridge, 1994. Corrected reprint of the 1991 original.

[37] I. C. F. Ipsen. Expressions and bounds for the GMRES residual. *BIT*, 40(3):524-535, 2000.

[38] C. T. Kelley. *Iterative methods for linear and nonlinear equations*, volume 16 of *Frontiers in Applied Mathematics*. Society for Industrial and Applied Mathematics, Philadelphia, 1995.

[39] W. Kheriji, R. Masson, and A. Moncorgé. Nearwell local space and time refinement in reservoir simulation. *Math. Comput. Simulation*, 118:273-292, 2015.

[40] D. Kim, E.-J. Park, and B. Seo. Space-time adaptive methods for the mixed formulation of a linear parabolic Problem. *J. Sci. Comput.*, 74(3):1725-1756, 2018.

[41] D. Krause and R. Krause. Enabling local time stepping in the parallel implicit solution of reaction-diffusion equations via space-time finite elements on shallow tree meshes. *Appl. Math. Comput.*, 277:164-179, 2016.

[42] J.-L. Lions, Y. Maday, and G. Turinici. Résolution d'EDP par un schéma en temps “pararéel”. *C. R. Acad. Sci., Paris, Sér. I, Math.*, 332(7):661-668, 2001.

[43] J.-L. Lions and E. Magenes. *Non-homogeneous boundary value problems and applications. Vol. I*. Springer-Verlag, New York-Heidelberg, 1972.

[44] C. Makridakis and R. H. Nochetto. A posteriori error analysis for higher order dissipative methods for evolution problems. *Numer. Math.*, 104:489-514, 2006.

[45] P. B. Nakshatrala, K. B. Nakshatrala, and D. A. Tortorelli. A time-staggered partitioned coupling algorithm for transient heat conduction. *Internat. J. Numer. Methods Engrg.*, 78(12):1387-1406, 2009.
[46] G. Pencheva and I. Yotov. Balancing domain decomposition for mortar mixed finite element methods. *Numer. Linear Algebra Appl.*, 10(1-2):159–180, 2003.

[47] I. Rybak and J. Magiera. A multiple-time-step technique for coupled free flow and porous medium systems. *J. Comput. Phys.*, 272:327–342, 2014.

[48] L. R. Scott and S. Zhang. Finite element interpolation of nonsmooth functions satisfying boundary conditions. *Math. Comput.*, 54(190):483–493, 1990.

[49] G. Starke. Field-of-values analysis of preconditioned iterative methods for nonsymmetric el-
liptic problems. *Numer. Math.*, 78(1):103–117, 1997.

[50] V. Thomée. Galerkin finite element methods for parabolic problems, volume 25 of Springer Series in Computational Mathematics. Springer-Verlag, Berlin, second edition, 2006.

[51] H. Yu. A local space-time adaptive scheme in solving two-dimensional parabolic problems based on domain decomposition methods. *SIAM J. Sci. Comput.*, 23(1):304–322, 2001.