In this work, we consider a data assimilation problem for a stationary convection diffusion equation

when convection dominates, that is $0<\mu\ll|\beta|$, and complement the diffusion dominated case discussed in the first part [BNO20]. We assume that $\Omega\subset\mathbb{R}^{n}$ is open, bounded and connected, and there exists a solution $u\in H^{2}(\Omega)$ to (1). The problem under study is to approximate the solution $u$ given the source $f$ in $\Omega$ and the perturbed restriction $\tilde{U}_{\omega}=\left.u\right|_{\omega}+\delta$ of the solution to an open subset $\omega\subset\Omega$. The perturbation $\delta$ is taken in $L^{2}(\omega)$. Notice that we consider no boundary conditions on $\partial\Omega$. This is a linear ill-posed problem also known as unique continuation.

To start with, let us briefly recall the main results obtained in the first part. Consider an open bounded set $B\subset\Omega$ that contains the data region $\omega$ such that $B\backslash\omega$ does not touch the boundary of $\Omega$. For $u\in H^{1}(\Omega)$, the following conditional stability estimate was proven for $\mu>0$ and $\beta\in L^{\infty}(\Omega)^{n}$,

where the Hölder exponent $\kappa\in(0,1)$ depends only on the geometric setting. In the case of simple geometric configurations, e.g. when $\omega,B,\Omega$ are three concentric balls, the exponent $\kappa\in(0,1)$ can be given explicitly, see [BNO20, Remark 1]. The stability constant $C_{st}$ is given explicitly in terms of the physical parameters

with constants $C_{1},C_{2}>0$ depending only on the geometry. Note that the continuum estimate (2) is valid in both the diffusion-dominated and convection-dominated regimes, and that the stability constant $C_{st}$ is uniformly bounded when diffusion dominates. However, when convection dominates $C_{st}$ grows exponentially, rendering the stability estimate ineffective in practice. We also recall that for global unique continuation from $\omega$ to the entire $\Omega$ the stability is no longer Hölder, but logarithmic, that is the modulus of continuity for the given data is not $|\cdot|^{\kappa}$ any more, but $|\log(\cdot)|^{-\kappa}$.

On the discrete level, the continuum estimate (2) was combined with a stabilized linear finite element method to obtain convergence orders for the approximate solution. More precisely, for a mesh size $h$, and defining the Péclet number

the following error bound [BNO20, Theorem 1] was proven for the approximation $u_{h}$ in the diffusive regime $Pe(h)<1$,

where the convergence order $\kappa\in(0,1)$ is the same as the Hölder exponent in (2) and the stability constant $C_{st}$ is proportional to the one in (3). Under an additional assumption on the divergence of the convective field $\beta$, similar error bounds were also proven in the $H^{1}$-norm, see [BNO20, Theorem 2].

The prototypical effect of dominating diffusion is shown in Figure 1, where the problem is set in the unit square and contour error plots are shown for data assimilation from a centered disk of radius 0.1 . One can notice oscillating errors that grow in size away from the data region towards the boundary. The exact solution in this example is $u=2\sin(5\pi x)\sin(5\pi y)$ where the factor 2 is taken to make its $L^{2}$-norm unitary. For the computation we used an unstructured mesh with 512 elements on a side and mesh size $h\approx 0.0025$.
1.1. Objective and main results. We consider a stabilized finite element method for data assimilation subject to the convection-diffusion equation in the convection dominated regime. Since the behaviour of the physical system changes fundamentally when convection dominates and

the goal of this second part paper is to reconsider the numerical method proposed in the first part [BNO20] and develop an error analysis that captures and exploits the governing transport phenomenon. This is illustrated in Figure 2 where the transition to the convection-dominated regime through an intermediate regime is made by decreasing the diffusion coefficient $\mu$. We aim to obtain sharper local error estimates along the characteristics of the convective field through the data region. Even though the error analysis that we perform herein is different in nature to the one in the first part, the numerical method itself is only slightly changed (see Remark 1 below). For the error localization technique we draw on ideas used for the streamline diffusion method in [JNP84], continuous interior penalty in [BGL09], and non-coercive hyperbolic problems in [Bur14].

From the definition of the Péclet number we see that the regime will also depend on the resolution of the computation besides the physical parameters. We can therefore expect the method to change behaviour as the resolution increases and $Pe(h)$ decreases. This phenomenon was already observed computationally in [Bur13] and can now be explained theoretically.

To make the presentation as simple as possible we consider a model case in the unit square $\Omega$ with constant convection

and the solution given in the subset

For the subset $\omega_{\beta}\subset\Omega$ covered by the characteristics of $\beta$ that go through $\omega$, we introduce the stability region $\stackrel_{\beta}\subset\omega_{\beta}$ by cutting off a crosswind layer of width $\mathcal{O}\left(h^{\frac{1}{2}}|\ln(h)|\right)$ (see Section 2.1 for more details). We separate the convection-dominated and diffusion dominated regimes by introducing a constant $Pe_{\text{lim }}>1$ such that

To reduce the number of constants appearing in the analysis, we will write this as $Pe(h)\gtrsim 1$. As suggested by Figure 2, we expect different results for data assimilation downstream vs upstream in an intermediate regime. We prove in Theorem 1 weighted error estimates that for $\beta_{1}>0$ essentially take the following form

This is similar to the typical error estimates for piecewise linear stabilized FEMs for convection-dominated well-posed problems, such as local projection stabilization, dG methods, continuous interior penalty or Galerkin least squares. On general shape-regular meshes these methods have an $\mathcal{O}\left(h^{\frac{1}{2}}\right)$ gap to the best approximation convergence order. Taking this into account, our result is thus quasi-optimal. For a recent overview of challenges and open problems in the well-posed case, see e.g. [JKN18] and [RS15].

When going against the characteristics, i.e. $\beta_{1}<0$, we prove in Theorem 2 first the pre-asymptotic bound

followed by

It follows that when solving the data assimilation problem against the flow, the diffusivity reduces the convergence order in an intermediate regime. Only for very small diffusion coefficients $\mu<|\beta|h^{2}$ do we get quasi-optimal bounds. This asymmetry of the error distribution for moderate Péclet numbers is clearly visible in the left plot of Figure 2.

Previous results on optimal control for stabilized convection-diffusion equations include [BV07, DQ05, HYZ09, YZ09]. We refer to the first part [BNO20] for a more detailed discussion.

In terms of notation, above and throughout the paper $C$ denotes a generic positive constant, not necessarily the same at each occurrence, that is independent of the coefficients $\mu,\beta$ and the mesh size $h$.

Let $V_{h}\subset H^{1}(\Omega)$ be the conforming finite element space of piecewise affine $\mathbb{P}_{1}$ functions defined on a computational mesh $\mathcal{T}_{h}$ that consists of shape-regular triangular elements $K$ with diameter $h_{K}$. The mesh size $h$ is the maximum over $h_{K}$ and we will assume that $h<1$. The interior faces of all the elements are collected in the set $\mathcal{F}_{i}$ and the jump of a quantity across such a face $F$ is denoted by $\llbracket\cdot\rrbracket_{F}$, omitting the subscript whenever the context is clear. We denote by $n$ the unit normal.

First we introduce the standard inner products with the induced norms

and the bilinear form in the weak formulation of (1)

We will make use of the stabilizing bilinear forms

which will act on the discrete solution penalizing the jumps of its gradient across interior faces, and

where $\gamma$ and $\gamma_{*}$ are positive constants that can be heuristically chosen at implementation. They do not play a role in the convergence of the method and most of the time we will include them in the generic constant $C$. For the data assimilation term we consider the scaled inner product in the data set $\omega$ given by

To this we add the stabilizing interior penalty term $s_{\Omega}$ to define for conciseness

The idea behind the computational method follows the discretize-then-optimize approach: we first discretize and then formulate the data assimilation problem as a PDEconstrained optimization problem with additional stabilizing terms. Apart from their stabilizing intake, these terms are also chosen such that they vanish at optimal rates. For an overview on this approach to ill-posed problems and more details on the desired properties of discrete stabilizers, we refer the reader to [Bur13]. To be more precise, for an approximation $u_{h}\in V_{h}$ and a discrete Lagrange multiplier $z_{h}\in V_{h}$, we consider the functional

where the first term is measuring the misfit between $u_{h}$ and the known perturbed restriction $\tilde{U}_{\omega}=\left.u\right|_{\omega}+\delta$, the next two terms are imposing the weak form of the PDE (1) as a constraint, and the last two terms have stabilizing role and act only on the discrete level.

We look for the saddle points of the Lagrangian $L_{h}$ and analyse their convergence to the exact solution. Using the optimality conditions we obtain the finite element method for data assimilation subject to (1), which reads as follows: find $\left(u_{h},z_{h}\right)\in\left[V_{h}\right]^{2}$ such that

Notice that the exact solution $u\in H^{2}(\Omega)$ (with noise $\delta\equiv 0$ ) and the dual variable $z\equiv 0$ satisfy (5) since the gradient of the exact solution has no jumps across interior faces. Hence the Lagrange multiplier $z_{h}$ should converge to zero.

Remark 1. The same finite element method (5) has been proposed in the first part [BNO20] for the diffusion-dominated case; $s_{\Omega}$ and $s_{*}$ are exactly the stabilizing terms introduced there. However, herein we have increased the penalty coefficient in the data term $s_{\omega}$ from $|\beta|h+\mu$ to $|\beta|h^{-1}+\mu h^{-\zeta}$. We note that the bounds in [BNO20] hold also for this stronger penalty term, but the sensitivity to perturbations in data increases by a factor of $h^{-1}$.

Proposition 1. The finite element method (5) has a unique solution $\left(u_{h},z_{h}\right)\in\left[V_{h}\right]^{2}$ and the Euclidean condition number $\mathcal{K}_{2}$ of the system matrix satisfies

Proof. The proof given in [BNO20, Proposition 2] holds verbatim. $\square$
2.1. Stability region and weight functions. We will exploit the convective term of the PDE to obtain stability in the zone that connects through characteristics to the data region $\omega$. The objective is to obtain weighted $L^{2}$-estimates in this region that are independent of $\mu$ (but not of the regularity of the exact solution) with the underlying assumption that $\mu\ll|\beta|$. To this end we first define the subdomain where we can obtain stability (see Figure 3 for a sketch) and some weight functions that will be used to define weighted norms. These can be given in explicit form in the simple framework where $\beta=\left(\beta_{1},0\right)$ and

Let $\omega_{\beta}$ denote the closed set of all the points $p\in\bar{\Omega}$ that can be reached through characteristics from $\omega$, i.e. for which there exists $s\in\mathbb{R}$ such that $p+s\beta\in\partial\omega$. Similarly to the classical work [JNP84], we define the stability region $\dot{\omega}_{\beta}$ by cutting off a crosswind layer from $\omega_{\beta}$, namely

with the constant $c_{\lambda}$ to be made precise in the following. In our setting, we simply have that $\stackrel_{\beta}=[0,1]\times\left[\stackrel^{-},\stackrel^{+}\right]$for some $\stackrel^{+}>\stackrel^{-}>0$. The crosswind layer and its width are motivated by the subsequent construction of weight functions with a specific decay outside $\omega_{\beta}$.

We will consider different weight functions depending on the direction of the convection field. In the downstream case we let

and in the upstream case

In both cases we have that $\nabla\psi_{1}=\left(-\psi_{1},0\right)$. Let then $\psi_{2}\in W^{1,\infty}(\Omega)$ be a function satisfying

Such a function can be obtained by taking a positive constant $\lambda$ that will be specified later and letting

Note that $\psi_{2}$ is only piecewise continuously differentiable. For $\psi_{2}$ to decrease sufficiently rapidly to $\mathcal{O}\left(h^{3}\right)$ outside of $\omega_{\beta}$, we can take

which corresponds to $c_{\lambda}=3\lambda$ in the definition of $\dot{\omega}_{\beta}$ given in (6). We thus have that

and in the subsequent proofs the constant $\lambda$ will be taken large enough.
We now define the weight function $\varphi\in W^{1,\infty}(\Omega)$ that will be used in the weighted norms. For the downstream case we take in Section 4.3

and for the upstream case in Section 4.4,

Using the product rule and the fact that $\beta\cdot\nabla\psi_{2}=0$, it follows that in both cases we have

and

We will denote the inflow and outflow boundaries by $\partial\Omega^{-}$and $\partial\Omega^{+}$, i.e. $\beta\cdot n<0$ on $\partial\Omega^{-}$and $\beta\cdot n>0$ on $\partial\Omega^{+}$. We will also assume that $\beta\cdot n=0$ can only hold on the boundary of $\Omega\backslash\omega_{\beta}$, and that $\mu\leq\mathrm{Pe}_{\text{lim }}^{-1}|\beta\cdot n|h$ when $\beta\cdot n\neq 0$. This is straightforward to verify in the model case of the unit square that we are considering.

We first collect several inequalities that will be used repeatedly. We recall the standard discrete inverse inequality

see e.g. [EG04, Lemma 1.138], the continuous trace inequality

see e.g. [MS99], and the discrete trace inequality

We will use standard estimates for the $L^{2}$-projection $\pi_{h}:L^{2}(\Omega)\mapsto V_{h}$, namely

Scaling the result in [BHL18, Lemma 2] we recall the Poincaré-type inequality

Using (13) and approximation estimates we also have the jump inequality

We also recall that for a Lipschitz domain $K$ - and hence for any element $K\in\mathcal{T}_{h}$ - a function $\varphi$ is Lipschitz continuous if and only if $\varphi\in W^{1,\infty}(K)$. This follows from the proof in [Eva10, Theorem 4, p. 294] where the extension operator in the third step of the proof is replaced by the extension operator in [Ste70, Theorem 5, p. 181]. This equivalence holds for more general domains satisfying the minimal smoothness property in [Ste70, p. 189]. The proof in [Eva10, Theorem 4, p. 294] also shows that the mean value theorem holds and for any $x,y\in K$,

where $|\varphi|_{W^{1,\infty}(K)}:=\|\nabla\varphi\|_{\infty,K}$ and the constant $C_{ex}>0$ is the norm of the extension operator used.

Lemma 1. For all $v_{h}\in V_{h}$ and $K\in\mathcal{T}_{h}$, the following inequalities hold

assuming that ( $h+\lambda^{-1}h^{\frac{1}{2}}$ ) is small enough.
Proof. Let $x^{*}\in K$ be such that $\left|\varphi\left(x^{*}\right)\right|=\|\varphi\|_{\infty,K}$. Using the triangle inequality we may write

By the mean value theorem (17) we have that

and by (11) together with the assumption that $C_{ex}\left(h+\lambda^{-1}h^{\frac{1}{2}}\right)<\frac{1}{2}$ we get

It follows that

from which the claim (18) is immediate. Considering now (19), using the standard element-wise trace inequality (13) we have

We bound the gradient term using (11) and the inverse inequality (12),

We conclude by collecting the terms and using (18).
3.1. Discrete commutator property. We denote by $i_{h}$ the Lagrange nodal interpolant. We herein consider a Lipschitz weight function and prove the following super approximation result, also known as the discrete commutator property. This result will be essential to derive local weighted estimates and is similar to the one proven in [Ber99] for smooth compactly supported weight functions. For an introduction to interior estimates we refer the reader to [NS74].
Lemma 2. Let $v_{h}\in V_{h}$ and $K\in\mathcal{T}_{h}$. Then for any weight function $\varphi\in W^{1,\infty}(K)$

Proof. We will first show the $L^{2}$-norm estimate

Let $x^{*}\in K$ be such that $\left|\varphi\left(x^{*}\right)\right|=\|\varphi\|_{\infty,K}$ and let $R_{\varphi}=\varphi-\varphi\left(x^{*}\right)$. Note that

Observe that $i_{h}\left(v_{h}\varphi\right)=i_{h}\left(v_{h}i_{h}\varphi\right)$ and therefore

By the triangle inequality

For the first term, since $v_{h}i_{h}R_{\varphi}\in H^{1}(K)$ we have by standard interpolation that

and then

By inserting $\nabla\varphi$ and $\varphi$, respectively, and using interpolation estimates in $W^{1,\infty}(K)$ [EG04, Theorem 1.103] and the mean value theorem (17) we have the following approximation

Combined with the previous estimate and the inverse inequality (12) this gives that

For the second term, using again interpolation [EG04, Theorem 1.103] we have

and thus we have shown that

The approximation estimate for the gradient follows by the same arguments. Indeed,

We first use interpolation and the inverse inequality (12) to get

Then we use an inverse inequality followed by interpolation and (20) to obtain

4.1. Consistency and continuity. The following consistency result holds exactly as in the diffusion-dominated case, see [BNO20, Lemma 4]. We give the proof for the sake of completeness.
Lemma 3 (Consistency). Let $u\in H^{2}(\Omega)$ be the solution to (1) and $\left(u_{h},z_{h}\right)\in\left[V_{h}\right]^{2}$ the solution to (5), then

and

for all $\left(v_{h},w_{h}\right)\in\left[V_{h}\right]^{2}$.
Proof. The first claim follows from the definition of $a_{h}$, since

where in the last equality we integrated by parts. The second claim follows similarly from

which combined with the fact that $s_{\Omega}\left(u,v_{h}\right)=0$ leads to

We now introduce the stabilization norm on $\left[V_{h}\right]^{2}$ by combining the primal and dual stabilizers

and the "continuity norm" defined on $H^{\frac{3}{2}+\varepsilon}(\Omega)$, for any $\varepsilon>0$,

From the jump inequality (16), standard approximation bounds for $\pi_{h}$ and the trace inequality (13), it follows that for $u\in H^{2}(\Omega)$