Data assimilation finite element method for the linearized Navier–Stokes equations in the low Reynolds regime

In this paper, we are interested in designing and analyzing a finite element data assimilation method for laminar steady flow described by the linearized incompressible Navier–Stokes equation. We propose a weakly consistent stabilized finite element method which reconstructs the whole fluid flow from noisy velocity measurements in a subset of the computational domain. Using the stability of the continuous problem in the form of a three balls inequality, we derive quantitative local error estimates for the velocity. Numerical simulations illustrate these convergence properties and we finally apply our method to the flow reconstruction in a blood vessel.


Introduction
The question of how to assimilate measured data into large scale computations of flow problems is receiving increasing attention from the computational community. There are several different situations where such data assimilation problems arise. One situation is when the data 4 Author to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. necessary to make the flow problem well posed is lacking, for instance when boundary data can not be obtained on parts of the boundary, but some other measured data on the boundary or in the bulk is available to make up for this shortfall. In such situation, the problem is ill-posed and numerical simulations are much more delicate to handle than for the wellposed flow equations. The traditional approach is to regularize the continuous problem to obtain a well-posed continuous problem, often using a variational framework, that can then be discretized using standard techniques. The regularization parameter then has to be tuned to have an optimal value with respect to noise in the data. The granularity of the computational mesh is chosen afterward to resolve all scales of the regularized problem. An example of this strategy is the quasi-reversibility method (see references [8][9][10][11]) which is applied to Stokes problem in [10] for the inverse identification of boundaries. Other examples can be found in [26], where additional measured data is used to compensate for a lack of knowledge of the boundary conditions in hemodynamics, or in [35], where a least squares method is proposed for combining and enhancing the results from an existing computational fluid dynamics model with experimental data. More generally the variational data assimilation method dates back to the work of Sasaki [38] where it was introduced in the context of meteorology. For further developments, we refer for example to [40,41]. Other approaches to data assimilation exist such as nudging [2,30,32] or minimax estimates [42]. The upshot in the present contribution is the design of a finite element method applied directly to the ill-posed variational data assimilation form. Regularization is then added on the discrete level, using methods from stabilized finite element methods allowing for a detailed analysis using conditional stability estimates.
A successful data assimilation method hinges on the existence of some stability properties of the ill-posed problem. Fortunately, it is known that a relatively large class of ill-posed problems has some conditional stability property. Stability estimates give a precise information on the effect of perturbations on the system. In particular, they imply that, if the data are compatible with the PDE, in the sense that there exists a solution in some suitable Sobolev space satisfying both the PDE and the data, then this solution is unique. For the Stokes equation, this unique continuation property was originally proven by Fabre and Lebeau [29]. The analysis of the stability properties of ill-posed problems based on the Navier-Stokes equations is a very active field of research and we refer to the works [3-5, 7, 33, 34, 36] for recent results. Stability estimates for inverse problems classically rely on Carleman inequalities or three-balls inequalities, two tools which are strongly related. The idea of applying Carleman estimates for the stability analysis of inverse problems is introduced in the seminal paper [16] by Bukhgeim and Klibanov.
There appears to be relatively few results in the literature discussing the combined error due to regularization, discretization and perturbations for inverse problems subject to the equations of fluid mechanics. To the best of our knowledge such a combined analysis has only been performed in the recent paper [22], where a nonconforming finite element method was used, together with regularization techniques developed in the context of discontinuous Galerkin methods, to analyze a data assimilation problem for Stokes problem. One of the reasons for this is that there is in general a gap between the stability estimates that can be proven analytically and the stability required to perform a numerical analysis. An approach allowing to bridge this gap was proposed in [19][20][21]23], drawing on earlier ideas for well-posed problems in [12,13]. This framework combines stabilized finite element methods designed for well-posed problems with variational formulations for data assimilation and sharp stability estimates for the continuous problems based on three balls inequalities or Carleman estimates. Recent developments include finite element data assimilation methods with optimal error estimates for the heat equation [24,25] and design of methods for indefinite or nonsymmetric scalar elliptic problems analyzed using Carleman estimates with explicit dependence on the physical parameters [17,18].
In this paper, our aim is to build on these results and use known techniques for the approximation of the (well-posed) Navier-Stokes equation in an optimization framework in order to assimilate data with computation. Contrary to the previous work [22], we consider the linearized Navier-Stokes equations and use standard H 1 -conforming, piecewise affine, finite element spaces. The key idea is that the ill-posed continuous problem is not regularized. Instead, we discretize the equation and set up a constrained optimization problem where we minimize the distance between the discrete solution and the measured data. To counter the instabilities in the discrete system, we introduce regularization terms. Taking advantage of the discretize-then-optimize perspective, we get an accuracy which is optimal with respect to the stability of the continuous problem through weakly consistent regularization terms. This property is illustrated by the error estimate of theorem 3.1 that is sharp relative to the stability estimate of corollary 2.1 (which is a consequence of a three balls inequality derived by Lin, Uhlmann and Wang [36]). Moreover, in our framework, the only regularization parameter, up to a constant scaling factor, is the mesh parameter.
We apply this method to an incompressible laminar steady flow in the low Reynolds regime. To be more specific, we are interested in a situation where a known laminar base flow U is available and we consider that a perturbation of the velocity of the base flow has been measured in some subset ω M of the computational domain Ω. Assuming that the perturbation is small, we then consider the linearized Navier-Stokes equation and our objective is to get quantitative local error estimates for the perturbation in some target subdomain ω T .
The rest of the paper is organized as follows. In section 2, we introduce the considered inverse problem and some related stability estimates. In section 3, we describe the proposed stabilized finite element approximation of the data assimilation problem and state the local error estimate. The numerical analysis of the method is carried out in section 4. Finally, section 5 presents a series of numerical examples which illustrate the performance of the proposed method. In particular, in section 5.3, we explore how the present approach can be applied to the estimation of relative pressure in blood flow from MRI velocity measurements.

Presentation of the inverse problem and stability results for the continuous problem
Let Ω be a bounded open polyhedral domain in R d with d = 2, 3. We denote by (U, P) a solution of the stationary incompressible Navier-Stokes equations and we consider some perturbation (u, p) of this base flow. It is then known that, if the quadratic term is neglected, the linearized Navier-Stokes equations for (u, p) may be written In what follows, we assume that U belongs to [W 1,∞ (Ω)] d and that (u, p) satisfies the regularity For this problem, we consider that measurements on u are available in some subdomain ω M ⊂ Ω having a nonempty interior and our objective is to reconstruct a fluid flow solution of system (1) from these measurements on the velocity.
Let us now introduce some useful notations. We consider the following spaces: where L 2 0 (Ω) = {p ∈ L 2 (Ω): Ω p = 0}. We also define the norms, for k = 1 or d, Let us notice that, in the first two definitions, with some abuse of notation, we use the same notation for k = 1 and k = d. For any subdomain ω ⊂ Ω, we set Besides, we introduce the bilinear forms: where H : The inverse problem we are interested in can be expressed in the following form: f ∈ V 0 being given, find (u, p) ∈ V × L 0 such that and Here, u M ∈ L 2 (ω M ) corresponds to the measurement of the velocity made on ω M and in what follows we will consider that this measurement corresponds to the exact velocity polluted by a small noise δu ∈ [L 2 (ω M )] d . For the analysis of our problem, we also introduce the linearized Navier-Stokes problem with a non-zero velocity divergence We assume that, if system (6) is completed by homogeneous Dirichlet boundary conditions, then it is well-posed. More precisely, we make the following assumption: Assumption A. For all f ∈ V 0 and g ∈ L 0 , we assume that system (6) admits a unique weak solution (u, p) ∈ V 0 × L 0 and that there exists a constant C S > 0 depending only on U, ν and Ω such that In particular, if ∇U [L ∞ (Ω)] d×d is small enough, it is straightforward to verify that Assumption A holds according to Lax-Milgram lemma. This assumption of smallness on ∇U however is a sufficient condition and there are reasons to believe that Assumption A holds in more general cases.
In the homogeneous case (which corresponds to f = 0 in (1) or to f = 0 and g = 0 in (6)), a solution (u, p) satisfies a three-balls inequality which only involves the L 2 norm of the velocity. This three-balls inequality result is stated in [36] (with their notations, A corresponds to U and B to ∇U ) and we recall it here. We emphasise the fact that, in this theorem and in its corollary, no boundary conditions are prescribed.
This theorem combined with assumption A implies a local stability inequality in the non-homogeneous case given by system (6). This stability property will be capital in the convergence study of the numerical method and is given by the following statement: Corollary 2.1 (Conditional stability for the linearized Navier-Stokes problem). Let f ∈ V 0 and g ∈ L 0 be given. For all ω T ⊂⊂ Ω, there exist C > 0 and 0 < τ < 1 such that The proof of corollary 2.1 is given in appendix A. In a classical way for ill-posed problems [1], corollary 2.1 gives a conditional stability result in the sense that, to be useful, this estimate has to be accompanied with an a priori bound on the solution on the global domain (due to the presence of u L in the right-hand side). Let us notice that corollary 2.1 implies in particular the uniqueness of a solution (u, p) in [H 1 (Ω)] d × H 1 (Ω) for problem (4) and (5), up to a constant for p. Moreover, in inequality (9), the exponent τ depends on the dimension d, the size of the measure domain ω M and the distance between the target domain ω T and the boundary of the computational domain Ω. Remark 2.1. Herein, we only consider error estimates for local L 2 -norms, but errors in local H 1 -norms of velocity together with L 2 -norms of pressure are also possible to analyze using three balls inequalities derived in [7], provided measurements of both velocities and pressures are available in ω M .
In what follows, we assume that f ∈ L 2 (Ω) and we introduce the operator A defined on where a and b are respectively defined by (2) and (3). Thus, we look for (u, p) ∈ V × L 0 such that and (4) holds.

Stabilized finite element approximation and error estimate
In this section, we first introduce a discretization of problem (11) using a standard finite element method. Then, the discrete inverse problem is reformulated as a constrained minimization problem in the discrete space where the regularization of the cost functional is achieved through stabilization terms. At last, the estimation of the error between the exact continuous solution and the discrete solution of our minimization problem is stated in theorem 3.1 which corresponds to our main theoretical result.
On the domain Ω, we consider a family {T h } h of shape regular, conforming, quasi-uniform meshes consisting of shape regular simplices K. This family is indexed by h defined as the maximum over the diameters h K of elements in the mesh. For a fixed h > 0, we denote by F i the set of interior faces of the mesh T h .
We will also use the notion of jump over a face F shared by the elements K and K (which means that F =K ∩K ) defined as follows: if ζ is a scalar, where n K denotes the outward pointing normal of the element K.
We denote by X h the standard H 1 -conforming finite element space of piecewise affine functions defined on T h and we introduce We may then write the finite element approximation of (11): To take into account the measurements on ω M given by (4), we introduce the measurement bilinear form where γ M > 0 will correspond to a free parameter representing the relative confidence in the measurements. The objective is then to minimize the functional under the constraint that (u h , p h ) satisfies (12). However, the discrete Lagrangian associated to this problem leads to an optimality system which is ill-posed. To regularize it, we introduce stabilization operators that will convexify the problem with respect to the direct variables u h , p h and the adjoint variables where γ u , γ p , γ * u and γ * p are positive user-defined parameters. Let us make some comments on these stabilization terms. The stabilization of the direct velocity acts on fluctuations of the discrete solution through a penalty on the jump of the solution gradient over element faces and has no equivalent on the continuous level. On the other hand, for the forward equation, the pressure stabilization s p coincides with a standard H 1 -Tikhonov regularization (let us notice that this stabilization term also appears the Brezzi-Pitkäranta method [14]). Since the exact solution of the adjoint equation is zero, consistency is ensured also for Tikhonov type regularization. For a more general discussion of the possible stabilization operators, we refer to [19,21].
For compactness, we introduce the primal and dual stabilizers: We may then write the discrete Lagrangian that will form the basis of our method as: If we differentiate with respect to (u h , p h ) and (z h , y h ), we get the following optimality system: for Let us prove that this discrete problem is well-posed. To do so, we recall a Poincaré inequality from [23, lemma 2] that will be crucial to get the stability of the method: If we take in the variational formulation (14) the test functions Then, according to (15), the left-hand side in equation (16) is the square of a norm in (V h × Q 0 h ) × (W h × Q h ) and, according to Babuska-Necas-Brezzi theorem (see [28]), we conclude that the square linear system defined by (14) admits a unique solution for all h > 0.
The following theorem is the main theoretical result of the paper and states an error estimate for this method.
The proof of this result will strongly rely on the conditional stability estimate for the continuous problem given by corollary 2.1 and the convergence of the residual quantities given by lemma 4.3.
According to this theorem, we get that, if the measurement noise δu is equal to 0, the error u − u h converges to 0 when h tends to 0 on any subset ω T ⊂⊂ Ω. Moreover, the convergence order τ corresponds to the Hölder coefficient of the continuous stability estimate (9). On the other hand, in the case of perturbed data, we notice that the accuracy of the error is limited and inequality (17) shows that the mesh size has to balance the error due to the discretization and the error due to the noise.

Remark 3.1.
For simplicity, we restrict the discussion to piecewise affine continuous approximation spaces, but the arguments can be extended to higher order finite element spaces, with the expected improvement of convergence order, following the ideas of [21]. It should however be noted that the system matrix becomes increasingly ill-posed as the polynomial order increases and the computation becomes more sensitive to noise in the measured data, so the practical interest in using high-order approximation spaces remains to be proven.

Stability and error analysis
In this section, we will present and prove several technical results and end with the proof of theorem 3.1.
Let us first notice that formulation (14) is weakly consistent in the sense that we have a modified Galerkin orthogonality relation with respect to the scalar product associated to A: Lemma 4.1 (Consistency). Let (u, p) satisfy (1) and (u h , p h ) be a solution of (14). Then there holds Proof. The result is immediate by taking the difference between (11) and the first equation of (14).
In the analysis below, we will use the following classical inverse and trace inequalities: • Inverse inequality (see [27,

section 1.4.3])
: Here P 1 (K) denotes the set of polynomials of degree less than or equal to 1 on the simplex K. • Trace inequalities (see [27, At last, by combining (20) and (19) v L 2 (∂K) Ch Let us now define the semi-norms associated to the stabilization operators defined on We also introduce the norm Let i h : L 2 (Ω) → X h be the Clément interpolant. We refer to [28] for the following results: for t = 1, 2, we have, for all u ∈ H t (Ω) and for all u ∈ H 2 (Ω) ⎛ Using the componentwise extension of i h to vectorial functions, we deduce the approximation bounds: The following continuity results for the bilinear form A motivate the definition of the triple norms. (24) and, for all

Lemma 4.2 (Continuity). For all ς ∈ V h + [H 2 (Ω)] d and ∈
Proof. The proof of (24) directly comes from the Cauchy-Schwarz inequality applied termwise in the definition (10) of A. For the second inequality (25), we setw = w − w h and y = y − y h and notice that For the first term in the right-hand side, an integration by parts in the viscous term gives, observing thatw| ∂Ω = 0, Using Cauchy-Schwarz inequality with the right scaling in h, we get Applying the trace inequality (20), we notice that Then, if we take the square, sum over K and use the H 1 −stability inequality of i h (22) for t = 1, this inequality becomes As a consequence, Similarly, for the second term in (26), an integration by parts gives Finally, the bound for the third term in (26) is immediate by the Cauchy-Schwarz inequality and the L 2 -stability of i h Gathering these results, we get (25).

Remark 4.1.
The augmented Lagrangian stabilization on the divergence in the operator s u is used in the proof of lemma 4.2 to bound in a direct way the third term in (26) but it is not strictly necessary. Indeed, if y h is chosen as the L 2 -projection of y we see that for all Hence, the stabilization of the gradient jump is sufficient to bound this term. In practice however, it can be useful to add the stabilization term on the divergence since it allows to get a stronger coercivity estimate.
The above lemma asserts that some residual converges with optimal order if the exact solution is smooth enough.

Lemma 4.3.
We assume that the solution (u, p) of (11) (14). Then there holds Using inequalities (23) and (22) for t = 1, we can directly bound the first two terms in the right-hand side. For the last two terms, according to inequality (15), we have To estimate the right-hand side, we notice that, using the second equation of (14) with Next, according to (18) Thus, adding these two equalities, we get We bound the terms I-III term by term. By lemma 4.2 and the approximation bound (23), we have for term I For term II, we have Thus, using (23) with (v, q) = (u, p) for the first term and the H 1 -stability of i h for the second term, we get For term III, according to (22) with t = 2, we have Thus, collecting the above bounds, we get and we conclude by dividing by (|||(ξ h , η h ) We then deduce from the previous lemma a priori bounds on the finite element solution. and Proof. In an evident way, we have Thus, according to lemma 4.3 and the fact that |||(u, p)||| = h u V + s p (p, p) 1 2 , we get (27). Moreover, by definition of ||| · |||, we have Thus inequality (27) directly implies (28).
We are now ready to prove our main result stated in theorem 3.1.
Proof of Theorem 3.1. Let us first introduce the weak formulation of the problem satisfied by (ξ, η) := (u − u h , p − p h ). By equation (11), we have for all w ∈ V 0 and q ∈ L, We introduce, u h and p h being fixed, the linear forms r f and r g on V 0 and L respectively defined by: for all w ∈ V 0 and q ∈ L, It follows that (ξ, η) is solution of (6) with the functions f and g in the right-hand sides replaced respectively by r f and r g . Applying now corollary 2.1, we directly get Using the first equation of (14), we can write the residuals: for all (w h , We take w h = i h w and q h = i h q in this equality. For the first term, according to (22) for t = 1, we have The second term can be bounded by using the relations (25) and (27). For the last term, we have, according to lemma 4.3 We thus get Since this bound holds for all w ∈ V 0 and q ∈ L, we conclude that Thus, we can bound the terms in the right-hand side of (29) in the following way: according to inequalities (27) and (28) and according to lemma 4.3. Using these two bounds in (29), we conclude that which completes the proof of theorem 3.1.

Numerical simulations
In this section, we apply the method introduced in section 3 in different two-dimensional different numerical examples. The free parameters in (14) are set to in all the numerical examples. The numerical computations have been performed with FreeFEM++ (see [31]).

Convergence study: Stokes example
In order to illustrate the convergence behavior of the method introduced in section 3, we take up the test case presented in [22] for the Stokes system. In the unit square Ω = (0, 1) 2 , we consider the velocity and pressure fields given by It is straightforward to verify that (u, p) is a solution to the homogeneous Stokes problem with ν = 1, which corresponds to system (1) with U = 0 and f = 0. We hence consider the formulation (14)  First, we perform the computation with unperturbed data. In figure 1 (left), we report the velocity and pressure errors both in the global L 2 -norm and the local L 2 -norm in the subdomain ω T . We also provide the convergence history of the residual quantity for the velocity stabilization: ⎛ The observed global asymptotic behaviors of the local velocity error (filled squares) and residual (filled circles) are in agreement with the convergence rates obtained in theorem 3.1 and corollary 4.1 with |δu| ω M = 0. It should be noted that, for the finest grids, the local velocity error tends to stagnate or increase, which can be related either to the impact of the roundingoff errors or to ill-conditioning issues of the system matrix, so that |δu| ω M > 0. The other error quantities, global velocity error (empty squares) and local and global pressure errors (filled and empty triangles, respectively) show a convergent behavior which also tends to stagnate for the smallest values of h. Figure 1 (right) presents the convergence history of the same quantities with a 10% Gaussian noise. The impact of the noise is clearly visible. In particular, it is worth noting that the convergence history of the local and global velocity and pressure errors is not monotone anymore and the residual loses first-order convergence rate. This is also in agreement with theorem 3.1 and corollary 4.1 with |δu| ω M > 0.
As in the previous case, we have performed numerical tests for unperturbed data and for data perturbed with a 10% noise and we have studied the convergence of the method. The obtained results are illustrated in figure 2. The convergence curves present similarities with the a contraction of 60%. The radii of inlet and outlet are 1 cm, the length of the vessel is 6 cm and the dynamic viscosity is given by ν = 0.035 Poise. Synthetic measurements are first generated by numerically solving the incompressible Navier-Stokes system with the following boundary and initial conditions: This direct problem (31) and (32) is discretized in space by continuous piece-wise affine finite element approximations based on the SUPG/PSPG stabilization method. The time discretization consists in a backward Euler scheme with a semi-implicit treatment of the convective term. A standard backflow stabilization term is also applied on the outlet boundary Γ o in order to guarantee the overall stability of the numerical scheme (see [15]). The discretization parameters are set to Δt = 0.002 s and h = 0.01 cm. This space-time grid generates a set of synthetic velocity measurements which can be perturbed either by noise or by space-time subsampling. Figure 4 represents the estimate RPD with the same discretization parameters as for the direct problem (no subsampling). When the data are unperturbed, we see that the reconstructed curve is perfectly superimposed with the exact curve ( figure 4, left). With a 10% Gaussian noise, we observe that the reconstructed curve (which corresponds to the mean curve obtained from 30 tests with variable noises) succeeds in following accurately the variations of the reference data ( figure 4, right).
In figure 5, the measurements are perturbed by a subsampling both in time and in space (the time step is 10 times larger and the mesh size is 8 times larger). We then solve the data assimilation problem with the time step or the mesh size corresponding to this subsampling.   Figure 5 shows that the proposed approach is able to provide a reasonable estimation of the RPD with and without noise (10% Gaussian noise). In particular, we can clearly observe that the RPD peak is well captured in both cases. Moreover, we can see in figure 6 that the velocity is well-reconstructed in the two-dimensional domain.
Here and in the sequel, we will frequently use the notation a b for a Cb for some C > 0. For the second term in the right-hand side of (33), applying theorem 2.1, we get that (w,ỹ) satisfies (8) and thus: We now revert back to u in the right-hand side: We conclude the proof by collecting the estimates forũ andw.
If ω M and ω T do not satisfy the assumptions for the construction of the balls B R 1 (x 0 ) and B R 0 (x 0 ), we introduce a finite sequence of intermediate balls in order to link ω T to ω M (as it is done for instance in [37]) and we get again the estimate.