# Numerical approximation of the unique continuation problem enriched by a database for the Stokes equations

Muriel Boulakia, Corrie James, Damiano Lombardi

► To cite this version:

Muriel Boulakia, Corrie James, Damiano Lombardi. Numerical approximation of the unique continuation problem enriched by a database for the Stokes equations. ESAIM: Mathematical Modelling and Numerical Analysis, 2025, 59 (3), pp.1399-1435. 10.1051/m2an/2025024. hal-04721560v3

HAL Id: hal-04721560

https://inria.hal.science/hal-04721560v3

Submitted on 29 Aug 2025

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
NUMERICAL APPROXIMATION OF THE UNIQUE CONTINUATION PROBLEM
ENRICHED BY A DATABASE FOR THE STOKES EQUATIONS

MURIEL BOULAKIA<fnref n="1" />, CORRIE JAMES <fnref n="1" />,<fnref n="2" />,*, AND DAMIANO LOMBARDI <fnref n="2" />

ABSTRACT. This paper studies the unique continuation problem for the Stokes equations given a database of population measurements. The problem is set up as a minimization problem under a PDE constraint and discretized using the finite element method. It is then regularized by the population data, by imposing that the solution lives near a finite-dimensional subspace generated by the database. This study examines how the inclusion of population data in the resolution improves both theoretical and numerical results. Using the proposed method, global error estimates for the velocity and pressure are obtained at an improved rate of convergence. The inclusion of the population data also has a positive impact on the numerical test cases, in both 2D and 3D, and especially when the measurements are scarce.

***Keywords:*** data assimilation, unique continuation, Stokes

# 1. INTRODUCTION

When working with experimental data, researchers are often limited by the machines and techniques used to collect measurements. For example, measurements might be scarce or available only in a restricted subdomain. Being able to estimate what happens outside this region of observation is thus of interest in many fields. In a unique continuation formulation, the measurements are seen as observations of a solution to a Partial Differential Equation (PDE) in a subdomain. The idea is then to reconstruct the solution to the PDE on a larger domain. As, in general, we do not know certain information about the PDE, like parameters or boundary conditions, the measurements are combined with the PDE model in order to remove its underdetermined nature.

A realistic application of this unique continuation problem is reconstructing the velocity and pressure of blood flow in an artery given images of the flow inside an arterial segment. This problem can be modeled by PDE-based flow models with incomplete boundary data on the inlet and outlet, which correspond to the artificial boundaries obtained when the artery domain is truncated. In addition to finding an extension of the blood velocity outside the zone of observation, solving the unique continuation problem also provides an estimate of the vascular pressure. This estimation is of interest because the pressure cannot be directly measured without invasive procedures and can be an important indicator in assessing the severity of a pathology.

Our paper deals with a model problem related to this application: given a fluid flow described by the stationary Stokes equations in a bounded domain with missing boundary conditions and velocity measurements taken in a subdomain, we numerically reconstruct the velocity and pressure in the whole domain. In contrast with previous studies addressing this question, in this work, in addition to a single individual's measurements, we assume access to a database of measurements taken on a set of experiments or a population of individuals. In the present work, we will refer to these measurements as the database or the population data. This corresponds to a realistic framework since, in general, measuring devices are used repeatedly under similar experimental conditions. Thus, we will consider population data which are

<footnote n="1">¹UNIVERSITÉ PARIS-SACLAY, UVSQ, CNRS, LABORATOIRE DE MATHÉMATIQUES DE VERSAILLES, 78000 VERSAILLES, FRANCE.</footnote>

<footnote n="2">²SORBONNE UNIVERSITÉ, INRIA AND CNRS, UMR 7598 LABORATOIRE JACQUES-LOUIS LIONS, PARIS, FRANCE.</footnote>

<footnote>*CORRESPONDING AUTHOR: CORRIE.JAMES@INRIA.FR</footnote>

<footnote>Date: Published 27 May 2025.</footnote>
of the same nature as the individual's measurements: same level of noise and same number and location of points of measurement. In this framework, our main objective is to investigate to what extent taking into account the population information improves theoretical and numerical results.

The study of numerical methods for the unique continuation problem has received a lot of attention in the case of single measurements. A common approach is to view it as a PDE-constrained optimization problem, in which the distance between the velocity solution and the measurements in the region of observation is minimized under the constraint of the Stokes equations. Since the reconstruction problem is ill-posed, a regularization step is usually added to this formulation. As presented in the course notes [13], there are basically two main regularization families: first, a regularization can be applied at the continuous level before discretization, classically with Tikhonov or quasi-reversibility methods. In this case, the main drawback is that the presence of the regularization affects the accuracy of the reconstruction since it amounts to adding an a priori on the reconstructed solution. Recently, a discretize first, then regularize alternative has been proposed and analyzed in several papers like [7, 10, 12, 13, 14] since the seminal work [9]. In these papers, the design of the regularization terms is based on stabilization methods usually applied to the finite element (FE) approximation of well-posed problems, such as in [16] and [25]. In what follows, we will refer to the “discretize then regularize” method as the *stabilized FE method for inverse problems*. In this case, contrary to Tikhonov-type regularization, the reconstruction is unaffected by the regularization terms since they disappear as the mesh size tends to 0 allowing the numerical reconstructed solution to converge to the exact solution. In the papers cited above, this consistency property is highlighted by numerical reconstruction error estimates dependent on the mesh size and level of noise. In all these studies, these error estimates are fundamentally based on the stability estimates obtained for the unique continuation problem at the continuous level and consequently are of the same nature (for instance, global logarithmic stability estimates at the continuous level lead to global error estimates of the form $|\log(h)|^{-\alpha}$ for $\alpha > 0$).

More specifically, in the case of the Stokes equations, stability estimates quantifying the unique continuation property established in [19] have been studied in papers [22, 4, 8, 3]. These estimates can be classified into two categories of results: Hölderian inequalities for local estimates as in [22, 8] and logarithmic inequalities for global estimates, see [4, 8, 3]. The stabilized FE method for inverse problems was studied in [11] in the nonconforming case applied to the Stokes equations and in [7] in the conforming case applied to the linearized Navier-Stokes equations. In both papers, an error bound of the form $h^\tau$ for $\tau \in (0,1)$, with $h$ the mesh size, is proved. In [10], the authors generalize this result to higher order polynomial approximation, finding an error bound in $h^{k\tau}$ where $k$ is the polynomial order of the FE method. Both of these papers rely on a stability estimate from [22], which is a local estimate on the velocity (not the pressure). As a consequence, the error estimates in [7, 10] act locally on only the velocity. According to theoretical results proven in [8, 3], global error estimates both on the velocity and pressure are achievable, but the estimate would be only logarithmic, which for practical cases is inadequate. Moreover, local pressure measurements would be required, in addition to the velocity measurements, which is not realistic.

As pointed out in the recent paper [14], it is possible to improve the weak error estimates satisfied by the solution of the numerical approximation of the inverse problem if additional a priori information is available and incorporated into the optimization method. In [14], the authors study the unique continuation problem for the Laplacian under the additional hypothesis that the trace of the solution is finite. This specific hypothesis allows the estimates of the global reconstruction error to be improved from logarithmic, as in the general case, to Lipschitz. In the same vein, our work studies to what extent the numerical reconstruction is improved when a database of measurements is used, in addition to a single individual's measurements. More specifically, we propose a new numerical method for the unique continuation problem applied to the Stokes equation in the case when a database is available.

In our work, the database of measurements is incorporated into the minimization problem by adding a regularizing term that enforces the solution to be close to this database. Concretely, we first construct a finite-dimensional subspace generated by the population data through the Proper Orthogonal Decomposition (POD). The POD basis, which is defined on the subdomain where the population data
are measured, is then extended to the entire domain through the resolution of an optimal recovery problem, taking advantage of the POD's denoising properties. The finite-dimensional subspace encoding the population information is then the span of the extended POD modes, which is included in the PDE-constrained optimization problem through a term imposing that the velocity and pressure solutions be near, but not necessarily inside, this POD-generated subspace. The inclusion of population data in the resolution is our first main contribution in this work. As we will see, the inclusion of this term has a regularizing effect on the problem, rendering stabilization terms for the primal variable unnecessary. However, in numerical cases, we observe that the presence of these stabilization terms is sometimes beneficial, so we have included them in our study.

This new method is studied from both theoretical and numerical points of view in order to assess the improvements brought by the population information. In our theoretical study, we show global error estimates both on the velocity and pressure with an error bound composed of two terms: a term in $h^k$ where $k$ is the polynomial order of the FE method and a remainder term corresponding to the distance between the exact solution and the vector space spanned by the extended POD modes. This is our second main contribution. The only hypotheses required are a standard a priori bound on the solution, which is assumed in all the previously mentioned works, and the existence of a database of measurements, which is specific to our study.

The third main contribution is a thorough numerical study to assess the performance of our population-enriched method when there is a fixed number of scarce data measurements. In this scenario, numerical results improve across the board when the database is included in the resolution, in both 2D and 3D test cases. In particular, we present a critical 3D case where the improvement brought by the population-enriched method is especially noticeable for the reconstruction of both velocity and pressure.

Our paper is organized in four sections. In Section 2, the unique continuation problem is set up, and all the elements of the PDE-constrained optimization problem, such as the regularization by population data term, are introduced. In Section 3, we analyze our method, showing it is well-posed and studying its convergence, obtaining the global error estimates. A perturbation analysis and a study of the non-homogeneous Stokes problem are performed as well. Section 4 introduces the numerical method used to synthetically construct the database, generate the POD basis, obtain the extended POD modes, and solve the unique continuation problem. Finally, in Section 5, details about our specific test cases are given and numerical results are displayed. We consider two data regimes, sparse and non-sparse, study the reconstruction error as a function of the mesh-size, perform a statistical study of the errors, and finish with a critical 3D test case.

## 2. PROBLEM SETTING AND PRESENTATION OF THE NUMERICAL METHOD

Let $\Omega$ be a bounded, Lipschitzian domain in $\mathbb{R}^d$ where $d = 2,3$ and let $\omega \subsetneq \Omega$ be the region of observation. The Stokes unique continuation problem can be written as: find $(u,p) \in H^1(\Omega)^d \times L_0^2(\Omega)$ such that

$$ \left\{ \begin{array}{ll} -\mu \Delta u + \nabla p = 0 & \text{in } \Omega, \\ \nabla \cdot u = 0 & \text{in } \Omega, \\ u|_{\omega} = u_M & \text{in } \omega \end{array} \right. \qquad (1) $$

where $\mu > 0$ is the viscosity and $u_M$ corresponds to the velocity measured in the region of observation. The space $L_0^2(\Omega)$ designates the space of $L^2$ functions of zero mean in $\Omega$.

In this work, we place ourselves in the framework where a database of $L^2$ measurements in $\omega$, noted $\{u_M^{(j)}\}_{1 \le j \le N_{pop}}$, is available, with $N_{pop} > 0$ the size of the database.

We make the hypothesis that the measurements $u_M$ come from a family of unknown solutions to the Stokes problem. More precisely, we assume that, for all $1 \le j \le N_{pop}$ there exists:

$$ (\bar{u}^{(j)}, \bar{p}^{(j)}) \in H^1(\Omega)^d \times L_0^2(\Omega) \text{ such that } u_M^{(j)} = \bar{u}^{(j)}|_{\omega} \text{ and } (\bar{u}^{(j)}, \bar{p}^{(j)}) \text{ solution to the Stokes equations.} \qquad (2) $$
Additionally, we assume that the {$u_M^{(j)}$}$_{1\le j\le N_{pop}}$ come from the same population, of which we have some samples. By “same population”, we mean that the measurements are observations of solutions to parametric PDEs where the parameters are sampled from the same distribution. In our numerical tests, the population database is synthetically generated by parametrizing the boundary conditions with parameters drawn from uniform distributions and solving forward problems, see Section 5. We assume that the measurements $u_M$ and the population database {$u_M^{(j)}$}$_{1\le j\le N_{pop}}$ are contained in a bounded subset of $L^2(\omega)^d$. We also assume that the measurements $u_M$ are close to a subspace generated by the measurements of the population database. Moreover, we assume that we have a global a priori bound: there exists $C_{pop} > 0$ such that for all $1 \le j \le N_{pop}$,

$$ \|\bar{u}^{(j)}\|_{H^1(\Omega)^d} + \|\bar{p}^{(j)}\|_{L^2(\Omega)} \le C_{pop}. \quad (3) $$

In a classical way, problem (1) is reformulated as an optimization problem, where the data misfit functional is minimized under the PDE constraint of satisfying the Stokes problem. The discrete Lagrangian associated with this optimization problem is then regularized by incorporating the information coming from the population data, in addition to the stabilization terms which are usually present in the stabilized FE method for inverse problems and which are optional in our context.

In what follows, we present the different building blocks that make up the discrete Lagrangian associated with our problem. First, the formulation as a constrained optimization problem is presented in Subsection 2.1. In Subsection 2.2, the method used to include the population data in the problem is described. Through the use of a POD basis, the population data are approximated by a finite dimensional subspace and, in the Lagrangian, we will impose that the solution $(u, p)$ lives near, but not necessarily inside, this space. Next, some stabilization terms are presented in Subsection 2.3. We observed that the stabilization terms on the primal variables, although not necessary in the theoretical study, improve some of the numerical results in a non-negligible way, so they are included. Finally, Subsection 2.4 summarizes all these elements by introducing the Lagrangian associated with our discrete reconstruction.

## 2.1. Variational Formulation and Finite Element Approximation

We first introduce the standard variational formulation associated with the Stokes equations: find $(u, p) \in H^1(\Omega)^d \times L_0^2(\Omega)$ such that

$$ A[(u, p), (v, q)] = 0 $$

for all $(v, q) \in H_0^1(\Omega)^d \times L^2(\Omega)$ where $A$ is the operator defined in $(H^1(\Omega)^d \times L_0^2(\Omega)) \times (H_0^1(\Omega)^d \times L^2(\Omega))$ by

$$ A[(u, p), (v, q)] = \mu \int_{\Omega} \nabla u : \nabla v - \int_{\Omega} p \nabla \cdot v + \int_{\Omega} q \nabla \cdot u. $$

On the domain $\Omega$, we consider a family of conforming meshes ($\mathcal{T}_h$) consisting of shape regular simplices $K$. We denote $h$ the mesh size and $k \ge 1$ the polynomial order of the finite element space and introduce the following standard $\mathbb{P}_k$ Finite Element spaces:

$$ X_h^k = \{v_h \in H^1(\Omega) : v_{h|_K} \in \mathbb{P}_k, \forall K \in \mathcal{T}_h\} $$

$$ V_h^k = [X_h^k]^d, \quad W_h^k = V_h^k \cap H_0^1(\Omega)^d, \quad Q_h^{0,k} = X_h^k \cap L_0^2(\Omega). $$

The finite element approximation then reads: find $(u_h, p_h) \in V_h^k \times Q_h^{0,k}$ such that for all $(v_h, q_h) \in W_h^k \times X_h^k$

$$ A[(u_h, p_h), (v_h, q_h)] = 0. $$

We notice that in the absence of boundary conditions, the spaces of the solutions and of the test functions do not coincide. This is related to the under-determined nature of the problem in the absence of measurements. This weak formulation of the Stokes equations will be included in the Lagrangian below, see (9), as a constraint on the problem of minimizing the data misfit functional

$$ \frac{\gamma_M}{2} \|u_h - u_M\|_{L^2(\omega)^d}^2 $$

where $\gamma_M > 0$ is a constant.
**Remark 2.1.** Should boundary conditions be partially known, the variational formulation would be slightly different. An example of this is found in the numerical sections 4 and 5. Denoting $\Gamma_D$ the boundary where there are known Dirichlet conditions, $\Gamma_N$ the boundary where there are known Neumann conditions, and $\Gamma_U$ the boundary where the conditions are unknown, the test function $v$ would be taken in the set of $H^1(\Omega)^d$ functions such that $w = 0$ on $\Gamma_D \cup \Gamma_U$, and the variational formulation would take into account the known boundary conditions.

## 2.2. Regularization from the Measurement Database: Use of the POD Method

Throughout the paper, we will denote the $L^2$ norm and $L^2$ scalar product over some domain $D \subset \mathbb{R}^d$ respectively by $\|\cdot\|_D$ and $\langle\cdot, \cdot\rangle_D$.

In this subsection, we describe the design of the population data-based regularization terms that will appear in our Lagrangian formulation. To do so, we first introduce the POD basis associated with the population data $\{u_M^{(j)}\}_{1\le j\le N_{pop}}$. POD is a widely applied reduced order method which combines Singular Value Decomposition with a reduced basis approach. The appeal of this method is that the POD basis generates a vector space optimally approximating the population data in the least squares sense. The construction of the POD basis presented in the following paragraphs is based on the work of Kunisch and Volkwein in [21].

Let us first define

$$ \mathcal{Y}_{N_{pop}} : \mathbb{R}^{N_{pop}} \to L^2(\omega)^d \\ (v_j)_{1\le j\le N_{pop}} \longmapsto \sum_{j=1}^{N_{pop}} v_j u_M^{(j)}. $$

Its adjoint is given by:

$$ \mathcal{Y}_{N_{pop}}^* : L^2(\omega)^d \to \mathbb{R}^{N_{pop}} \\ z \mapsto (\langle z, u_M^{(1)} \rangle_{\omega}, \dots, \langle z, u_M^{(N_{pop})} \rangle_{\omega})^T. $$

Let us pose

$$ \mathcal{R}_{N_{pop}} = \mathcal{Y}_{N_{pop}} \mathcal{Y}_{N_{pop}}^* \in \mathcal{L}(L^2(\omega)^d) $$

and

$$ \mathcal{K}_{N_{pop}} = \mathcal{Y}_{N_{pop}}^* \mathcal{Y}_{N_{pop}} \in \mathbb{R}^{N_{pop} \times N_{pop}}. $$

We denote by $k^* \le N_{pop}$ the rank of $\mathcal{Y}_{N_{pop}}$. We see that $\mathcal{R}_{N_{pop}}$ is bounded, auto-adjoint, positive, with an image of finite dimension $k^*$, and therefore compact. By the Hilbert-Schmidt theorem, there exists a Hilbert basis $(\varphi_i)_{i \in \mathbb{N}} \in L^2(\omega)^d$ and a sequence of non-negative real numbers $(\lambda_i)_{i \in \mathbb{N}}$ such that

$$ \lambda_1 \ge \dots \ge \lambda_{k^*} > 0, \quad \lambda_i = 0 \text{ for } i > k^* $$

and

$$ \mathcal{R}_{N_{pop}} \varphi_i = \lambda_i \varphi_i \quad \forall i \in \mathbb{N}^*. $$

For $1 \le i \le k^*$, we define $v_i = \frac{1}{\sqrt{\lambda_i}} \mathcal{Y}_{N_{pop}}^* \varphi_i$. We remark that

$$ \mathcal{K}_{N_{pop}} v_i = \lambda_i v_i, \quad 1 \le i \le k^*. $$

The family $(v_i)_{1\le i\le k^*}$ is a basis of eigenfunctions of $\mathcal{K}_{N_{pop}}$ for its image, and a simple calculation shows that it is orthonormal. Moreover, we notice that the $k^*$ eigenfunctions of $\mathcal{R}_{N_{pop}}$ associated with the nonzero eigenvalues can be expressed as:

$$ \varphi_i = \frac{1}{\sqrt{\lambda_i}} \mathcal{Y}_{N_{pop}} v_i \quad 1 \le i \le k^*, $$

which can also be rewritten

$$ \varphi_i = \frac{1}{\sqrt{\lambda_i}} \sum_{j=1}^{N_{pop}} v_{i,j} u_M^{(j)}. \quad (4) $$
The order reduction is carried out by disregarding the POD basis functions $\varphi_i$ whose associated singular value $\sqrt{\lambda_i}$ is smaller than a cut-off threshold. The data space is therefore approximated by the vector space spanned by the first $n > 1$ basis functions $(\varphi_i)_{1\le i\le n}$ which contain the most relevant information in the least square sense (which can be interpreted as an optimal energy criteria).

The population information, as described by the POD basis $(\varphi_i)_{1\le i\le n}$, is still constrained to the region of observation $\omega$. In our method, the information provided by this basis will be incorporated, not as local information on the measurement domain $\omega$, but on the global domain $\Omega$. To do so, the following proposition introduces an extension of the POD basis functions to the entire domain $\Omega$ and an upper bound.

**Proposition 2.2.** Let $1 \le i \le n$. For each POD basis function $\varphi_i$, there exists a unique, extended POD mode pair $(\xi_i, \xi_i^p) \in H^1(\Omega)^d \times L_0^2(\Omega)$ verifying the following unique continuation problem:

$$ \left\{ \begin{array}{ll} -\mu \Delta \xi_i + \nabla \xi_i^p & = 0 & \text{in } \Omega, \\ \nabla \cdot \xi_i & = 0 & \text{in } \Omega \\ \xi_i & = \varphi_i & \text{in } \omega. \end{array} \right. \qquad (5) $$

The pair also verifies,

$$ \sum_{i=1}^{n} \|\xi_i\|_{H^1(\Omega)^d} + \|\xi_i^p\|_{L^2(\Omega)} \le \hat{C}, \qquad (6) $$

where $\hat{C}$ depends on $N_{pop}$, $\lambda_i$ for $1 \le i \le n$, and $C_{pop}$ which corresponds to the a priori bound in (3).

*Proof.* The uniqueness of the pair $(\xi_i, \xi_i^p)$ is a direct consequence of the unique continuation result for the Stokes equations, proven in [19] in a more general setting. We remark that, since this is the stationary Stokes problem, the uniqueness also follows from the unique continuation result for the Laplacian equation. Furthermore, using hypothesis (2) and the linearity of Stokes equations, we deduce from formula (4) that $\xi_i$ and $\xi_i^p$ are given by

$$ \xi_i = \frac{1}{\sqrt{\lambda_i}} \sum_{j=1}^{N_{pop}} v_{i,j} \bar{u}^{(j)} \quad \text{and} \quad \xi_i^p = \frac{1}{\sqrt{\lambda_i}} \sum_{j=1}^{N_{pop}} v_{i,j} \bar{p}^{(j)}. $$

Using the Cauchy-Schwarz inequality and the fact that the $(v_i)_{1\le i\le n}$ are orthonormal in the $\ell^2$ norm, we see that

$$ \begin{aligned} \|\xi_i\|_{H^1(\Omega)^d} &\le \frac{1}{\sqrt{\lambda_i}} \sum_{j=1}^{N_{pop}} |v_{i,j}| \|\bar{u}^{(j)}\|_{H^1(\Omega)^d} \\ &\le \frac{1}{\sqrt{\lambda_i}} \underbrace{\left(\sum_{j=1}^{N_{pop}} |v_{i,j}|^2\right)}_{=1}^{1/2} \left(\sum_{j=1}^{N_{pop}} \|\bar{u}^{(j)}\|_{H^1(\Omega)^d}^2\right)^{1/2} \\ &= \frac{1}{\sqrt{\lambda_i}} \left(\sum_{j=1}^{N_{pop}} \|\bar{u}^{(j)}\|_{H^1(\Omega)^d}^2\right)^{1/2} \end{aligned} $$

and we conclude by using (3). The calculation is similar for the estimate of $\xi_i^p$ in $L^2(\Omega)$. So, we get (6) with $\hat{C} = C_{pop} \sqrt{N_{pop}} \sum_{i=1}^{n} \frac{1}{\sqrt{\lambda_i}}$. $\square$

**Remark 2.3.** Let us remark that the constant $\hat{C}$ depends on $\sum_{i=1}^n \frac{1}{\sqrt{\lambda_i}}$. If we did not truncate the POD modes at rank $n$, then this upper bound could be quite large since the $\lambda_i$ could become quite small
as $i \to N$. However, because we truncate at an $n$ for which the singular value $\sqrt{\lambda_n}$ is greater than a cut-off threshold, we can conclude that the constant $\hat{C}$'s dependence on $\sum_{i=1}^n \frac{1}{\sqrt{\lambda_i}}$ is not concerning.

The population information is now expressed over the entirety of $\Omega$, and the span of these extended POD modes generates a finite-dimensional subspace, to which we will impose our solution to be close. We define the following bilinear forms: for $(u_h, p_h)$, $(v_h, q_h) \in V_h^k \times Q_h^{k,0}$ let

$$r_{\text{POD}}(u_h, v_h) := \gamma_{\text{POD}} \int_{\Omega} (u_h - \sum_{i=1}^{n} \langle u_h, \varphi_i \rangle_{\omega} \xi_i)(v_h - \sum_{i=1}^{n} \langle v_h, \varphi_i \rangle_{\omega} \xi_i) \quad (7)$$

and

$$r_{\text{POD}}^{\mathbf{p}}((u_h, p_h), (v_h, q_h)) := \gamma_{\text{POD}} \int_{\Omega} (p_h - \sum_{i=1}^{n} \langle u_h, \varphi_i \rangle_{\omega} \xi_i^{\mathbf{p}})(q_h - \sum_{i=1}^{n} \langle v_h, \varphi_i \rangle_{\omega} \xi_i^{\mathbf{p}}). \quad (8)$$

We denote $R_{\text{POD}}$ the quadratic form associated with $r_{\text{POD}}$. It describes the distance from the velocity to the finite dimensional space generated by the extended POD modes $(\xi_i)_{1\le i\le n}$. Similarly, the quadratic form associated with $r_{\text{POD}}^{\mathbf{p}}$ is denoted $R_{\text{POD}}^{\mathbf{p}}$. As seen in the definition (9), these two quadratic forms $R_{\text{POD}}$ and $R_{\text{POD}}^{\mathbf{p}}$ appear in the Lagrangian formulation. They act as regularizing terms by penalizing the case where the solution is far from the POD space.

In what follows, the projection of $u_h$ into the POD space and its complement will be respectively denoted as

$$\mathcal{P}(u_h) := \sum_{i=1}^{n} \langle u_h, \varphi_i \rangle_{\omega} \xi_i, \quad \mathcal{Q}(u_h) := (I - \mathcal{P})(u_h).$$

For the pressure, we define

$$\mathcal{P}^{\mathbf{p}}(u_h) := \sum_{i=1}^{n} \langle u_h, \varphi_i \rangle_{\omega} \xi_i^{\mathbf{p}}$$

and

$$\mathcal{Q}^{\mathbf{p}}(u_h, p_h) := p_h - \sum_{i=1}^{n} \langle u_h, \varphi_i \rangle_{\omega} \xi_i^{\mathbf{p}} = p_h - \mathcal{P}^{\mathbf{p}}(u_h).$$

**Remark 2.4.** The exact extended POD modes $(\xi_i, \xi_i^{\mathbf{p}})_{1\le i\le n}$ are used in the theoretical analysis of the unique continuation problem in Sections 2 and 3, whereas their approximation is used for the numerical study in Sections 4 and 5. This approximation is obtained by an optimal recovery method, which is detailed in Section 4.2. A complete theoretical analysis would use the approximate extended POD modes and is a potential avenue for future work.

## 2.3. Stabilization Terms

In our approach, we have two types of stabilization terms: the Continuous Interior Penalization (CIP), which was introduced by Douglas and Dupont in [16], and the Galerkin Least Squares (GLS) method. We define: for $(u_h, p_h)$, $(v_h, q_h) \in V_h^k \times Q_h^{0,k}$,

$$s[(u_h, p_h), (v_h, q_h)] = s_{\text{GLS}}[(u_h, p_h), (v_h, q_h)] + s_{\text{CIP}}[(u_h, v_h)]$$

where $s_{\text{GLS}}$ is the GLS stabilization given by:

$$s_{\text{GLS}}[(u_h, p_h), (v_h, q_h)] = \frac{\gamma_{\text{GLS}}}{\mu} \sum_{K \in \mathcal{T}_h} h^2 \langle -\mu \Delta u_h + \nabla p_h, -\mu \Delta v_h + \nabla q_h \rangle_{L^2(\mathcal{K})}$$

and $s_{\text{CIP}}$ is the CIP stabilization given by:

$$s_{\text{CIP}}[(u_h, v_h)] = \gamma_{\text{CIP}} \mu \sum_{K \in \mathcal{T}_h} h \langle [[\nabla u_h]], [[\nabla v_h]] \rangle_{L^2(\partial K \setminus \partial \Omega)^d}.$$

Here, $[[\cdot]]$ represents a jump across a face $F$ of two adjacent triangles in the mesh:

$$[[\nabla v_h]] = \nabla v_{h|_{K_1}} \cdot n_1 + \nabla v_{h|_{K_2}} \cdot n_2$$
where $K_1, K_2 \in \mathcal{T}_h$ are such that $K_1 \cap K_2 = F$, and $n_1$ and $n_2$ the exterior pointing normal vectors to $K_1$ and $K_2$. The CIP and GLS parameters are $\gamma_{CIP}, \gamma_{GLS} \ge 0$. Let us remark that $\gamma_{CIP} = \gamma_{GLS} = 0$ is possible because, as we will see, due to the presence of the POD regularization, these primal stabilization terms are optional in the analysis of our method. Several stabilization regimes are studied numerically in Figures 6, 7, and 8, namely the case without GLS and CIP stabilization terms, only a little GLS stabilization, and both GLS and CIP terms with the standard values of $\gamma_{GLS}$ and $\gamma_{CIP}$ chosen.

On the other hand, in order to have some control of the dual variables, it is necessary to include dual stabilization terms: for all $(z_h, y_h), (w_h, x_h) \in W_h^k \times X_h^k$, let

$$s^*[(z_h, y_h), (w_h, x_h)] = s_u^*(z_h, w_h) + s_p^*(y_h, x_h)$$

where

$$s_u^*(z_h, w_h) = \gamma_u^* \langle \nabla z_h, \nabla w_h \rangle_{\Omega}, \quad s_p^*(y_h, x_h) = \gamma_p^* \langle y_h, x_h \rangle_{\Omega}$$

and $\gamma_u^*, \gamma_p^* > 0$. The dual stabilization terms are a regularization in the Tikhonov sense and do not disappear when the mesh size goes to 0, unlike the primal stabilization terms. However, as the dual solution is null at the continuous level, adding these stabilization terms does not change the original problem.

Associated with each of these bilinear forms, $s, s_{GLS}, s_{CIP}, s^*, s_u^*$ and $s_p^*$, we will represent the quadratic forms by $S, S_{GLS}, S_{CIP}, S^*, S_u^*$ and $S_p^*$ respectively.

## 2.4. Lagrangian Formulation

Gathering the different elements introduced in the previous paragraphs, we can now introduce the discrete Lagrangian $\mathcal{L}$: for all $(u_h, p_h) \in V_h^k \times Q_h^{0,k}$, $(z_h, y_h) \in W_h^k \times X_h^k$, let:

$$\mathcal{L}[(u_h, p_h), (z_h, y_h)] = \frac{\gamma_M}{2} \|u_h - u_M\|_{L^2(\Omega)^d}^2 + A[(u_h, p_h), (z_h, y_h)] \\ + \frac{1}{2} S(u_h, p_h) - \frac{1}{2} S^*(z_h, y_h) + \frac{1}{2} R_{POD}(u_h) + \frac{1}{2} R_{POD}^p(u_h, p_h). \quad (9)$$

Computing the Euler-Lagrange equations associated with this Lagrangian leads us to the following discrete problem: find $(u_h, p_h) \in V_h^k \times Q_h^{0,k}$ and $(z_h, y_h) \in W_h^k \times X_h^k$ such that

$$g[(u_h, p_h), (z_h, y_h), (v_h, q_h), (w_h, x_h)] = \gamma_M \langle u_M, v_h \rangle_{\omega} \quad (10)$$

for all $(v_h, q_h) \in V_h^k \times Q_h^{0,k}$ and $(w_h, x_h) \in W_h^k \times X_h^k$ and where

$$g[(u_h, p_h), (z_h, y_h), (v_h, q_h), (w_h, x_h)] := A[(u_h, p_h), (w_h, x_h)] + A[(v_h, q_h), (z_h, y_h)] + s[(u_h, p_h), (v_h, q_h)] \\ - s^*[(z_h, y_h), (w_h, x_h)] + r_{POD}(u_h, v_h) + r_{POD}^p[(u_h, p_h), (v_h, q_h)] + \gamma_M \langle u_h, v_h \rangle_{\omega}.$$

## 3. ANALYSIS

In this section, we show that the discrete problem (10) is well-posed and analyze the convergence of our method. In particular, in Theorem 3.8, we provide bounds on the $L^2$ norms of $u - u_h$ and $p - p_h$ in the whole domain $\Omega$. In contrast with the convergence results for the stabilized FE method for inverse problems, the proof of the convergence result is self-contained in the sense that it does not rely on stability estimates at the continuous level. As we can see in estimate (19), the bound includes, in addition to a term of the form $h^k$, remainder distance terms $\|Q(u)\|_{L^2(\Omega)^d}$ and $\|Q^p(u, p)\|_{L^2(\Omega)}$. Though these terms are expected to be small, the POD method only allows the quantification of the norm of $Q(u)$ in $\omega$. Therefore, in subsection 3.3, we focus on these remainder terms and estimate them with respect to $Q(u)$ in $\omega$. This estimate of logarithmic type is derived from stability estimates proven in [3] for the continuous problem. Finally, the study of the convergence error is extended to the case of perturbed measurements and right hand side in Subsection 3.4 and to the nonhomogeneous Stokes problem in Subsection 3.5.
### 3.1. Well-posedness of the Discrete Problem

We start by defining the following quantity: for $(u_h, p_h) \in V_h^k \times Q_h^{0,k}$ and $(z_h, y_h) \in W_h^k \times X_h^k$,

$$ |||(u_h, p_h), (z_h, y_h)|||^2 := \gamma_M \|u_h\|_{L^2(\omega)^d}^2 + S(u_h, p_h) + R_{POD}(u_h) + R_{POD}^p(u_h, p_h) + S^*(z_h, y_h). \quad (11) $$

**Lemma 3.1.** The quantity given by (11) is a norm over $(V_h^k \times Q_h^{0,k}) \times (W_h^k \times X_h^k)$.

*Proof.* The positivity, homogeneity, and triangular inequality are all clearly verified. For the separation, we assume that (11) is equal to 0.

$$ \bullet S_p^*(y_h) = 0 \implies y_h = 0 \text{ in } \Omega. $$

$$ \bullet S_u^*(z_h) = 0 \implies \nabla z_h = 0 \text{ in } \Omega \text{ and } z_h \in W_h^k \implies z_h = 0 \text{ in } \Omega. $$

$$ \bullet \|u_h\|_{\omega}^2 = 0 \implies u_h = 0 \text{ in } \omega. $$

$$ \bullet R_{POD}(u_h) = 0. \text{ Therefore, } u_h = \sum_{i=1}^{n} \langle u_h, \varphi_i^h \rangle_{\omega} \xi_i = 0 \text{ using the previous point. This implies that } u_h = 0 \text{ in } \Omega. $$

$$ \bullet R_{POD}^p(u_h, p_h) = 0. \text{ So, } p_h = \sum_{i=1}^{n} \langle u_h, \varphi_i^h \rangle_{\omega} \xi_i^p = 0 \text{ in } \Omega. $$

□

In order to prove the existence and uniqueness of the discrete problem (10), we use the following lemma.

**Lemma 3.2.** Let $(u_h, p_h, z_h, y_h) \in V_h^k \times Q_h^{0,k} \times W_h^k \times X_h^k$. Then,

$$ |||(u_h, p_h), (z_h, y_h)||| \le \sup_{\substack{(v_h, q_h) \in V_h^k \times Q_h^{0,k}, \\ (w_h, x_h) \in W_h^k \times X_h^k}} \frac{g[(u_h, p_h), (z_h, y_h), (v_h, q_h), (w_h, x_h)]}{|||(v_h, q_h), (w_h, x_h)|||}. \quad (12) $$

*Proof.* Let $(u_h, p_h) \in V_h^k \times Q_h^{0,k}$ and $(z_h, y_h) \in W_h^k \times X_h^k$. We start by setting $(v_h, q_h) = (u_h, p_h)$ and $(w_h, x_h) = (-z_h, -y_h)$.

$$ g[(u_h, p_h), (z_h, y_h), (u_h, p_h), (-z_h, -y_h)] = S(u_h, p_h) + S^*(z_h, y_h) + R_{POD}(u_h) + R_{POD}^p(u_h, p_h) + \gamma_M \|u_h\|_{\omega}^2 \\ = |||(u_h, p_h), (z_h, y_h)|||^2. $$

We see that

$$ |||(u_h, p_h), (-z_h, -y_h)||| = |||(u_h, p_h), (z_h, y_h)|||. $$

So,

$$ \sup_{\substack{(v_h, q_h) \in V_h^k \times Q_h^{0,k}, \\ (w_h, x_h) \in W_h^k \times X_h^k}} \frac{g[(u_h, p_h), (z_h, y_h), (v_h, q_h), (w_h, x_h)]}{|||(v_h, q_h), (w_h, x_h)|||} \ge \frac{g[(u_h, p_h), (z_h, y_h), (u_h, p_h), (-z_h, -y_h)]}{|||(u_h, p_h), (-z_h, -y_h)|||} \\ = |||(u_h, p_h), (z_h, y_h)||| $$

□

The following proposition is the existence and uniqueness result for the solution to the discrete problem (10).

**Proposition 3.3.** *The discrete problem (10) is well-posed.*

*Proof.* To show the existence and uniqueness of the solution, we apply the Banach-Nečas-Babuška (BNB) Theorem, see [18]. As a direct consequence of Lemma 3.2, the inf-sup condition is verified. Since we are in a case where the solution space and the test space are identical, this property is sufficient to apply the BNB theorem and to deduce that the problem is well-posed. □
**Remark 3.4.** We remark that the well-posedness of the problem is proved without using the primal stabilization terms. In the rest of this section, we will see that they are neither necessary nor affect the stability results.

## 3.2. Error Estimation

In this section, we introduce several lemmas and propositions in the lead-up to our main theoretical result: a convergence error estimation given in Theorem 3.8.

**Remark 3.5.** In the rest of the paper, we use the notation $\lesssim$ to denote an inequality where some constant $C > 0$ is a multiplicative factor of the right hand side. It will change values and unless specified, does not depend on any parameter of the problem. When an inequality does depend on a parameter, we will denote it $\lesssim_{\text{param}}$. If the constant depends on the geometry, it will be denoted $\lesssim_D$ where $D$ is the domain in question.

**Lemma 3.6.** Let $(u, p) \in H^1(\Omega)^d \times L^2(\Omega)$ be the solution to (1) and $((u_h, p_h), (z_h, y_h)) \in (V_h^k \times Q_h^{0,k}) \times (W_h^k \times X_h^k)$ the solution to (10). Therefore,

$$ |||(u_h - u, p_h - p), (z_h, y_h)||| \le C_* \inf_{\substack{\tilde{u}_h \in V_h^k, \\ \tilde{p}_h \in Q_h^{0,k}}} \|\nabla(u - \tilde{u}_h)\|_{L^2(\Omega)^d} + \|p - \tilde{p}_h\|_{L^2(\Omega)} + |||(u - \tilde{u}_h, p - \tilde{p}_h), (0, 0)||| \\ + C_* (R_{\text{POD}}(u)^{1/2} + R_{\text{POD}}^p(u, p)^{1/2}) $$

where $C_*$ is a constant depending on $\gamma_u^*$, $\gamma_p^*$, $\gamma_M$, and $\mu$.

*Proof.* First, we observe that

$$ g[(u_h - u, p_h - p), (z_h, y_h), (v_h, q_h), (w_h, x_h)] = -r_{\text{POD}}(u, v_h) - r_{\text{POD}}^p((u, p), (v_h, q_h)). \quad (13) $$

Let $\tilde{u}_h \in V_h^k$ and $\tilde{p}_h \in Q_h^{0,k}$. By the triangular inequality and Lemma 3.2:

$$ |||(u_h - u, p_h - p), (z_h, y_h)||| \le |||(u_h - \tilde{u}_h, p_h - \tilde{p}_h), (z_h, y_h)||| + |||(\tilde{u}_h - u, \tilde{p}_h - p), (0, 0)||| \\ \le \sup_{\substack{(v_h, q_h) \in V_h^k \times Q_h^{0,k}, \\ (w_h, x_h) \in W_h^k \times X_h^k}} \frac{g[(u_h - \tilde{u}_h, p_h - \tilde{p}_h), (z_h, y_h), (v_h, q_h), (w_h, x_h)]}{|||(v_h, q_h), (w_h, x_h)|||} \\ + |||(\tilde{u}_h - u, \tilde{p}_h - p), (0, 0)|||. \quad (14) $$

Using the bilinearity of $g$ and the Galerkin Relationship (13), we have

$$ g[(u_h - \tilde{u}_h, p_h - \tilde{p}_h), (z_h, y_h), (v_h, q_h), (w_h, x_h)] = -r_{\text{POD}}(u, v_h) - r_{\text{POD}}^p((u, p), (v_h, q_h)) \\ + g[(u - \tilde{u}_h, p - \tilde{p}_h), (0, 0), (v_h, q_h), (w_h, x_h)]. \quad (15) $$

We can bound the two first terms as follows:

$$ r_{\text{POD}}(u, v_h) \lesssim R_{\text{POD}}(u)^{1/2} |||(v_h, q_h), (w_h, x_h)||| \quad (16) $$

and

$$ r_{\text{POD}}^p((u, p), (v_h, q_h)) \lesssim R_{\text{POD}}^p(u, p)^{1/2} |||(v_h, q_h), (w_h, x_h)|||. \quad (17) $$

Now, looking at the last term in (15), we have

$$ g[(u - \tilde{u}_h, p, -\tilde{p}_h), (0, 0), (v_h, q_h), (w_h, x_h)] = A[(u - \tilde{u}_h, p - \tilde{p}_h), (w_h, x_h)] + s[(u - \tilde{u}_h, p - \tilde{p}_h), (v_h, q_h)] \\ + r_{\text{POD}}(u - \tilde{u}_h, v_h) + r_{\text{POD}}^p((u - \tilde{u}_h, p - \tilde{p}_h), (v_h, q_h)) \\ + \langle u - \tilde{u}_h, v_h \rangle_{\omega}. $$

Let us denote by $I_i$, for $1 \le i \le 5$, the five terms in the right hand side. We will estimate them successively. First, developing $I_1$, we obtain:

$$ I_1 = \langle \mu \nabla(u - \tilde{u}_h), \nabla w_h \rangle_{\Omega} - \langle p - \tilde{p}_h, \text{div}(w_h) \rangle_{\Omega} + \langle x_h, \text{div}(u - \tilde{u}_h) \rangle_{\Omega}. $$
Using the Cauchy-Schwarz inequality for these three terms, we find

$$I_1 \lesssim_{\mu, \gamma_u^*, \gamma_p^*} (\|\nabla(u - \tilde{u}_h)\|_{\Omega} + \|p - \tilde{p}_h\|_{\Omega}) \|(v_h, q_h), (w_h, x_h)\|\|\|$$

Moving on to $I_2$, by the Cauchy-Schwarz inequality,

$$I_2 = s[(u - \tilde{u}_h, p - \tilde{p}_h), (v_h, q_h)] \le S[(u - \tilde{u}_h, p - \tilde{p}_h)]^{1/2} \cdot S[(v_h, q_h)]^{1/2} \\ \lesssim \|(u - \tilde{u}_h, p - \tilde{p}_h), (0, 0)\|\|\| \|(v_h, q_h), (w_h, x_h)\|\|\|$$

Similarly, for $I_3$,

$$I_3 = r_{POD}(u - \tilde{u}_h, v_h) \le R_{POD}(u - \tilde{u}_h)^{1/2} \cdot R_{POD}(v_h)^{1/2} \\ \lesssim \|(u - \tilde{u}_h, p - \tilde{p}_h), (0, 0)\|\|\| \|(v_h, q_h), (w_h, x_h)\|\|\|$$

In the same way for $I_4$,

$$I_4 = r_{POD}^p((u - \tilde{u}_h, p - \tilde{p}_h), (v_h, q_h)) \le R_{POD}^p((u - \tilde{u}_h, p - \tilde{p}_h)^{1/2} \cdot R_{POD}^p(v_h, q_h)^{1/2} \\ \lesssim \|(u - \tilde{u}_h, p - \tilde{p}_h), (0, 0)\|\|\| \|(v_h, q_h), (w_h, x_h)\|\|\|$$

Lastly,

$$I_5 \lesssim_{\gamma_M} \|(u - \tilde{u}_h, p - \tilde{p}_h), (0, 0)\|\|\| \|(v_h, q_h), (w_h, x_h)\|\|\|$$

Therefore, gathering (16), (17) and the estimates on $I_i$ for $1 \le i \le 5$ to bound (15), we deduce

$$\sup_{\substack{(v_h, q_h) \in V_h^k \times Q_h^{0,k} \\ (w_h, x_h) \in W_h^k \times X_h^k}} \frac{g[(u_h - \tilde{u}_h, p_h - \tilde{p}_h), (z_h, y_h), (v_h, q_h), (w_h, x_h)]}{\|\|\|(v_h, q_h), (w_h, x_h)\|\|\|} \\ \lesssim_{\mu, \gamma_u^*, \gamma_p^*, \gamma_M} R_{POD}(u)^{1/2} + R_{POD}^p(u, p)^{1/2} + \|\nabla(u - \tilde{u}_h)\|_{\Omega} + \|p - \tilde{p}_h\|_{\Omega} + \|\|(u - \tilde{u}_h, p - \tilde{p}_h), (0, 0)\|\|\|$$

Finally, we use this estimate in (14) and, as $\tilde{u}_h \in V_h^k$ and $\tilde{p}_h \in Q_h^{0,k}$ were chosen arbitrarily, we can pass to the infimum over $(\tilde{u}_h, \tilde{p}_h)$ in $V_h^k \times Q_h^{0,k}$ and conclude. $\square$

Using this lemma, we obtain a first error estimation.

**Proposition 3.7.** Let $(u, p) \in H^1(\Omega)^d \times L^2(\Omega)$ be the solution to (1) and $((u_h, p_h), (z_h, y_h)) \in (V_h^k \times Q_h^{0,k}) \times (W_h^k \times X_h^k)$ the solution to (10). Therefore,

$$\|u - u_h\|_{L^2(\Omega)^d} + \|p - p_h\|_{L^2(\Omega)^d} \le C_* \inf_{\substack{\tilde{u}_h \in V_h^k, \\ \tilde{p}_h \in Q_h^{0,k}}} \|\nabla(u - \tilde{u}_h)\|_{L^2(\Omega)^d} + \|p - \tilde{p}_h\|_{L^2(\Omega)} + \|\|(u - \tilde{u}_h, p - \tilde{p}_h), (0, 0)\|\|\| \\ + C_*(R_{POD}(u)^{1/2} + R_{POD}^p(u, p)^{1/2}) \quad (18)$$

where $C_*$ is a constant that depends on $\gamma_u^*, \gamma_p^*, \gamma_M, \gamma_{POD}$, and $\hat{C}$, which appears in Proposition 2.2.

**Proof.** Using the projection into the space generated by the extended POD modes, we have

$$u - u_h = \mathcal{P}(u - u_h) + \mathcal{Q}(u - u_h) \quad \text{and} \quad p - p_h = \mathcal{P}^p(u - u_h) + \mathcal{Q}^p(u - u_h, p - p_h).$$

So, by the triangular inequality,

$$\|u - u_h\|_{\Omega} \le \|\mathcal{P}(u - u_h)\|_{\Omega} + \|\mathcal{Q}(u - u_h)\|_{\Omega}$$

and

$$\|p - p_h\|_{\Omega} \le \|\mathcal{P}^p(u - u_h)\|_{\Omega} + \|\mathcal{Q}^p(u - u_h, p - p_h)\|_{\Omega}.$$

We notice that:

$$\|\mathcal{Q}(u - u_h)\|_{\Omega} = \sqrt{\frac{1}{\gamma_{POD}}} R_{POD}(u - u_h)^{1/2}$$

and

$$\|\mathcal{Q}^p(u - u_h, p - p_h)\|_{\Omega} = \sqrt{\frac{1}{\gamma_{POD}}} R_{POD}^p(u - u_h, p - p_h)^{1/2}.$$
Therefore,

$$ \|\mathcal{Q}(u - u_h)\|_{\Omega} + \|\mathcal{Q}^{p}(u - u_h, p - p_h)\|_{\Omega} \lesssim_{\gamma_{POD}} |||(u - u_h, p - p_h), (z_h, w_h)|||. $$

Bounding the terms $\mathcal{P}(u - u_h)$ and $\mathcal{P}^p(u - u_h)$ using the Cauchy-Schwarz inequality leads to:

$$ \|\mathcal{P}(u - u_h)\|_{\Omega} \le \|u - u_h\|_{\omega} \sum_{i=1}^{n} \|\varphi_i\|_{\omega} \|\xi_i\|_{\Omega} \le \|u - u_h\|_{\omega} \sum_{i=1}^{n} \|\xi_i\|_{\Omega} $$

and

$$ \|\mathcal{P}^p(u - u_h)\|_{\Omega} \le \|u - u_h\|_{\omega} \sum_{i=1}^{n} \|\varphi_i\|_{\omega} \|\xi_i^p\|_{\Omega} \le \|u - u_h\|_{\omega} \sum_{i=1}^{n} \|\xi_i^p\|_{\Omega} $$

since $(\varphi_i)_{1 \le i \le n}$ is an orthonormal basis with respect to the $L^2(\omega)^d$ norm. So, using Proposition 2.2, we obtain

$$ \|\mathcal{P}(u - u_h)\|_{\Omega} + \|\mathcal{P}^p(u - u_h)\|_{\Omega} \lesssim_{\hat{C}, \gamma_M} |||(u - u_h, p - p_h), (z_h, w_h)|||, $$

which implies

$$ \|u - u_h\|_{\Omega} + \|p - p_h\|_{\Omega} \lesssim_{\hat{C}, \gamma_M, \gamma_{POD}} |||(u - u_h, p - p_h), (z_h, w_h)||| $$

To conclude, we apply Lemma 3.6. □

With the help of interpolations, we can now make a good choice of $\tilde{u}_h$ and $\tilde{p}_h$ to deduce from Proposition 3.7 the following convergence result.

**Theorem 3.8.** For $k \ge 1$, we assume that $(u, p)$ the solution to (1) belongs to $H^{k+1}(\Omega)^d \times H^k(\Omega)$. Let $((u_h, p_h), (z_h, y_h)) \in (V_h^k \times Q_h^{0,k}) \times (W_h^k \times X_h^k)$ be the solution to (10). Then,

$$ \|u - u_h\|_{L^2(\Omega)^d} + \|p - p_h\|_{L^2(\Omega)} \le C_* \left( h^k \|D^{k+1}u\|_{L^2(\Omega)^d} + h^k \|D^k p\|_{L^2(\Omega)} + \|\mathcal{Q}(u)\|_{L^2(\Omega)^d} + \|\mathcal{Q}^p(u, p)\|_{L^2(\Omega)} \right) \quad (19) $$

where $C_*$ is a constant depending on $\gamma_u^*$, $\gamma_p^*$, $\gamma_M$, $\gamma_{POD}$, $\gamma_{GLS}$, $\gamma_{CIP}$, $\mu$, $\hat{C}$, and $\tilde{C}$.

*Proof.* Let $\mathcal{I}_h^k$ be the interpolation operator associated with $\mathbb{P}_k$ finite elements, as defined in [18] in Proposition 1.61. According to the global interpolation theorem given in [18] (Corollary 1.109), we have: for $u \in H^{k+1}(\Omega)$,

$$ \|u - \mathcal{I}_h^k u\|_{L^2(\Omega)^d} + \sum_{m=1}^{k+1} h^m \left( \sum_{K \in \mathcal{T}_h} |u - \mathcal{I}_h^k u|^2_{m,2,K} \right)^{1/2} \le \tilde{C} h^{k+1} \|D^{k+1} u\|_{L^2(\Omega)^d}, \quad (20) $$

where $\tilde{C} > 0$ is a constant. In what follows, we will extensively use inequality (20) in our estimations.

Estimate (19) will be proven by Proposition 3.7 by bounding the infimum term in (18). To do so, we set $\tilde{u}_h = \mathcal{I}_h^k(u) \in V_h^k$ and $\tilde{p}_h = \mathcal{I}_h^k(p) \in Q_h^{0,k}$ and estimate the terms in the infimum for this specific choice of $\tilde{u}_h$ and $\tilde{p}_h$.

First, we have

$$ \|\nabla(u - \tilde{u}_h)\|_{\Omega} + \|p - \tilde{p}_h\|_{\Omega} \le \tilde{C} h^k (\|D^{k+1}u\|_{\Omega} + \|D^k p\|_{\Omega}). \quad (21) $$

Then, for the triple norm, we have

$$ |||(u - \tilde{u}_h, p - \tilde{p}_h), (0, 0)||| = (\gamma_M \|u - \tilde{u}_h\|_{\omega}^2 + S(u - \tilde{u}_h, p - \tilde{p}_h) + R_{POD}(u - \tilde{u}_h) + R_{POD}^{p}(u - \tilde{u}_h, p - \tilde{p}_h))^{1/2}. $$

For the first term, we can write

$$ \|u - \tilde{u}_h\|_{\omega} \lesssim_{\hat{C}} h^{k+1} \|D^{k+1} u\|_{\Omega}. $$

For the global stabilization term $S(u - \tilde{u}_h, p - \tilde{p}_h)$, let us consider the CIP stabilization and the GLS stabilization terms separately. To bound the CIP stabilization term, we introduce the local trace inequality with scaling: for $K \in \mathcal{T}_h$ and $u \in H^1(K)$, we have

$$ h^{1/2} \|u\|_{\partial K} \lesssim \|u\|_{K} + h \|\nabla u\|_{K}. \quad (22) $$
Applying this inequality to $\nabla(u - \tilde{u}_h)$ which belongs to $H^1(K)$, we obtain

$$ 
\begin{aligned}
S_{CIP}(u - \tilde{u}_h) &\lesssim_{\gamma_{CIP}, \mu} \sum_{K \in \mathcal{T}_h} h \|\nabla(u - \tilde{u}_h)\|_{\partial K \setminus \partial \Omega}^2 \\
&\lesssim_{\gamma_{CIP}, \mu} \sum_{K \in \mathcal{T}_h} \|\nabla(u - \tilde{u}_h)\|_K^2 + h^2 \|D^2(u - \tilde{u}_h)\|_K^2.
\end{aligned}
$$

From this inequality, we deduce that

$$ S_{CIP}(u - \tilde{u}_h)^{1/2} \lesssim_{\gamma_{CIP}, \mu} \tilde{C} h^k \|D^{k+1} u\|_{\Omega}. \quad (23) $$

For the GLS stabilization term, we have

$$ 
\begin{aligned}
S_{GLS}(u - \tilde{u}_h, p - \tilde{p}_h)^{1/2} &\lesssim_{\gamma_{GLS}, \mu} \sqrt{\sum_{K \in \mathcal{T}_h} h^2 \| -\Delta(u - \tilde{u}_h) + \nabla(p - \tilde{p}_h) \|_K^2} \\
&\lesssim_{\gamma_{GLS}, \mu} \sqrt{\sum_{K \in \mathcal{T}_h} h^2 \| \Delta(u - \tilde{u}_h) \|_K^2 + h^2 \| \nabla(p - \tilde{p}_h) \|_K^2} \\
&\lesssim_{\gamma_{GLS}, \mu, \tilde{C}} h^k \| D^{k+1} u \|_{\Omega} + h^k \| D^k p \|_{\Omega}.
\end{aligned}
$$

Using (6), the POD terms can be bounded as follows:

$$ 
\begin{aligned}
R_{POD}(u - \tilde{u}_h)^{1/2} &\le \sqrt{\gamma_{POD}} \left( \|u - \tilde{u}_h\|_{\Omega} + \|u - \tilde{u}_h\|_{\omega} \sum_{i=1}^{n} \|\varphi_i\|_{\omega} \|\xi_i\|_{\Omega} \right) \\
&\lesssim_{\hat{C}, \gamma_{POD}} \|u - \tilde{u}_h\|_{\Omega} \\
&\lesssim_{\hat{C}, \gamma_{POD}, \tilde{C}} h^{k+1} \|D^{k+1} u\|_{\Omega}.
\end{aligned}
$$

In the same way,

$$ 
\begin{aligned}
R_{POD}^{p}(u - \tilde{u}_h, p - \tilde{p}_h)^{1/2} &\le \sqrt{\gamma_{POD}} \left( \|p - \tilde{p}_h\|_{\Omega} + \|u - \tilde{u}_h\|_{\omega} \sum_{i=1}^{n} \|\varphi_i\|_{\omega} \|\xi_i^p\|_{\Omega} \right) \\
&\lesssim_{\hat{C}, \gamma_{POD}, \tilde{C}} h^k \|D^k p\|_{\Omega} + h^{k+1} \|D^{k+1} u\|_{\Omega}.
\end{aligned}
$$

Combining all these inequalities leads to the estimate of the triple norm:

$$ |||(u - \tilde{u}_h, p - \tilde{p}_h), (0, 0)||| \lesssim_{\gamma_{POD}, \gamma_{CIP}, \gamma_{GLS}, \mu, \hat{C}, \tilde{C}} h^k \|D^{k+1} u\|_{\Omega} + h^{k+1} \|D^k p\|_{\Omega}. $$

We can then conclude thanks to this estimate and (21). □

**Remark 3.9.** In the proof of Theorem 3.8, we see that if $\gamma_{GLS} = \gamma_{CIP} = 0$, then the results are unchanged and the proof is simplified. The constant $C_*$ no longer depends on $\gamma_{GLS}$ and $\gamma_{CIP}$.

Let us consider the specific case where $u$ can be expressed as its decomposition in the POD space in $\omega$, i.e.

$$ u = \mathcal{P}(u) \text{ on } \omega. \quad (24) $$

According to the unique continuation property for Stokes equations, $Q(u) = 0$ and $Q^p(u, p) = 0$ in $\Omega$. Therefore, by Theorem 3.8 above, the error bound is in $h^k$. This result is therefore of the same nature as the result obtained in [14] for the Laplace equation, when the unknown trace was assumed to belong to a finite dimensional space. In the general case when hypothesis (24) is not satisfied, the error bound still retains a $h^k$ element but adds a corrective term that describes the distance between the solution and the space generated by the extended POD modes. These two “distance terms”, $\|Q(u)\|_{L^2(\Omega)^d}$ and $\|Q^p(u, p)\|_{L^2(\Omega)}$, will be bounded in the next section.
### 3.3. Bounds of Remainder Terms

In this section, we provide upper bounds on the terms in $Q(u)$ and $Q^p(u, p)$ which appear in the right hand side of (19). By the construction of the extended modes $(\xi_i, \xi_i^p)$, these terms are expected to be small. To estimate them quantitatively, we give in the following proposition an upper bound with respect to the local norm $\|Q(u_M)\|_{L^2(\omega)^d}$. Having this local term instead of global norms involving $u$ and $p$ is useful for two reasons. First, if $u_M$ belongs to the database, this local term can be theoretically estimated thanks to classical properties satisfied by the POD modes on $\omega$ (see Remark 3.12). Next, for practical applications, since the measurement $u_M$ is given, it is possible to numerically compute the value of $\|Q(u_M)\|_{L^2(\omega)^d}$.

**Proposition 3.10.** Let us assume that $(u, p)$, the solution to (1), belongs to $H^2(\Omega)^d \times H^1(\Omega)$ and that the database $(\bar{u}^{(j)}, \bar{p}^{(j)})_{1 \le j \le N_{pop}}$ is in $H^2(\Omega)^d \times H^1(\Omega)$.

We also assume that there exists $M > 0$ such that

$$ \|Q(u)\|_{H^2(\Omega)^d} + \|Q^p(u, p)\|_{H^1(\Omega)} \le M. \qquad (25) $$

Then, the pair $(Q(u), Q^p(u, p))$ verifies

$$ \|Q(u)\|_{L^2(\Omega)^d} \lesssim_{\Omega, \omega} \frac{M}{\ln\left(1 + \frac{M}{\|Q(u_M)\|_{L^2(\omega)^d}}\right)} \qquad (26) $$

and

$$ \|Q^p(u, p)\|_{L^2(\Omega)} \lesssim_{\Omega, \omega} \frac{M}{\left(\ln\left(1 + \frac{M}{\|Q(u_M)\|_{L^2(\omega)^d}}\right)\right)^{1/2}}. \qquad (27) $$

**Remark 3.11.** We notice that, as soon as the solution $(u, p)$ and the database $(\bar{u}^{(j)}, \bar{p}^{(j)})$, for $1 \le j \le N_{pop}$, are bounded in $H^2(\Omega)^d \times H^1(\Omega)$, the existence of the bound $M$ in (25) is guaranteed. The bounds of $(\xi_i, \xi_i^p)$ in $H^2(\Omega)^d \times H^1(\Omega)$ are proven by arguments similar to those in the proof of Proposition 2.2.

**Remark 3.12.** According to [21] (Proposition 3.3), if the measurement $u_M$ belongs to the database $\{u_M^{(j)}\}_{1 \le j \le N_{pop}}$, we deduce from the classical mean square error estimate satisfied by the POD basis that:

$$ \sum_{j=1}^{N_{pop}} \|Q(u_M^{(j)})\|_{L^2(\omega)^d}^2 \le \sum_{i=n+1}^{k^*} \lambda_i. $$

where the positive real numbers $(\lambda_i)_{1 \le i \le k^*}$ have been defined in Subsection 2.2.

In several situations, such as certain parametric systems of parametric elliptic PDEs, the eigenvalues might exhibit a fast (sometimes exponential) decay. This was first studied in [23] for a Reduced Basis approximation of an elliptic PDE with one parameter. In a related way, a general study of the Kolmogorov widths decay under holomorphic mappings has been proposed in [15] and an investigation of the Kolmogorov widths decay in parametrized elliptic PDEs has been described in [2].

Thus, in cases where the eigenvalues have an exponential decay, the logarithmic nature of the estimate given by Proposition 3.10 can be compensated for by choosing to increase the number $n$ of POD modes used to approximate our database.

*Proof of Proposition 3.10.* As $(u, p)$ and $(\xi_i, \xi_i^p)$, for all $1 \le i \le n$, are solutions to the Stokes problem, by linearity, we deduce that $(Q(u), Q^p(u, p)) \in H^2(\Omega)^d \times H^1(\Omega)$ also satisfies the Stokes equations.

We can now use results obtained in [3] which studies the unique continuation problem for the Stokes equations at the continuous level. Applying Theorem 1.2 stated in this paper, we directly get the logarithmic bound of $Q(u)$ given in equation (26).

For the $Q^p(u, p)$ term, since $Q^p(u, p)$ belongs to $L_0^2(\Omega)$ and according to Proposition 1.2 from [26], we find

$$ \|Q^p(u, p)\|_{L^2(\Omega)} \le C(\Omega) \|\nabla Q^p(u, p)\|_{H^{-1}(\Omega)^d}. $$
Since $(u, p)$ and $(\xi_i, \xi_i^p)$, for $1 \le i \le n$, satisfy the Stokes equations, we notice that $\nabla Q^p(u, p) = \Delta Q(u)$ in $H^{-1}(\Omega)$. Then,

$$ \|Q^p(u, p)\|_{L^2(\Omega)} \lesssim_{\Omega} \|\Delta Q(u)\|_{H^{-1}(\Omega)}. $$

By definition,

$$ \begin{aligned} \|\Delta Q(u)\|_{H^{-1}(\Omega)} &= \sup_{\varphi \in H_0^1(\Omega)} \frac{\langle \Delta Q(u), \varphi \rangle}{\|\varphi\|_{H_0^1(\Omega)}} \\ &\le \sup_{\varphi \in H_0^1(\Omega)} \frac{\|\nabla Q(u)\|_{\Omega} \|\nabla \varphi\|_{\Omega}}{\|\varphi\|_{H_0^1(\Omega)}} \quad \text{by Green's theorem and the Cauchy-Schwarz inequality} \\ &\le \sup_{\varphi \in H_0^1(\Omega)} \frac{\|\nabla Q(u)\|_{\Omega} \|\varphi\|_{H_0^1(\Omega)}}{\|\varphi\|_{H_0^1(\Omega)}} \\ &= \|\nabla Q(u)\|_{L^2(\Omega)^d}. \end{aligned} $$

Therefore,

$$ \|Q^p(u, p)\|_{L^2(\Omega)^d} \lesssim_{\Omega} \|\nabla Q(u)\|_{L^2(\Omega)^d} \lesssim_{\Omega} \|Q(u)\|_{L^2(\Omega)^d}^{1/2} \|Q(u)\|_{H^2(\Omega)^d}^{1/2}. $$

The first term in the right-hand side, $\|Q(u)\|_{L^2(\Omega)^d}^{1/2}$, can be bounded as a result of (26), while the second term is bounded by $M^{1/2}$. This leads us to (27). $\square$

Theorem 3.8 and the bounds given by Proposition 3.10 bring us to the following corollary.

**Corollary 3.13.** Let $k \ge 1$ and $(u, p) \in H^{k+1}(\Omega)^d \times H^k(\Omega)$ be the solution to (1) and $((u_h, p_h), (z_h, y_h)) \in (V_h^k \times Q_h^{0,k}) \times (W_h^k \times X_h^k)$ the solution to (10). We assume that the database $(\tilde{u}^{(j)}, \bar{p}^{(j)})_{1 \le j \le N_{pop}}$ is in $H^2(\Omega)^d \times H^1(\Omega)$ and that there exists $M > 0$ such that (25) holds. Then,

$$ \begin{aligned} & \|u - u_h\|_{L^2(\Omega)^d} + \|p - p_h\|_{L^2(\Omega)} \\ & \le \begin{cases} C_* (h^k \|D^{k+1}u\|_{L^2(\Omega)^d} + h^k \|D^k p\|_{L^2(\Omega)}) & \text{under hypothesis (24)} \\ C_* (h^k \|D^{k+1}u\|_{L^2(\Omega)^d} + h^k \|D^k p\|_{L^2(\Omega)}) \\ + C_* \left( \frac{M}{\ln\left(1+\frac{M}{\|\nabla Q(u_M)\|_{L^2(\Omega)^d}}\right)} + \frac{M}{\left(\ln\left(1+\frac{M}{\|\nabla Q(u_M)\|_{L^2(\Omega)^d}}\right)\right)^{1/2}} \right) & \text{otherwise.} \end{cases} \end{aligned} $$

where $C_*$ is a constant depending on $\gamma_u^*, \gamma_p^*, \gamma_M, \gamma_{POD}, \gamma_{GLS}, \gamma_{CIP}, \mu, \hat{C}, \tilde{C}, \Omega$ and $\omega$.

Therefore, if the exact solution is not in the space spanned by the extended POD modes, an additional remaining term of logarithmic type appears in the estimate.

## 3.4. Perturbation Analysis

Now let us assume that the measurements in $\omega$ are affected by some perturbation $\delta u$ in $L^2(\omega)^d$. We write this measurement under the form:

$$ u_M^\delta = u + \delta u \text{ in } \omega. $$

We also assume that in the discrete problem, the second member is not exactly equal to 0 but to some small perturbation $\delta f \in L^2(\Omega)^d$. The discrete problem is the following: find $(u_h^\delta, p_h^\delta) \in V_h^k \times Q_h^{0,k}$ and $(z_h^\delta, y_h^\delta) \in W_h^k \times X_h^k$ such that

$$ \begin{aligned} g[(u_h^\delta, p_h^\delta), (z_h^\delta, y_h^\delta), (v_h, q_h), (w_h, x_h)] &= \gamma_M \langle u_M^\delta, v_h \rangle_\omega + \langle \delta f, w_h \rangle_\Omega \\ &= \gamma_M \langle u, v_h \rangle_\omega + \gamma_M \langle \delta u, v_h \rangle_\omega + \langle \delta f, w_h \rangle_\Omega \end{aligned} \quad (28) $$
for all $(v_h, q_h) \in V_h^k \times Q_h^{0,k}$ and $(w_h, x_h) \in W_h^k \times X_h^k$. Adapting the propositions and lemmas from earlier in this section leads us to the convergence error for the perturbed case.

**Theorem 3.14.** Let $k \ge 1$ and $(u, p) \in H^{k+1}(\Omega)^d \times H^k(\Omega)$ be the solution to (1) and $((u_h^\delta, p_h^\delta), (z_h^\delta, y_h^\delta)) \in (V_h^k \times Q_h^{0,k}) \times (W_h^k \times X_h^k)$ the solution to (28). Then,

$$ \|u - u_h^\delta\|_{L^2(\Omega)^d} + \|p - p_h^\delta\|_{L^2(\Omega)} \le C_* (h^k \|D^{k+1}u\|_{L^2(\Omega)^d} + h^k \|D^k p\|_{L^2(\Omega)} + \|\delta u\|_{L^2(\omega)^d} + \|\delta f\|_{L^2(\Omega)^d} \\ + \|Q(u)\|_{L^2(\Omega)^d} + \|Q^p(u, p)\|_{L^2(\Omega)}) $$

where $C_*$ is a constant depending on $\gamma_u^*$, $\gamma_p^*$, $\gamma_M$, $\gamma_{POD}$, $\gamma_{GLS}$, $\gamma_{CIP}$, $\mu$, $\hat{C}$, and $\tilde{C}$.

As with Theorem 3.8, the two “distance terms”, $\|Q(u)\|_{L^2(\Omega)^d}$ and $\|Q^p(u, p)\|_{L^2(\Omega)}$, can be bounded by Proposition 3.10.

**Remark 3.15.** In our analysis of the theoretical problem, we only consider the case when the population data $\{u_M^{(j)}\}_{1 \le j \le N_{pop}}$ are noiseless. The noisy case would make an interesting study in future theoretical works, since the $(\xi_i, \xi_i^p)$ pairs would not necessarily satisfy the Stokes equations and Remark 3.12 would not be applicable. However, in practice, as we will see later in Sections 4.2 and A.3.2, the POD has a denoising effect on the population data, meaning that the first and most significant $\varphi_i$'s and extended POD modes can be considered to be unaffected by noise.

### 3.5. The Case of the Nonhomogeneous Stokes System

In this part, we assume that in the original problem (1) we have a right hand side $f \in H^{k-1}(\Omega)^d$. Thus, we want to find $(u^f, p^f) \in H^{k+1}(\Omega)^d \times (H^k(\Omega) \cap L_0^2(\Omega))$ such that

$$ \left\{ \begin{array}{ll} -\Delta u^f + \nabla p^f = f & \text{in } \Omega, \\ \nabla \cdot u^f = 0 & \text{in } \Omega, \\ u^f|_{\omega} = u_M^f & \text{in } \omega. \end{array} \right. \qquad (29) $$

To solve this problem, we exploit its linearity and decompose it into two problems. First, let $(u_0, p_0) \in H^{k+1}(\Omega)^d \times (H^k(\Omega) \cap L_0^2(\Omega))$ be the solution to the well-posed forward problem:

$$ \left\{ \begin{array}{ll} -\Delta u_0 + \nabla p_0 = f & \text{in } \Omega, \\ \nabla \cdot u_0 = 0 & \text{in } \Omega, \\ u_0 = 0 & \text{in } \partial\Omega. \end{array} \right. \qquad (30) $$

Let $(u^\dagger, p^\dagger) \in H^{k+1}(\Omega)^d \times (H^k(\Omega) \cap L_0^2(\Omega))$ be given by $u^\dagger = u^f - u_0$ and $p^\dagger = p^f - p_0$. By linearity of the Stokes problem, this “corrective solution” $(u^\dagger, p^\dagger)$ must satisfy:

$$ \left\{ \begin{array}{ll} -\Delta u^\dagger + \nabla p^\dagger = 0 & \text{in } \Omega, \\ \nabla \cdot u^\dagger = 0 & \text{in } \Omega, \\ u^\dagger|_{\omega} = u_M^\dagger - u_0 & \text{in } \omega. \end{array} \right. \qquad (31) $$

The database associated with the inverse problem (31) is obtained by the following translation:

$$ u_M^{(k)} = u_M^{f,(k)} - u_{0|\omega} $$

for $1 \le k \le N_{pop}$, where $N_{pop}$ is the number of population data and $\{u_M^{f,(k)}\}_{1 \le k \le N_{pop}}$ is the family of population data associated with the problem (29). We have thus returned to the homogeneous case considered earlier, and can therefore use the method developed in Section 2 to approximate the solution of problem (31) by the solution $(u_h^\dagger, p_h^\dagger)$ of (10) where $u_M$ is replaced by $u_M^f - u_0$ and where the POD basis is associated with the translated measurements $\{u_M^{(k)}\}_{1 \le k \le N_{pop}}$. In addition, the forward problem (30)
is approximated by the $\mathbb{P}_k/\mathbb{P}_{k-1}$ FE method to obtain the discrete solution $(u_{0,h}, p_{0,h}) \in W_h^k \times Q_h^{0,k-1}$.
By [18] (see Section 4.2.5), the error estimate is bounded:

$$ \|u_{0,h} - u_0\|_{L^2(\Omega)^d} + h\|p_{0,h} - p_0\|_{L^2(\Omega)} \le C(\Omega)h^{k+1}(\|u\|_{H^{k+1}(\Omega)^d} + \|p\|_{H^k(\Omega)}). $$

By regularity of the Stokes problem, the right hand side can be further bounded by $\|f\|_{H^{k-1}(\Omega)^d}$. Dividing by $h < 1$ leads to

$$ \|u_{0,h} - u_0\|_{L^2(\Omega)^d} + \|p_{0,h} - p_0\|_{L^2(\Omega)} \le C(\Omega)h^k \|f\|_{H^{k-1}(\Omega)^d}. $$

To return to an approximated solution of the unique continuation problem (29), we define

$$ (u_h^f, p_h^f) := (u_h^\dagger + u_{0,h}, p_h^\dagger + p_{0,h}). $$

Using the estimate above and Theorem 3.8 to estimate $u^\dagger - u_h^\dagger$ and $p^\dagger - p_h^\dagger$, we then obtain the following theorem.

**Theorem 3.16.** Using the notations and definitions from the previous paragraph, we find the following estimation error:

$$ \|u^f - u_h^f\|_{L^2(\Omega)^d} + \|p^f - p_h^f\|_{L^2(\Omega)} \le C_* (h^k \|D^{k+1}u^\dagger\|_{L^2(\Omega)^d} + h^k \|D^k p^\dagger\|_{L^2(\Omega)} + \|Q(u^\dagger)\|_{L^2(\Omega)^d} + \|Q^P(u^\dagger, p^\dagger)\|_{L^2(\Omega)}) \\ + C(\Omega)h^k \|f\|_{H^{k-1}(\Omega)^d} $$

where $C_*$ is a constant depending on $\gamma_u^*$, $\gamma_p^*$, $\gamma_M$, $\gamma_{POD}$, $\gamma_{GLS}$, $\gamma_{CIP}$, $\mu$, $\hat{C}$, and $\tilde{C}$.

Let us notice that, as in the previous theorems, the last two terms of the estimate can be bounded using Proposition 3.10.

## 4. IMPLEMENTATION OF THE NUMERICAL METHOD

The method used to numerically solve the unique continuation problem can be broken down into three steps. The first is the generation of a POD basis in $\omega$ from a database. This database is given, and it is a datum of the present formulation. In our work, it is synthetically generated by performing direct numerical simulations. Next, the population information has to be extended to the entirety of $\Omega$ by solving the inverse problem (5). Finally, the unique continuation problem can be set up and solved. We choose to work with $\mathbb{P}_1$ finite elements in all numerical sections.

### 4.1. Setup and Data Generation

We start by constructing the domain $\Omega$ and its subdomain $\omega$. In our numerical tests, which will be discussed in Section 5, we work in 2 and 3 dimensional space, letting $\Omega$ be a tube and a rectangular box respectively. Let $\Gamma_{in}$ and $\Gamma_{out}$ denote the inlet and outlet boundaries, and let $\Gamma_{wall}$ be the walls of the domain such that $\partial\Omega = \Gamma_{in} \cup \Gamma_{wall} \cup \Gamma_{out}$. An example of the tube and region of observation in 2D is seen in Figure 1.

<figure><img src="image_2.png" /><figcaption>FIGURE 1. Example of 2D Domain</figcaption></figure>
The unique continuation problem to be solved numerically is written

$$ 
\left\{ 
\begin{array}{lll}
-\mu\Delta u + \nabla p & = 0 & \text{in } \Omega, \\
\nabla \cdot u & = 0 & \text{in } \Omega, \\
u & = 0 & \text{on } \Gamma_{\text{wall}}, \\
\mu \frac{\partial u}{\partial n} - p \vec{n} & = 0 & \text{on } \Gamma_{\text{out}}, \\
u|_{\omega} & = u_M & \text{in } \omega.
\end{array}
\right. 
\qquad (32) 
$$

In regard to the unique continuation problem studied in Sections 2 and 3, we notice that we are in a more favorable case here since boundary data are partially known (no-slip conditions along the wall $\Gamma_{\text{wall}}$ and homogeneous Neumann conditions on the outlet $\Gamma_{\text{out}}$) and the only unknown boundary condition is on the inlet $\Gamma_{\text{in}}$. For this problem, the variational formulation can be easily adapted (we refer to Remark 2.1) and the theoretical results previously obtained remain valid.

Before numerically solving this problem, we first need to generate synthetic data measurements for the database in $\omega$. The synthetic data are obtained by solving a family of direct problems over $\Omega$ and storing the values of the velocity at the nodes in $\omega$. In these direct problems, boundary conditions on the inlet $\Gamma_{\text{in}}$ are provided and denoted $b_{\alpha}^{\text{pop}}$ (more details on these conditions are given in the numerical results section 5). Since we are considering mixed boundary conditions in non smooth domains, the solution pair $(u, p)$ will be sought in the space $H^1(\Omega)^d \times L^2(\Omega)$ and, to avoid discontinuities on the prescribed boundary data, we will impose the inlet boundary condition $b_{\alpha}^{\text{pop}}$ to be equal to zero at the corners in the 2D case and equal to zero along the edges of the inlet face in 3D.

The forward problem can thus be written:

$$ 
\left\{ 
\begin{array}{lll}
-\mu\Delta u + \nabla p & = 0 & \text{in } \Omega, \\
\nabla \cdot u & = 0 & \text{in } \Omega, \\
u & = 0 & \text{on } \Gamma_{\text{wall}}, \\
\mu \frac{\partial u}{\partial n} - p \vec{n} & = 0 & \text{on } \Gamma_{\text{out}}, \\
u(x) & = b_{\alpha}^{\text{pop}}(x) & \text{on } \Gamma_{\text{in}}.
\end{array}
\right. 
\qquad (33) 
$$

The forward problem is solved numerically, with the inclusion of the primal stabilization terms, as described in the presentation of the problem in Section 2, because $\mathbb{P}_1/\mathbb{P}_1$ finite elements are used. We choose $\gamma_{\text{GLS}} = \gamma_{\text{CIP}} = 0.1$, which is standard in literature.

Once the direct problem is solved for different inlet boundary conditions and the velocity solutions in $\omega$ are stored, the next step is the computation of the POD. We compute the Singular Value Decomposition (SVD) of the measurements with respect to the norm given by the mass matrix of the mesh. As mentioned in Subsection 2.2, we only keep a certain number of the POD basis vectors, whose associated singular values are not too small. By taking the logarithm of the singular values, we can see how many basis vectors to keep. This is demonstrated graphically in Figure 2 for the 2D case.

## 4.2. Extension of the POD Basis Functions from $\omega$ to $\Omega$

One of the challenging aspects of the implementation of our method is the numerical computation of the $n$ extended modes $(\xi_i, \xi_i^p)_{1 \le i \le n}$, which are defined as the solutions to the unique continuation problem (5). What might first come to mind would be to generate them by implementing the classical stabilized FE method for inverse problems (which is associated with the Lagrangian (9) without the last two terms coming from the population data). However, as is theoretically shown in, e.g. [13] and [7], the accuracy of the extended modes would be degraded near the boundary. As a consequence, since these errors propagate to the global optimization problem, the benefit brought by the information that the solution is close to a finite space is lost by the fact that the finite dimensional space is not precisely enough described.

In this subsection, we present an alternative method for generating the extended modes $(\xi_i, \xi_i^p)_{1 \le i \le n}$. It relies on an optimal recovery formulation which, in very general terms, consists in the estimation of
features of an unknown function from limited information. The general framework of the method is presented in [27]. Optimal recovery methods can be applied to PDEs, as in [6], where it was applied to a data assimilation problem for parametric PDEs or in [5], where it was used to solve a Poisson problem with missing boundary conditions. No numerical analysis is available for the optimal recovery strategy we implement, but numerically, we observe that it accurately generates extended modes.

Our optimal recovery method is based on the following approach. We assume that we have access to a finite number, $N_\omega$, of point evaluations of each POD basis function $\varphi_i$. We will denote $\{X_j\}_{1\le j\le N_\omega}$ the set of points in $\omega \subset \mathbb{R}^d$ where these measurements are taken. For $1 \le i \le n$, we will construct the extended POD mode pair $(\xi_i, \xi_i^P)$ as minimizing the Stokes energy functional

$$ \mathcal{E}(u) = \frac{\mu}{2} \int_{\Omega} |\nabla u|^2 $$

and such that $\xi_i$ is equal to $\varphi_i$ in the points $\{X_j\}_{1\le j\le N_\omega}$. We define the spaces

$$ V = \{\xi \in H^1(\Omega)^d, \nabla \cdot \xi = 0 \text{ in } \Omega\} \times L^2(\Omega) $$

and

$$ \kappa_i = \{(\xi, \xi^P) \in V : \xi(X_j) = \varphi_i(X_j), \forall 1 \le j \le N_\omega\}. $$

We remark that by considering a finite number of point evaluations, we lose the unique continuation framework and its properties. The i-th extended POD mode pair is now given by:

$$ (\xi_i, \xi_i^P) = \underset{(\xi, \xi^P) \in \kappa_i}{\arg\inf} \mathcal{E}(\xi). $$

Let us comment on this strategy. First, in this method, we impose $\xi_i(X_j)$ to coincide exactly with $\varphi_i(X_j)$ for all $\{X_j\}_{1\le j\le N_\omega}$. This approach can be hazardous if the function $\varphi_i$ is affected too much by noise. Here, we take advantage of the fact that the POD method acts like a filter for noise. In our problem, we observe that the first singular values are separated from the rest by a spectral gap (this is illustrated, for instance, in Figure 2), and this allows us to distinguish the first few modes with the highest energy from the lower energy modes, which are affected by noise and thrown out. This is another way in which this optimal recovery method distances itself from the original unique continuation problem (1), as in the latter, the measurements would have to be noiseless, which is not always the case. In Appendix A.3.2, we show that when the POD basis functions are calculated using noisy measurements, for the first few (and most important) extended modes, there is very little difference between them and the modes generated from noiseless measurements.

In what follows, to simplify the notation, we omit the index $i$ of the POD mode. To numerically solve the extension problem, we will introduce a family of solutions to the Stokes problem and construct the extended POD modes as linear combinations of these solutions. To do so, for $N_b \in \mathbb{N}^*$, we consider a set of functions $\{b_k\}_{1\le k\le N_b}$ in $H_0^1(\Gamma_{in})^d$ and we define $(u_k^b, p_k^b)$ as the solution to:

$$ \left\{ \begin{array}{ll} -\mu \Delta u_k^b + \nabla p_k^b = 0 & \text{in } \Omega, \\ \nabla \cdot u_k^b = 0 & \text{in } \Omega, \\ u_k^b = 0 & \text{on } \Gamma_{wall}, \\ \mu \frac{\partial u_k^b}{\partial n} - p_k^b \vec{n} = 0 & \text{on } \Gamma_{out}, \\ u_k^b = b_k & \text{on } \Gamma_{in}. \end{array} \right. \qquad (34) $$

To discretize this problem, let $(u_{h,k}^b, p_{h,k}^b)$ be the $\mathbb{P}_1/\mathbb{P}_1$ FE approximated solution of the forward problem, for $1 \le k \le N_b$. We construct the discretized, extended POD modes $(\xi_h, \xi_h^P)$ by:

$$ \xi_h = \sum_{k=1}^{N_b} a_k u_{h,k}^b \quad \text{and} \quad \xi_h^P = \sum_{k=1}^{N_b} a_k p_{h,k}^b $$

where $a = (a_1, \dots, a_{N_b})$ is a vector of coefficients to be determined. By linearity of the Stokes problem, we know that $(\xi_h, \xi_h^P)$ is also a solution to the Stokes problem with homogeneous Dirichlet boundary
conditions on $\Gamma_{wall}$ and Neumann conditions on $\Gamma_{out}$. We take for measurement points the nodes of the mesh in $\omega$. Letting $V_{N_b} = \text{Span}\{(u_{h,k}^b, p_{h,k}^b)\}, 1 \le k \le N_b\}$, the discretized $\kappa_{N_b}$ associated with some approximated POD mode $\varphi_h$ is given by:

$$ \kappa_{N_b} = \{(\xi_h, \xi_h^p) \in V_{N_b} : \xi_h(X_j) = \varphi_h(X_j), \forall 1 \le j \le N_\omega\}. \quad (35) $$

The resulting minimization problem is then:

$$ (\xi_h^*, \xi_h^{p,*}) = \arg\inf_{(\xi_h, \xi_h^p) \in \kappa_{N_b}} \mathcal{E}_h(\xi_h), \quad (36) $$

where $\mathcal{E}_h$ is the discretized Stokes energy functional. In this problem, the number of constraints, which we will denote $N_{data}$, is twice the number of measurement points in the 2D case and 3 times the number of measurement points in the 3D case, because the $\xi_h$'s are vectoriel, like the velocity. In other words, $N_{data} = 2N_\omega$ in 2D and $N_{data} = 3N_\omega$ in 3D. As a result, in order to be able to determine the solution of the minimization problem, we need to be in a case where $N_b \ge N_{data}$. Since the value of $N_b$ is directly linked to the number of degrees of freedom along the inlet boundary, if the mesh of the inlet boundary $\Gamma_{in}$ is too coarse, we refine it, calculate the extended POD modes on the mesh with the refined boundary, and then interpolate them back to the original mesh. In Appendix A.3.1, we provide practical details on the choice of the family of the inlet functions $\{b_k\}_{1 \le k \le N_b}$.

In the discrete Lagrangian associated with this problem, we include the Galerkin least squares and CIP stabilization terms. In order to write this discretized Lagrangian, let us introduce some notation. Let $a = (a_1, \dots, a_{N_b})^T$ be the unknown. We write $\xi_h = U_h^b a$ where $U_h^b$ is the matrix of columns $u_{h,k}^b$. Similarly, let the extended POD mode for the pressure be $\xi_h^p = P_h^b a$ where $P_h^b$ is the matrix of column $p_{h,k}^b$. Let $K_{u,h}$ be the stiffness matrix associated with the velocity, multiplied by the constant $\mu$, $K_{GLS,h}$ the matrix corresponding to the Galerkin Least Squares stabilization term, and finally $K_{CIP,h}$ the CIP stabilization matrix. We introduce the matrix $H_h$ to map the restriction from $\Omega$ to $\omega$. In Subsection 5.1, we will consider numerical test cases where the mesh on $\omega$ is coarser than the mesh on $\Omega$, in which case $H_h$ will also be an interpolation matrix bringing the nodes on $\Omega$ to the nodes on $\omega$ (see also Remark 5.1). The Lagrangian is then for all $(a, \lambda) \in \mathbb{R}^{N_b} \times \mathbb{R}^{N_{data}}$:

$$ \mathcal{L}(a, \lambda) = \frac{1}{2} a^T [\mu U_h^{bT} K_{u,h} U_h^b + P_h^{bT} K_{GLS,h} P_h^b + U_h^{bT} K_{CIP,h} U_h^b] a - \lambda^T (\varphi_h - H_h U_h^b a). \quad (37) $$

Let us notice that the Galerkin least squares term used above is $P_h^{bT} K_{GLS,h} P_h^b$, though by definition, the velocity $U_h^b$ should also appear. In fact, as we are working with $\mathbb{P}_1$ elements, the Laplacian of the velocity is equal to zero. For ease of notation, we only write the pressure terms. Also, the stiffness matrix $K_{u,h}$ used above is obtained with Neumann boundary conditions. This is not exactly the same stiffness matrix used in the stabilization of the unique continuation problem, as this one requires some Dirichlet conditions.

To simplify the next few paragraphs, we introduce some additional notations.

Let $A^* = \mu U_h^{bT} K_{u,h} U_h^b + P_h^{bT} K_{GLS,h} P_h^b + U_h^{bT} K_{CIP,h} U_h^b$, and $B = H_h U_h^b$. Differentiating the discrete Lagrangian $\mathcal{L}$ leads to the following linear system:

$$ \begin{bmatrix} A^* & -B^T \\ -B & 0 \end{bmatrix} \begin{bmatrix} a \\ \lambda \end{bmatrix} = \begin{bmatrix} 0 \\ -\varphi_h \end{bmatrix}. \quad (38) $$

This linear system needs some additional attention as, in practice, it is extremely ill-conditioned. In order to solve the optimization problem (37), we write its dual formulation and use a proximal iteration, as described in [24]. This is detailed in Appendix A.1.

Finally, denoting $a^*$ the solution of this optimization problem, the extended POD mode is given by $(\xi_h, \xi_h^p) = (\mathcal{I}U_h^b a^*, \mathcal{I}P_h^b a^*)$ where $\mathcal{I}$ is the interpolation operator from the FE space defined on the mesh with a refined inlet boundary to the FE space defined on the current mesh without the refinement.
### 4.3. Unique Continuation

Once equipped with the extended POD modes and an individual's measurement, $u_M$, we are able to solve the discrete problem (10) for that individual and implement it in FreeFem++.

In order to facilitate the implementation of the POD terms, we orthonormalize the extended modes using the Gram-Schmidt method to obtain $(\tilde{\xi}_h^i)_{1\le i\le n}$ and $(\tilde{\xi}_h^{p,i})_{1\le i\le n}$. To the system matrix comprised of the variational formulation and stabilization terms, we add the matrix

$$ \gamma_{\text{POD}}(M_u - M_u \Xi_h^T \Xi_h^M u) + \gamma_{\text{POD}}(M_p - M_p \Xi_h^p \Xi_h^{pM}), \quad (39) $$

which is the discretization of the modified version of terms $r_{\text{POD}}(u_h, v_h) + r_{\text{POD}}^p((u_h, p_h), (v_h, q_h))$, see (7) and (8). For numerical reasons, we prefer to work with the orthonormalized modes, replacing $\langle u_{M,h}, \varphi_h^i \rangle_{\omega} \xi_h^i$ by $\langle u_h, \tilde{\xi}_h^i \rangle_{\Omega} \tilde{\xi}_h^i$ in (7) and $\langle u_{M,h}, \varphi_h^i \rangle_{\omega} \xi_h^{p,i}$ by $\langle u_h, \tilde{\xi}_h^i \rangle_{\Omega} \tilde{\xi}_h^{p,i}$ in (8). $M_u$ and $M_p$ are the mass matrices for $u$ and $p$ respectively. $\Xi_h$ and $\Xi_h^p$ are the matrices whose columns are the discrete, orthonormalized, extended POD modes of the velocity and the pressure.

Additionally, we use the POD modes to denoise and interpolate the measurement $u_M$, replacing it by

$$ u_{M,h}^{\text{POD}} = \sum_{i=1}^{n} \langle u_{M,h}, \xi_h^i \rangle_{\omega} \xi_h^i. \quad (40) $$

The interpolation by the POD modes is especially useful when the measurements are scarce, meaning the spatial resolution of the measurements is coarser than the desired resolution of the final solution. This interpolation allows us to solve the unique continuation problem over the resolution of the $\xi_h^i$'s, which is finer than that of the measurements. A discussion of this can be found in Section 5.1.

**Remark 4.1.** In equation (40), we are denoising the measurements through an approximation by their truncated POD decomposition. This truncation does not lose any valuable information due to a well-chosen $n$, the number of POD modes. In fact, we choose a big enough $n$ so that the POD decomposition is a good approximation of the measurements $u_M$ (the underlying hypothesis being that the measurements are close enough to the POD space, otherwise this method would not be appropriate), but not too big, otherwise the noise is amplified. In practice, by calculating the singular values of the POD decomposition, see Figures 2 and 12, we fix $n$ as the number of singular values above the spectral gap.

The unique continuation problem can be solved directly for small systems with built-in solvers (for example, UMFPACK or MUMPS) and, for larger systems, we considered a strategy based on Woodbury matrix identity in order to reduce the computational cost both in terms of memory and number of operations. This method is detailed in Appendix A.2.

## 5. NUMERICAL RESULTS

In this work, we study two test cases. First, we reconstruct the velocity and pressure in a 2D tube and then, in a 3D rectangular box in which there is a spherical obstacle. In both cases, the viscosity is chosen to be $\mu = 0.035$ Poise, and $P_1$ finite elements are used. Two types of data regimes are considered: one where the data measurements are available on each node of the mesh in $\omega$ and one where there is a fixed number of scarce data measurements. In addition to the noiseless case, two levels of noise are applied to the data: 1%, and 5% (relative to the $L^2(\omega)^d$ norm). This noise is additive and unbiased, calculated by drawing samples from normal distributions and then added point by point to each component of the data measurements. The numerical results obtained are compared with those of the “Classical Method”: the stabilized FE method for inverse problems in which no population data are used.

### 5.1. 2D Case

Starting with the 2D case, let $\Omega = [0, 6] \times [-1, 1]$ be a rectangle and $\omega = [1, 3] \times [-\frac{1}{2}, \frac{1}{2}]$ the region of observation, as shown in Figure 1. We construct the synthetic data by solving the direct problem
(33) over $\Omega$ and saving the values of the velocity at the nodes in $\omega$. The non-homogeneous Dirichlet conditions on the inlet, $b_{\alpha}^{pop}$, are of the form

$$b_{\alpha}^{pop}(x, y) = \begin{bmatrix} (1 - y^2)(\alpha_0 + \alpha_1 y + \alpha_2 y^2 + \alpha_3 y^3) \\ 0 \end{bmatrix}$$

with $\alpha = (\alpha_0, \alpha_1, \alpha_2, \alpha_3)^T$. We randomly generate the coefficients following the Uniform Law, taking $\alpha_0 \sim U[1, 2]$ and $\alpha_1, \alpha_2, \alpha_3 \sim U[-0.4, 0.4]$. Solving this direct problem for 100 $\alpha$'s gives us the database of size 100.

Noise is added to the database by drawing from the distribution $N(0, 0.01)$ in the 1% case and from $N(0, 0.05)$ in the 5% case. A randomly generated sample is then added to each measurement, to both the $x$ and $y$ components. Noise is also added to the individual measurements $u_M$ by generating samples from the distribution $N(0, 0.01)$ and then scaling to either 1% or 5% with respect to the $L^2(\omega)^d$ norm.

In Figure 2, the logarithm of the singular values is shown for the noiseless, 1% noise, and 5% noise cases in 2D. We observe that it makes sense to retain the first four POD modes in the noiseless case. We keep the first three in the 1% noise case and the first two in the 5% noise case.

<figure><img src="image_3.png" /><figcaption>FIGURE 2. Logarithm of Singular Values</figcaption></figure>

The POD basis vectors on $\omega$ are then extended to $\Omega$ using the method presented in Subsection 4.2. Some details on the implementation of this step, in particular, the choice of the inlet boundary conditions $b_k$ to generate the solutions $(u_k^b, p_k^b)$ of (34), for $1 \le k \le N_b$, are specified in Appendix A.3.1.

In our tests, we consider two data regimes: a fine case and a coarse case. In the fine case, the data measurements are taken on the nodes of the mesh in $\omega$. This means that if the mesh size decreases, the number of measurements increases. In the coarse case, however, the spatial resolution of the data is fixed at a coarse resolution to mimic a realistic scenario. When the problem is solved over a mesh finer than that of the data, the measurements must be interpolated to the finer mesh. We remark that the coarse data case can thus be seen as a perturbation of the fine data case, as the interpolation to the finer mesh introduces errors. Figure 3 is an example of the fine and coarse data cases. On the left, in Figure 3a, representing the fine case, we see a mesh over $\omega$, where each node corresponds to a data point. On the right, Figure 3b is the coarse data case.
<figure><img src="image_4.png" /><figcaption>FIGURE 3. Comparison Between Fine and Coarse Meshes on $\omega$</figcaption></figure>

As previously mentioned, we work with the projection of the data measurements $u_M$ in the POD space, as given by equation (40). This not only denoises the measurements, but when dealing with coarse data, it interpolates the measurements from the coarse resolution to the finer resolution of the $\xi_h^i$ modes. Thus, even though the scalar product $\langle u_{M,h}, \xi_h^i \rangle_\omega$ is calculated over the coarse mesh, the measurements used to solve the unique continuation problem are defined on a refined mesh.

**Remark 5.1.** When the POD modes are extended following the method presented in Subsection 4.2, if we are in the coarse data regime, then $\varphi_h$ is defined on the coarse mesh whereas $\xi_h$ is defined on the fine mesh. For this reason, the mapping $H_h$ was introduced in the discrete Lagrangian for the extension problem in Equation (37). The mapping brings $\xi_h$ defined on the fine mesh to the coarse mesh so that the condition $\xi_{h|\omega} = \varphi_h$ holds on the coarse mesh.

Our numerical tests will study the relative errors corresponding to the estimation errors of the velocity, pressure, and on occasion, the velocity and pressure combined: $\|u - u_h\|_{L^2(\Omega)^d}$, $\|p - p_h\|_{L^2(\Omega)}$, and $\|u - u_h\|_{L^2(\Omega)^d} + \|p - p_h\|_{L^2(\Omega)}$. As our method has several free parameters, such as $\gamma_M$ and $\gamma_{POD}$, we will start by analyzing their effect on the estimation errors in the following paragraphs. We will also study how the errors behave as the mesh-size decreases in both the fine and coarse data regimes and for different regularization techniques. These tests will be performed using the Poiseuille analytical solution as the target solution, which is given by

$$u(x, y) = \begin{bmatrix} 1 - y^2 \\ 0 \end{bmatrix} \quad \text{and} \quad p(x, y) = \mu(L - 2x),$$

where $L$ corresponds to the length of our tube ($L = 6$). Letting the measurements $u_M$ be the values of the velocity solution at nodes in $\omega$ and choosing the dual stabilization parameters consistent with literature ($\gamma_u^* = \gamma_p^* = 0.1$), we begin the tests, starting with a study of $\gamma_M$ and $\gamma_{POD}$.

### 5.1.1. Effects of $\gamma_M$ and $\gamma_{POD}$

Two free parameters playing an important role in our method are $\gamma_M$ and $\gamma_{POD}$. The data penalty constant $\gamma_M$ can be interpreted as a measure of how much we trust our data: the higher the $\gamma_M$ chosen, the more we impose the solution to match the data measurements in $\omega$. The POD parameter, $\gamma_{POD}$, can be similarly interpreted: a higher value indicates that we want the solution to be closer to the space generated by the population data. In these tests, we fix $\gamma_M$ to be equal to 10, 100, or 1000 and let $\gamma_{POD}$ vary in [0.05, 0.5, 5., 50., 500, 5000], plotting the errors in a log-log plot in Figures 4 and 5. Three levels of noise are compared: a noiseless case, a 1%, and a 5% noise case. We let $\gamma_{GLS} = 0.001$ and $\gamma_{CIP} = 0$ (see Subsection 5.1.2 for a discussion of this choice of stabilization) and choose a medium-sized mesh, see Mesh 3 in Table 1, and the fine data regime.
<figure><img src="image_5.png" /><figcaption>FIGURE 4. $\|u - u_h\|_{L^2(\Omega)^d}$ Relative Error vs $\gamma_{POD}$</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIGURE 5. $\|p - p_h\|_{L^2(\Omega)}$ Relative Error vs $\gamma_{POD}$</figcaption></figure>

Figure 4 shows that when $\gamma_M$ is high (equal to 100 or 1000), varying $\gamma_{POD}$ has little effect on the error, for all levels of noise. On the other hand, when $\gamma_M = 10$, increasing $\gamma_{POD}$ has a detrimental effect on the error in the noiseless case: the errors reach the same height as in the 5% noise case. In fact, if $\gamma_M = 10$, then intuitively, we are not placing much trust in our individual data measurements. Since the subspace generated by the population data is constructed using data of the same type (same noise, distribution, and points of measurements), it would not make sense to choose a high value of $\gamma_{POD}$, i.e. placing more trust in these measurements than in the individual's. We see the same phenomenon in Figure 5, this time for the pressure. However, the pressure errors show less stability than the velocity errors, and $\gamma_{POD} = 5$ seems to be a sweet spot across all values of $\gamma_M$ tested. The parameter values that suit all levels of noise are $\gamma_M = 10$ and $\gamma_{POD} = 5$, which will be used later in Subsection 5.1.3.

### 5.1.2. Effects of the Mesh-size and Regularization Techniques

In the main theoretical result, Theorem 3.8, the rate of convergence of our method for $P_1/P_1$ FE is in $h$ with a corrective term denoting the distance between the solution and the space generated by population data. In this section, we solve the unique continuation problem for different mesh-sizes to study the numerical rate of convergence. Additionally, we compare different regularization techniques. First is the “Classical Method” baseline, in which the unique continuation problem is solved without using any population data and with the GLS and CIP stabilization terms. The parameters $\gamma_{GLS} = \gamma_{CIP} = 0.1$ are chosen consistent with what is the standard in literature. Let us specify that this baseline, otherwise known as the stabilized FE method for inverse problems, uses no population information, neither a soft penalty in the Lagrangian formulation, nor to interpolate the measurements. The next regularization technique is the “With POD, Standard Stabilization” case, where the population data is included by adding the term (39) to the Lagrangian formulation and by using the projected measurements (40).
In this case, we keep the GLS and CIP stabilization terms with the standard values of parameters ($\gamma_{GLS} = \gamma_{CIP} = 0.1$). Next, we remove all the traditional stabilization terms (GLS and CIP), just leaving the POD as the only regularization: the “With POD, No Stabilization” method. Finally, we include the population data as before and add in a little GLS stabilization, taking $\gamma_{GLS} = 0.001$ and $\gamma_{CIP} = 0$ for the “With POD, Some Stabilization” case.

We compare the errors in the fine and coarse data regimes in the noiseless case. For the “Classical Method”, the parameter $\gamma_M = 1000$ is chosen, which is standard in literature ($\gamma_{POD}$ does not appear). In all subsequent tests, both in 2D and 3D, this will be the value chosen for the “Classical Method”. For the three “With POD” methods, we keep $\gamma_M = 1000$, which is consistent with our study of $\gamma_M$ in Subsection 5.1.1. For the “With POD, Standard Stabilization” method, we let $\gamma_{POD} = 5000$; due to the heavier stabilization, we find a higher value of $\gamma_{POD}$ performs better. For the “With POD, Some Stabilization” and “With POD, No Stabilization” techniques, we let $\gamma_{POD} = 5$. The Figures 6, 7, and 8 show the relative errors in a log-log plot as the mesh-size decreases. We start with the $\|u - u_h\|_{L^2(\Omega)^d}$ and $\|p - p_h\|_{L^2(\Omega)}$ error in the fine and coarse data regimes in Figures 6a and 6b.

<figure><img src="image_7.png" /><figcaption>FIGURE 6. Relative Errors of $\|u - u_h\|_{L^2(\Omega)^d}$ and $\|p - p_h\|_{L^2(\Omega)}$ as a Function of Mesh-Size</figcaption></figure>

In Figures 7a and 7b, the separate relative errors of $\|u - u_h\|_{L^2(\Omega)^d}$ are shown in both data regimes.
<figure><img src="image_8.png" /><figcaption>FIGURE 7. Relative Errors of $||u - u_h||_{L^2(\Omega)^d}$ as a Function of Mesh-Size</figcaption></figure>

Finally, in Figures 8a and 8b, the separate relative errors of $||p - p_h||_{L^2(\Omega)}$ can be found.

<figure><img src="image_9.png" /><figcaption>FIGURE 8. Relative Errors of $||p - p_h||_{L^2(\Omega)}$ as a Function of Mesh-Size</figcaption></figure>

In all the cases tested, we see that the population data improve results, no matter the regularization technique chosen, both for the velocity and the pressure. In particular, the “With POD, No Stabilization” and “With POD, Some Stabilization” methods have basically the same rate of convergence and the lowest errors. The “With POD, Standard Stabilization” case has the fastest rate of convergence for the pressure, interestingly enough, but its errors are higher than the other two “With POD” methods. Due to the pressure error being higher than the velocity error, this regularization regime also provides the fastest convergence for the combined velocity and pressure case. More generally, we observe that, when the relative error is below a certain threshold, the curves corresponding to the “With POD” method have a downward trend when the mesh size decreases that appear to be of logarithmic nature. We understand this behavior to be related to the mixed rate of convergence obtained in Corollary 3.13. In this theoretical result, a rate of convergence in $h$ is expected if the solution is in the population-generated POD space.
Otherwise, the $h$ rate is mixed with a logarithmic rate which depends on $\|Q(u_M)\|_{L^2(\omega)^d}$, where $u_M$ are the continuous measurements in $\omega$. Due to the different discretization errors, in particular in the approximation of the extended POD modes, it is reasonable to infer that this logarithmic term leads to a term of the form $|\log(h)|^{-\alpha}$ for some $\alpha > 0$. It is then very likely that, when the error goes below a certain level, the logarithmic rate takes over.

**Remark 5.2.** Though it seems that the “With POD, No Stabilization” and “With POD, Some Stabilization” test cases perform exactly alike, when noise is added, we observe that it is the “With POD, Some Stabilization” method that performs the best. For this reason, in Subsection 5.1.1 and in the rest of the numerical tests, this will be the population-enriched method used to solve the unique continuation problem.

Next, let us compare the number of measurement points available in $\omega$ in the fine and coarse data regimes. In Table 1, we show the percentage of measurements points in $\omega$ out of the total number of nodes in $\Omega$, with Mesh 1 denoting the coarsest mesh and Mesh 6 the finest.

<table><caption>TABLE 1. % of dof in $\omega$ out of $\Omega$</caption><thead><tr><th>Mesh #</th><th>Fine Data Regime</th><th>Coarse Data Regime</th></tr></thead><tbody><tr><td>1</td><td>22.4%</td><td>22.4%</td></tr><tr><td>2</td><td>21.3%</td><td>6.0%</td></tr><tr><td>3</td><td>21.0%</td><td>2.7%</td></tr><tr><td>4</td><td>20.8%</td><td>1.6%</td></tr><tr><td>5</td><td>20.7%</td><td>1.0%</td></tr><tr><td>6</td><td>20.6%</td><td>0.70%</td></tr></tbody></table>

In Table 1, as the mesh is refined, the amount of information we have access to through the measurement points in $\omega$ in the coarse data regime is significantly less than in the fine data regime. In figures 6b, 7b, and 8b, we see the effects. The “Classical Method” generally shows a decrease in error as the mesh size decreases, but for the velocity, this decrease is very slight. We see that the population-enriched method helps to work around this lack of information by the inclusion of the POD terms (39) as well as the interpolation of the coarse data measurements to the fine mesh by (40), thereby reducing interpolation errors. Due to the database, there is much improvement in the errors and even some improvement in the rate of convergence.

The error curves for the Poiseuille pressure solutions 8a and 8b show the relative errors over the entirety of $\Omega$. We are also interested in the performance on certain subdomains of $\Omega$, in particular how the inclusion of the population information improves the pressure estimation near the inlet boundary $\Gamma_{in}$. In Figure 9, we plot the error in absolute value between the Poiseuille analytical pressure solution and the approximate pressure. On the left, in 9a, is the solution given by the “With POD, Some Stabilization” method and on the right, in 9b is the solution from the “Classical Method”. We place ourselves in the coarse data regime without noise, choosing $\gamma_M = 1000$ and $\gamma_{POD} = 5$ as before. The errors are scaled to the “with population data” case.
<figure><img src="image_10.png" /><figcaption>FIGURE 9. Pressure Error, Coarse Data, 0% Noise</figcaption></figure>

We see that most of the error in both cases is concentrated along $\Gamma_{in}$. With population information, the maximum error in absolute value is 0.0036, compared to 0.017 in the “Classical Method” case. Figure 10 shows the same comparison, but in the 5% noise case, using the same $\gamma_M$ and $\gamma_{POD}$ values as in the noiseless case.

<figure><img src="image_11.png" /><figcaption>FIGURE 10. Pressure Error, Coarse Data, 5% Noise</figcaption></figure>

In the “With POD, Some Stabilization” case, the error is concentrated along the inlet boundary and reaches at most 0.029 in absolute value. In the classical case, however, we see the more of the effect of the noise on the measurements in $\omega$, as well as errors along $\Gamma_{in}$, reaching 0.051 in absolute value. This illustrates the de-noising capabilities of the POD.

### 5.1.3. Statistical Study of the Errors

In this section, we solve the unique continuation problem for multiple sets of measurements in order to perform a statistical analysis of the relative estimation errors. The sets of measurements used are “out of database, in range”. When generating the database, we generate some extra solutions, drawing their inlet polynomial coefficients from the same distribution used for the population data. These extra solutions, not being part of the database, are called “out of database, in range”. In our computation of means and standard deviations, we use 32 sets of individual measurements in $\omega$. In the noisy cases, we generate 32 perturbations and solve the unique continuation problem for each individual set of measurements and each perturbation. This signifies that the means and standard deviations of $32^2 = 1024$ velocity and pressure errors are computed in the 1% and 5% noise cases, while in the noiseless case, the means and standard deviations are of 32 total errors. Table 2 compares the “With POD, Some Stabilization” method with the “Classical Method” in the fine data regime. Table 3 compares the same two methods, but in the coarse data regime. An intermediate-sized mesh is used, called Mesh 3 in Table 1. In all tests, the parameters $\gamma_M = 10$ and $\gamma_{POD} = 5$ are chosen, which are values that suit all three levels of noise by Subsection 5.1.1.
<table><caption>TABLE 2. Mean and Standard Deviations of Relative $L^2(\Omega)$ Errors for the Fine Data Regime</caption><thead><tr><th rowspan="2">Noise</th><th colspan="2">||u - u<sub>h</sub>||<sub>L<sup>2</sup>(&Omega;)<sup>d</sup></sub></th><th colspan="2">||p - p<sub>h</sub>||<sub>L<sup>2</sup>(&Omega;)</sub></th></tr><tr><th>Classical Method</th><th>With Population Data</th><th>Classical Method</th><th>With Population Data</th></tr></thead><tbody><tr><td>0%</td><td>1.20% ± 0.379</td><td>0.258% ± 0.103</td><td>2.44% ± 0.477</td><td>0.676% ± 0.179</td></tr><tr><td>1%</td><td>1.23% ± 0.368</td><td>0.682% ± 0.319</td><td>2.45% ± 0.477</td><td>1.102% ± 0.382</td></tr><tr><td>5%</td><td>1.70% ± 0.327</td><td>2.128% ± 1.032</td><td>2.530% ± 0.555</td><td>2.140% ± 0.862</td></tr></tbody></table>

<table><caption>TABLE 3. Mean and Standard Deviations of Relative $L^2(\Omega)$ Errors for the Coarse Data Regime</caption><thead><tr><th rowspan="2">Noise</th><th colspan="2">||u - u<sub>h</sub>||<sub>L<sup>2</sup>(&Omega;)<sup>d</sup></sub></th><th colspan="2">||p - p<sub>h</sub>||<sub>L<sup>2</sup>(&Omega;)</sub></th></tr><tr><th>Classical Method</th><th>With Population Data</th><th>Classical Method</th><th>With Population Data</th></tr></thead><tbody><tr><td>0%</td><td>1.46% ± 0.286</td><td>0.308% ± 0.095</td><td>2.925% ± 0.427</td><td>0.738% ± 0.159</td></tr><tr><td>1%</td><td>1.559% ± 0.31</td><td>0.878% ± 0.328</td><td>2.957% ± 0.453</td><td>1.477% ± 0.504</td></tr><tr><td>5%</td><td>2.913% ± 0.518</td><td>1.786% ± 0.65</td><td>3.253% ± 0.785</td><td>2.163% ± 0.822</td></tr></tbody></table>

Tables 2 and 3 reveal that, for the noiseless and 1% noise cases, using the population-enriched method leads to average velocity and pressure errors that are much lower than those found with the “Classical Method”. The errors in the 1% noise case are reduced by about half, but in the noiseless case, by well over half. Even using coarse data, we are able to outperform the “Classical Method” in the fine data regime. As for the 5% noise case, the average velocity error in the fine data regime is not lower than in the “Classical Method”. However, the population-enriched method does yield lower errors on average for both velocity and pressure in the coarse data regime and for pressure in the fine data regime for the 5% noise case.

## 5.2. 3D Case

Let $R = [0, 6] \times [0, 4]$, $\times [0, 4]$ be a rectangular box and let $S$ be a sphere inside of radius 1 and center $(2, 2, 2)$. We define $\Omega = R \setminus S$. The inlet face is $\{0\} \times [0, 4] \times [0, 4]$ and the region of observation $\omega$ is $[0.5, 2] \times [0, 4] \times [0, 4]$, which splits the sphere in half. For this geometry, we only consider a coarse data regime which puts us in a situation where the velocity is only known at 3.7% of the total nodes of the fine mesh. From this, we reconstruct the velocity and pressure everywhere over a finer mesh. Figure 11a shows the fine mesh of $\Omega$ and Figure 11b, one individual’s sparse measurements in $\omega$, with the coarse resolution.
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>FIGURE 11. Domain in 3D</figcaption></figure>

The generation of the database has a similar setup as in the 2D case, except now, the polynomial inlet boundary condition is of the form

$$b_{\alpha}^{\text{pop}}(x, y, z) = \begin{bmatrix} yz(4-y)(4-z)(\alpha_0 + \alpha_1 y + \alpha_2 z + \alpha_3 yz + \alpha_4 y^2 + \alpha_5 z^2 + \alpha_6 zy^2 + \alpha_7 yz^2 + \alpha_8 y^3 + \alpha_9 z^3) \\ 0 \\ 0 \end{bmatrix}$$

The coefficients were drawn with the following Uniform Laws: $\alpha_0, \alpha_1, \alpha_2 \sim U([0, 10])$, $\alpha_3, \alpha_4, \alpha_5 \sim U([0, 5])$, and $\alpha_6, \alpha_7, \alpha_8, \alpha_9 \sim U([0, 2.5])$. Once the database is generated, we compute the POD. The number of POD modes to be kept is 10 for all noise cases, as seen in Figure 12. The noise is generated in the same way as in the 2D case.

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

<figure><img src="image_15.png" /><figcaption>FIGURE 12. Logarithm of Singular Values in 3D</figcaption></figure>

As for the 2D case, details on the extension of the POD modes are given in Appendix A.3.1. The unique continuation problem is solved with a little GLS stabilization, taking $\gamma_{GLS} = 0.001$ and $\gamma_{CIP} = 0$, as well as $\gamma_M = 1000$ and $\gamma_{POD} = 5000$. Figures 13 and 14 compare the velocity and pressure approximations with the exact solutions. Table 4 shows the relative errors for the methods with and without population data.
<figure><img src="image_16.png" /><figcaption>FIGURE 13. Velocity, 0% Noise</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>FIGURE 14. Pressure, 0% Noise</figcaption></figure>

<table><caption>TABLE 4. Relative Errors for 0% Noise Case</caption><thead><tr><th>L<sup>2</sup>(Ω) Error</th><th>Classical Method</th><th>With Population Data</th></tr></thead><tbody><tr><td>Velocity u</td><td>26.55%</td><td>0.49%</td></tr><tr><td>Pressure p</td><td>41.66%</td><td>7.18%</td></tr></tbody></table>

The velocity error decreases by over a factor of 10 with the inclusion of population data, and the pressure decreases by about a factor of 6. We remark that there is no zone, like the inlet, where the
pressure error is higher than elsewhere - it is roughly 7% over all of $\Omega$. When noise is added, we obtain the errors displayed in Table 5.

<table><caption>TABLE 5. Relative Errors for 1% and 5% Noise Cases</caption><thead><tr><th rowspan="2">L<sup>2</sup>(Ω) Error</th><th colspan="2">1% Noise</th><th colspan="2">5% Noise</th></tr><tr><th>Classical</th><th>With Population Data</th><th>Classical</th><th>With Population Data</th></tr></thead><tbody><tr><td>Velocity u</td><td>26.55%</td><td>0.52%</td><td>26.55%</td><td>1.90%</td></tr><tr><td>Pressure p</td><td>41.68%</td><td>7.12%</td><td>41.68%</td><td>7.11%</td></tr></tbody></table>

The 1% noise case performs similarly to the noiseless case, whereas the velocity error increases when the 5% noise is added. To explain why the error shows little change overall when noise is added for the “with population data” method, Table 6 in Appendix A.3.2 shows how close the extended POD modes in noisy cases are to the extended modes in the noiseless case. The first few modes have relative errors well under 1%. Since they are the modes that carry the most information, it makes sense that the unique continuation solutions are similar with and without noise. The 9th and 10th noisy modes are significantly different to the noiseless ones, but as they do not bring a lot of information to the span, they do not really change the solutions.

## 6. CONCLUSION

In this work, we include population data in a unique continuation problem, which to our best of knowledge, is novel. Though the population data are not any richer than the individual data used in the classical approach to the problem, their inclusion brings improvements both theoretically and numerically. Our method only depends on one additional, but realistic hypothesis: the existence of a database. This work also introduces a method to extend population information from a subdomain to the entire domain using optimal recovery.

Specifically, the inclusion of population data allows us to obtain a global error estimation for the discrete velocity and pressure over the entire domain of resolution. We also find a faster rate of convergence than in previous works, which is supported numerically by the 2D test case. Additionally, the population data improve results when the measurements are scarce. The best example of this is in the 3D test case when given only 3.7% of the velocity values as measurements, the velocity error decreases by over a factor of 10 and the pressure by about a factor of 6 when the population data are included in the resolution.

Here, we work with a linear, stationary model - the Stokes equations. In the future, it would be interesting to consider non-stationary models or nonlinear models, such as the time-dependent Navier-Stokes equations. A few starting points for future work are [1], [17], and [20], which study data assimilation for the Navier-Stokes equations in the case of single measurements. Both [1] and [17] focus on steady-state blood flow models with non-trivial geometries, such as a carotid bifurcation in [1]. This paper [1] incorporates velocity measurements in the resolution of the forward problem using a least-squares finite element formulation and presents numerical tests in 3D. In [17], the authors seek to minimize a data misfit functional under the constraint of the PDE discretized by the FEM and Tikhonov-regularized. The influence of noisy and sparse data measurements are both studied, including the case when the measurements are not taken on the nodes of the mesh. Finally, [20] provides an error analysis for the velocity when coarse measurements are incorporated into the Navier-Stokes equations (with missing initial conditions) through a nudging parameter.

In the future, we would like to expand on our study of the numerical method for the POD mode extension. To do so, we will first focus on the Poisson equation, which will make a comparison and analysis of different approximation methods more accessible. Additionally, given that the extension method in our work exploits the linearity of the Stokes problem, further work would need to be done to extend the POD modes in nonlinear settings. Regarding the analysis of the convergence rate, a next
step could also be to consider approximate, discrete POD modes instead of the exact, continuous ones that we use here. Moreover, it would be interesting to consider multiple databases. For example, in the practical application of determining blood velocity and pressure in an artery, we could realistically consider two groups of patients: one healthy and one sick. In short, there are several potential avenues for future work studying unique continuation problems enriched by population data.

## REFERENCES

[1] Solveigh Averweg, Alexander Schwarz, Carina Schwarz, and Jörg Schröder. 3d modeling of generalized newtonian fluid flow with data assimilation using the least-squares finite element method. *Computer Methods in Applied Mechanics and Engineering*, 392:114668, 2022.

[2] Markus Bachmayr and Albert Cohen. Kolmogorov widths and low-rank approximations of parametric elliptic pdes. *Mathematics of Computation*, 86(304):701–724, 2017.

[3] Mehdi Badra, Fabien Caubet, and Jérémi Dardé. Stability estimates for navier-stokes equations and application to inverse problems. *Discrete and Continuous Dynamical Systems - Series B*, 21(8):2379–2407, September 2016.

[4] Andrea Ballerini. Stable determination of an immersed body in a stationary Stokes fluid. *Inverse Problems*, 26(12):125015, 25, 2010.

[5] Peter Binev, Andrea Bonito, Albert Cohen, Wolfgang Dahmen, Ronald DeVore, and Guergana Petrova. Solving pdes with incomplete information, 2023.

[6] Peter Binev, Albert Cohen, Wolfgang Dahmen, Ronald DeVore, Guergana Petrova, and Przemysław Wojtaszczyk. Data assimilation in reduced modeling. *SIAM/ASA J. Uncertain. Quantif.*, 5(1):1–29, January 2017.

[7] Muriel Boulakia, Erik Burman, Miguel A. Fernández, and Colette Voisembert. Data assimilation finite element method for the linearized navier–stokes equations in the low reynolds regime. *Inverse Problems*, 36(8):085003, August 2020.

[8] Muriel Boulakia, Anne-Claire Egloffe, and Céline Grandmont. Stability estimates for the unique continuation property of the stokes system and for an inverse boundary coefficient problem. *Inverse Problems*, 29(11):115001, sep 2013.

[9] Erik Burman. Stabilized finite element methods for nonsymmetric, noncoercive, and ill-posed problems. part i: Elliptic equations. *SIAM Journal on Scientific Computing*, 35(6):A2752–A2780, 2013.

[10] Erik Burman, Deepika Garg, and Janosch Preuss. Data assimilation finite element method for the linearized navier-stokes equations with higher order polynomial approximation. *ESAIM: Mathematical Modelling and Numerical Analysis*, 58(1):223–245, January 2024.

[11] Erik Burman and Peter Hansbo. Stabilized nonconforming finite element methods for data assimilation in incompressible flows. *Math. Comp.*, 87(311):1029–1050, 2018.

[12] Erik Burman, Peter Hansbo, and Mats G Larson. Solving ill-posed control problems by stabilized finite element methods: an alternative to tikhonov regularization. *Inverse Problems*, 34(3):035004, January 2018.

[13] Erik Burman and Lauri Oksanen. *Weakly Consistent Regularisation Methods for Ill-Posed Problems*, pages 171–202. Springer International Publishing, Cham, 2018.

[14] Erik Burman and Lauri Oksanen. Finite element approximation of unique continuation of functions with finite dimensional trace. *Mathematical Models and Methods in Applied Sciences*, 34(10):1809–1824, July 2024.

[15] Albert Cohen and Ronald DeVore. Kolmogorov widths under holomorphic mappings. *IMA Journal of Numerical Analysis*, 36(1):1–12, 2016.

[16] Jim Douglas and Todd Dupont. *Interior Penalty Procedures for Elliptic and Parabolic Galerkin Methods*, page 207–216. Springer Berlin Heidelberg, 1976.

[17] Marta D'Elia, Mauro Perego, and Alessandro Veneziani. A variational data assimilation procedure for the incompressible navier-stokes equations in hemodynamics. *Journal of Scientific Computing*, 52(2):340–359, 2012.
[18] Alexandre Ern and Jean-Luc Guermond. *Theory and Practice of Finite Elements*. Springer New York, 2004.

[19] Caroline Fabre and Gilles Lebeau. Prolongement unique des solutions de l'équation de stokes. *Communications in Partial Differential Equations*, 21(3-4):573-596, January 1996.

[20] Bosco García-Archilla and Julia Novo. Error analysis of fully discrete mixed finite element data assimilation schemes for the navier-stokes equations. *Advances in Computational Mathematics*, 46(4), 2020.

[21] K. Kunisch and S. Volkwein. Galerkin proper orthogonal decomposition methods for a general equation in fluid dynamics. *SIAM Journal on Numerical Analysis*, 40(2):492-515, 2002.

[22] Ching-Lung Lin, Gunther Uhlmann, and Jenn-Nan Wang. Optimal three-ball inequalities and quantitative uniqueness for the Stokes system. *Discrete Contin. Dyn. Syst.*, 28(3):1273-1290, 2010.

[23] Yvon Maday, Anthony T Patera, and Gabriel Turinici. A priori convergence theory for reduced-basis approximations of single-parameter elliptic partial differential equations. *Journal of Scientific Computing*, 17:437-446, 2002.

[24] Neal Parikh and Stephen Boyd. *Proximal algorithms*, volume Vol. 1, No. 3 (2013) 123-231. Foundations and Trends in Optimization, 2013.

[25] Farzin Shakib, Thomas J.R. Hughes, and Zdeněk Johan. A new finite element formulation for computational fluid dynamics: X. the compressible euler and navier-stokes equations. *Computer Methods in Applied Mechanics and Engineering*, 89(1):141-219, 1991. Second World Congress on Computational Mechanics.

[26] Roger Temam. *On the Theory and Numerical Analysis of the Navier-Stokes Equations*. Universite Paris XI, 1973.

[27] Peter R Turner. *Numerical analysis, Lancaster 1984*. Lecture notes in mathematics. Springer, Berlin, Germany, 1985 edition, June 1985.

# APPENDIX A.

## A.1. Linear System (38) and Iterative Proximal Algorithm

The linear system (38) requires a little linear algebra help before it is solvable. To simplify notation, we write it as

$$MX = Y.$$

The block matrices that make up $M$ are $A^* \in \mathbb{R}^{N_b \times N_b}$ and $B \in \mathbb{R}^{N_{data} \times N_b}$. We observe that $A^*$ is a symmetric, positive definite matrix, therefore invertible. $B$, however, is not. In fact, in numerical tests, we found $rank(B) \ll N_{data}$ and most of the singular values to be extremely small. To fix this, we filter out the extremely small singular values in $B$. Taking the singular value decomposition $B = USV^T$, let

$$\nu = \#\{\sigma : \sigma > \varepsilon\}$$

where $\sigma$ are the singular values of $B$ and $\varepsilon > 0$ is a threshold parameter. We now truncate the SVD of $B$ to obtain $\hat{B} = \hat{U}\hat{S}\hat{V}^T$ where $\hat{U} \in \mathbb{R}^{N_{data} \times \nu}$ is the first $\nu$ columns of $U$, $\hat{S} \in \mathbb{R}^{\nu \times \nu}$ the diagonal matrix, whose elements are the first $\nu$ elements of $S$, and so on. To ensure that $\varphi_h \in \mathbb{R}^{N_{data}}$ is in the span of the columns of $\hat{U}$, we work with

$$\tilde{\varphi}_h = \hat{U}\hat{U}^T\varphi_h \in \mathbb{R}^{N_{data}}.$$

Numerically, we can verify that this modified POD mode is close to $\varphi_h$. This leads us to a modified Lagrangian defined in $\mathbb{R}^{N_b} \times \mathbb{R}^{N_{data}}$ by

$$\tilde{\mathcal{L}}(a, \lambda) = \frac{1}{2}a^T A^* a + \lambda^T (\tilde{\varphi}_h - \hat{B}a).$$

By differentiating with respect to the Lagrange multiplier $\lambda$, we obtain

$$\hat{B}a = \tilde{\varphi}_h.$$
Differentiating the Lagrangian with respect to a leads to the following relation

$$a = A^{-1} \hat{B}^T \lambda,$$

which we substitute for the value of a in the Lagrangian. This leads us to the dual formulation:

$$\tilde{\mathcal{L}}^*(\lambda) = \lambda^T \tilde{\varphi}_h - \frac{1}{2} \lambda^T \hat{B} A^{-1} \hat{B}^T \lambda.$$

Finding the critical point of the Lagrangian $\tilde{\mathcal{L}}$ corresponds to maximizing the dual functional $\tilde{\mathcal{L}}^*$. To solve this maximization problem, we use an iterative proximal algorithm. In the next few paragraphs, we introduce this algorithm and proximal operators, based on Sections 1.1, 2.3, and 4.1 of the monograph by Parikh and Boyd, see [24].

Let $f: \mathbb{R}^m \to \mathbb{R} \cup \{+\infty\}$ be a closed and proper convex function. The scaled proximal operator, defined between $\mathbb{R}^m \to \mathbb{R}^m$ is

$$\operatorname{prox}_{\alpha f}(v) = \arg \min_{x \in \mathbb{R}^m} \{f(x) + \frac{1}{2\alpha} \|x - v\|_2^2\} \quad \forall v \in \mathbb{R}^m,$$

where $\alpha > 0$ is the value of each step. An important property of the proximal operator is that its fixed point is a minimizer of $f$. The sequence $(x_k)_{k \ge 0}$, defined by

$$x_{k+1} = \operatorname{prox}_{\alpha f}(x_k),$$

converges to the minimizer of $f$.

Let us now calculate the proximal operator for our function $-\tilde{\mathcal{L}}^*$, which is quadratic with a positive leading coefficient, so closed and convex. For $v \in \mathbb{R}^{N_{\text{data}}}$, we have

$$\operatorname{prox}_{-\alpha \tilde{\mathcal{L}}^*}(v) = \arg \min_{\lambda \in \mathbb{R}^{N_{\text{data}}}} \left\{ \frac{1}{2} \lambda^T \hat{B} A^{-1} \hat{B}^T \lambda - \lambda^T \tilde{\varphi}_h + \frac{1}{2\alpha} \|\lambda - v\|_2^2 \right\}.$$

Taking the derivative of the terms in the brackets, the minimizer, which we denote $\lambda^*$, must verify:

$$\hat{B} A^{-1} \hat{B}^T \lambda^* - \tilde{\varphi}_h + \frac{1}{\alpha}(\lambda^* - v) = 0.$$

This is equivalent to

$$(I + \alpha \hat{B} A^{-1} \hat{B}^T) \lambda^* = v + \alpha \tilde{\varphi}_h$$

The proximal operator is thus

$$\operatorname{prox}_{-\alpha \tilde{\mathcal{L}}^*}(v) = (I + \alpha \hat{B} A^{-1} \hat{B}^T)^{-1} (v + \alpha \tilde{\varphi}_h)$$

for all $v \in \mathbb{R}^{N_{\text{data}}}$ as $I + \alpha \hat{B} A^{-1} \hat{B}^T$ is invertible. Given the solution for the $k$-th iteration, $\lambda_k$, the next solution in the sequence is the solution to the linear problem,

$$(I + \alpha \hat{B} A^{-1} \hat{B}^T) \lambda_{k+1} = \lambda_k + \alpha \tilde{\varphi}_h,$$

which can be solved directly with a Python solver. Returning to the relationship between $\lambda$ and $a$, we find

$$a_k = A^{-1} \hat{B}^T \lambda_k.$$

We continue iterating the proximal algorithm until $\frac{\|\tilde{\varphi}_h - \hat{B} a_k\|}{\|\tilde{\varphi}_h\|} < \varepsilon$, where $\varepsilon > 0$ is a certain tolerance, or until a maximum number of iterations has been completed.
## A.2. A Linear System Solving Method Based on Woodbury Matrix Identity

To supplement Subsection 4.3 on the numerical resolution of the unique continuation problem, let us recall the terms making up the linear system. The right hand side system matrix is composed of two elements, the variational formulation with the stabilization terms, as in (37), and the POD matrix, as in (39). Letting $\Lambda$ denote the matrix composed of the variational formula and stabilization terms, $X$ the solution, and $b$ the right hand side, the system is:

$$[\Lambda + \gamma_{\text{POD}}(M_u - M_u \Xi_h \Xi_h^T M_u) + \gamma_{\text{POD}}(M_p - M_p \Xi_h^p \Xi_h^{pT} M_p)]X = b. \quad (41)$$

Remark that, by construction, the system matrix is invertible. The matrix $\Lambda$ is symmetric, positive definite, and sparse. The mass matrices $M_u$ and $M_p$ are also symmetric, positive definite, and sparse. The other two matrices, constructed with the extended POD modes, are given in low-rank format. For large systems (as, for instance, when dealing with 3D test cases), it is not feasible to store them as full matrices. Instead, we would like to devise a resolution strategy keeping them in low-rank format, in order to reduce the computational burden both in terms of memory and number of operations. The method presented hereafter makes use of the Woodbury matrix identity.

Let $\tilde{A} = \Lambda + \gamma_{\text{POD}}(M_u + M_p)$ be the sparse part of the left hand side matrix. Let $U = M_u \Xi_h$, $V = M_p \Xi_h^p$, and $W = [U, V]$. We notice that $W$ has $2n$ columns, where $n$ is the number of POD modes. The extended POD matrix terms can be noted $\gamma_{\text{POD}}WW^T$. The system is now written

$$(\tilde{A} - \gamma_{\text{POD}}WW^T)X = b$$

or equivalently

$$(\tilde{A} + WCW^T)X = b$$

where $C = -\gamma_{\text{POD}}I$, with $I$ the identity matrix. We solve the above system using the Woodbury matrix identity:

$$(\tilde{A} + WCW^T)^{-1} = \tilde{A}^{-1} - \tilde{A}^{-1}W(C^{-1} + W^T\tilde{A}^{-1}W)^{-1}W^T\tilde{A}^{-1}.$$

The solution is then

$$X = \tilde{A}^{-1}b - \tilde{A}^{-1}W(C^{-1} + W^T\tilde{A}^{-1}W)^{-1}W^T\tilde{A}^{-1}b.$$

From a numerical point of view, the above-written solution is obtained as follows:

(1) Solve the linear system $\tilde{A}y^{(0)} = b$.

(2) Solve the $2n$ linear systems $\tilde{A}y^{(i)} = w^{(i)}$, $1 \le i \le 2n$, where $w^{(i)}$ is the $i$-th column of the matrix $W$.

(3) Assemble the matrix $D = C^{-1} + W^T\tilde{A}^{-1}W$, whose components are:

$$D_{ij} = -\frac{1}{\gamma_{\text{POD}}} \delta_{ij} + [w^{(i)}]^T y^{(j)}, \quad 1 \le i, j \le 2n.$$

(4) Compute the vector $\beta \in \mathbb{R}^{2n}$, whose components are defined as:

$$\beta_i = [w^{(i)}]^T y^{(0)}, \quad 1 \le i \le 2n.$$

(5) Solve the linear system of size $2n$:

$$D\eta = \beta.$$

(6) Assemble the system solution:

$$X = y^{(0)} + \sum_{i=1}^{2n} \eta_i y^{(i)}.$$

Remark that the computational cost of the method is dominated by the resolution of the $2n+1$ sparse linear systems in the two first steps of the method. If the resolution of the sparse linear systems is performed by using a direct method, the proposed method is a direct method.
### A.3. Numerical Resolution Supplements

#### A.3.1. Extension of the POD Basis: Choice of the Family of Solutions $(u_k^b, p_k^b)$

In this paragraph, we give details of the implementation of the extension of the POD basis, in both the 2D and 3D cases. In particular, we specify the choice of the boundary functions $b_k$ used to construct the solutions $\{(u_k^b, p_k^b)\}_{1 \le k \le N_b}$ of problem (34). To have a good approximation of arbitrary inlet boundary data, we will select these functions as part of the Fourier basis on the inlet boundary.

In the 2D case, we consider the following functions which belong to $H_0^1(\Gamma_{in})^2$:

$$b_k^{(1)}(y) = \left[ \sin\left(\frac{\pi k(y+1)}{2}\right) \right] \quad \text{and} \quad b_k^{(2)}(y) = \left[ \sin\left(\frac{\pi k(y+1)}{2}\right) \right]$$

for $k$ varying in $\{1, \dots, N_{\Gamma_{in}}\}$ where $N_{\Gamma_{in}}$ corresponds to number of internal degrees of freedom on $\Gamma_{in}$. As we have both $b_k^{(1)}$ and $b_k^{(2)}$ to deal with, this means that $N_b = 2N_{\Gamma_{in}}$.

With the above inlet boundary condition, we can generate the family of solutions and solve the optimal recovery problem, taking the threshold parameter $\epsilon = 10^{-3}$ for every test case except the 5% noise case with coarse data, where $\epsilon = 0.25$ is needed. The step size in the iterative proximal algorithm $\alpha = 1$ is chosen.

In the 3D case, we consider the inlet boundary conditions given by the $H_0^1(\Gamma_{in})^3$ functions:

$$b_{i,j}^{(1)}(y, z) = \begin{bmatrix} \sin\left(\frac{\pi iy}{4}\right) \cdot \sin\left(\frac{\pi jz}{4}\right) \\ 0 \\ 0 \end{bmatrix}, \quad b_{i,j}^{(2)}(y, z) = \begin{bmatrix} 0 \\ \sin\left(\frac{\pi iy}{4}\right) \cdot \sin\left(\frac{\pi jz}{4}\right) \\ 0 \end{bmatrix},$$

and

$$b_{i,j}^{(3)}(y, z) = \begin{bmatrix} 0 \\ 0 \\ \sin\left(\frac{\pi iy}{4}\right) \cdot \sin\left(\frac{\pi jz}{4}\right) \end{bmatrix}$$

for $i, j$ varying in $\{1, \dots, N_{\Gamma_{in}}\}$ where $N_{\Gamma_{in}}$ corresponds to the approximate number of internal degrees of freedom along an edge of the inlet boundary face. Then, $N_b = 3N_{\Gamma_{in}}^2$. In the 3D case, we are working on a single fixed mesh and take $N_{\Gamma_{in}} = 25$. Since the number of degrees of freedom in $\omega$ is equal to 591, the necessary inequality $N_b = 3 \cdot 25^2 = 1875 > N_{data} = 3 \cdot 591 = 1773$ holds and we do not need to introduce a refined mesh to tackle the minimization problem (36).

Once the family of solutions is generated, the POD basis modes are extended to the entire domain by the same method as in the 2D case. The parameter $\epsilon = 10^{-3}$ and the step-size $\alpha = 1$ are selected. Once the extended modes are obtained, the unique continuation problem is solved using the Woodbury matrix identity as detailed in Appendix A.2, with $\gamma_{POD} = 5000$ and $\gamma_M = 1000$.

#### A.3.2. 3D Numerical Results

To understand the similarity in errors between noiseless and noisy cases, let us compare the extended POD modes. Table 6 shows the $L^2(\Omega)^3$ errors between the $\xi$ mode in the 0% Noise case, our baseline, and the $\xi$ mode in either the 1% or 5% Noise case. The errors for the $\xi^p$ modes are comparable to the errors for $\xi$.
<table><caption>TABLE 6. Relative Errors for ξ</caption><thead><tr><th>Mode</th><th>1% Noise Compared with 0% Noise</th><th>5% Noise Compared with 0% Noise</th></tr></thead><tbody><tr><td>ξ<sub>0</sub></td><td>1.22e-4 %</td><td>6.16e-4%</td></tr><tr><td>ξ<sub>1</sub></td><td>3.76e-3%</td><td>0.03%</td></tr><tr><td>ξ<sub>2</sub></td><td>9.29e-3%</td><td>0.25%</td></tr><tr><td>ξ<sub>3</sub></td><td>0.07%</td><td>1.18%</td></tr><tr><td>ξ<sub>4</sub></td><td>1.05%</td><td>1.57%</td></tr><tr><td>ξ<sub>5</sub></td><td>1.21%</td><td>2.29%</td></tr><tr><td>ξ<sub>6</sub></td><td>1.10%</td><td>4.11%</td></tr><tr><td>ξ<sub>7</sub></td><td>2.43%</td><td>13.00%</td></tr><tr><td>ξ<sub>8</sub></td><td>2.52%</td><td>16.47%</td></tr><tr><td>ξ<sub>9</sub></td><td>5.21%</td><td>203%</td></tr></tbody></table>