A stabilized finite element method for inverse problems subject to the convection–diffusion equation. II: convection-dominated regime

We consider the numerical approximation of the ill-posed data assimilation problem for stationary convection–diffusion equations and extend our previous analysis in Burman et al. (Numer. Math. 144:451–477, 2020) to the convection-dominated regime. Slightly adjusting the stabilized finite element method proposed for dominant diffusion, we draw upon a local error analysis to obtain quasi-optimal convergence along the characteristics of the convective field through the data set. The weight function multiplying the discrete solution is taken to be Lipschitz continuous and a corresponding super approximation result (discrete commutator property) is proven. The effect of data perturbations is included in the analysis and we conclude the paper with some numerical experiments.


Introduction
In this work, we consider a data assimilation problem for a stationary convectiondiffusion equation ( 1) when convection dominates, that is 0 < µ |β|, and complement the diffusiondominated case discussed in the first part [BNO20].We assume that Ω ⊂ R n is open, bounded and connected, and there exists a solution u ∈ H 2 (Ω) to (1).The problem under study is to approximate the solution u given the source f in Ω and the perturbed restriction Ũω = u| ω + δ of the solution to an open subset ω ⊂ Ω.The perturbation δ is taken in L 2 (ω).Notice that we consider no boundary conditions on ∂Ω.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 ⊂ Ω that contains the data region ω such that B \ ω does not touch the boundary of Ω.For u ∈ H 1 (Ω), the following conditional stability estimate was proven for µ > 0 and β ∈ L ∞ (Ω) n , (2) , where the Hölder exponent κ ∈ (0, 1) depends only on the geometric setting.In the case of simple geometric configurations, e.g. when ω, B, Ω are three concentric balls, the exponent κ ∈ (0, 1) can be given explicitly, see [BNO20, Remark 1].The stability constant C st is given explicitly in terms of the physical parameters (3) 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 ω to the entire Ω the stability is no longer Hölder, but logarithmic, that is the modulus of continuity for the given data is not | • | κ any more, but | log(•)| −κ .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 P e(h) := |β|h µ , the following error bound [BNO20, Theorem 1] was proven for the approximation u h in the diffusive regime P e(h) < 1, where the convergence order κ ∈ (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 β, similar error bounds were also proven in the 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πx) sin(5π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 ≈ 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 convectiondominated regime.Since the behaviour of the physical system changes fundamentally when convection dominates and P e(h) 1, 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 µ.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 P e(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 Ω with constant convection and the solution given in the subset ω := (0, x) × (y − , y + ) with x > h and y + − y − > h.
For the subset ω β ⊂ Ω covered by the characteristics of β that go through ω, we introduce the stability region ωβ ⊂ ω β by cutting off a crosswind layer of width O(h 1 2 | ln(h)|) (see Section 2.1 for more details).We separate the convection-dominated and diffusiondominated regimes by introducing a constant P e lim > 1 such that P e(h) > P e lim > 1.
To reduce the number of constants appearing in the analysis, we will write this as P e(h) 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 β 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 O(h 1 2 ) 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. β 1 < 0, we prove in Theorem 2 first the pre-asymptotic bound 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 µ < |β|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 µ, β and the mesh size h.

Discrete setting
Let V h ⊂ H 1 (Ω) be the conforming finite element space of piecewise affine P 1 functions defined on a computational mesh 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 F i and the jump of a quantity across such a face F is denoted by • 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 γ and γ * 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 ω given by To this we add the stabilizing interior penalty term s Ω 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 ∈ V h and a discrete Lagrange multiplier z h ∈ V h , we consider the functional where the first term is measuring the misfit between u h and the known perturbed restriction Ũω = u| ω + δ, 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 Notice that the exact solution u ∈ H 2 (Ω) (with noise δ ≡ 0) and the dual variable z ≡ 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 Ω and s * are exactly the stabilizing terms introduced there.However, herein we have increased the penalty coefficient in the data term s ω from |β|h + µ to |β|h −1 + µh −ζ .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 . Ω Figure 3. Data set ω (gray) and the stability region ωβ (hatched).
Proposition 1.The finite element method (5) has a unique solution and the Euclidean condition number K 2 of the system matrix satisfies Proof.The proof given in [BNO20, Proposition 2] holds verbatim.
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 ω.The objective is to obtain weighted L 2 -estimates in this region that are independent of µ (but not of the regularity of the exact solution) with the underlying assumption that µ |β|.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 β = (β 1 , 0) and ω := (0, x) × (y − , y + ) with x > h and y + − y − > h.
Let ω β denote the closed set of all the points p ∈ Ω that can be reached through characteristics from ω, i.e. for which there exists s ∈ R such that p + sβ ∈ ∂ω.Similarly to the classical work [JNP84], we define the stability region ωβ by cutting off a crosswind layer from ω β , namely , with the constant c λ to be made precise in the following.In our setting, we simply have that ωβ = [0, 1] × [ẙ − , ẙ+ ] for some ẙ+ > ẙ− > 0. The crosswind layer and its width are motivated by the subsequent construction of weight functions with a specific decay outside ω β .
We will consider different weight functions depending on the direction of the convection field.In the downstream case we let ψ 1 (x, y) := e −x , when β 1 > 0, and in the upstream case ψ 1 (x, y) := −e −x , when β 1 < 0.
In both cases we have that ∇ψ 1 = (−ψ 1 , 0).Let then ψ 2 ∈ W 1,∞ (Ω) be a function satisfying (7) Such a function can be obtained by taking a positive constant λ that will be specified later and letting Note that ψ 2 is only piecewise continuously differentiable.For ψ 2 to decrease sufficiently rapidly to O(h 3 ) outside of ω β , we can take which corresponds to c λ = 3λ in the definition of ωβ given in (6).We thus have that and in the subsequent proofs the constant λ will be taken large enough.We now define the weight function ϕ ∈ W 1,∞ (Ω) 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, (9) Using the product rule and the fact that β • ∇ψ 2 = 0, it follows that in both cases we have We will denote the inflow and outflow boundaries by ∂Ω − and ∂Ω + , i.e. β • n < 0 on ∂Ω − and β • n > 0 on ∂Ω + .We will also assume that β • n = 0 can only hold on the boundary of Ω \ ω β , and that µ ≤ Pe −1 lim |β • n|h when β • n = 0.This is straightforward to verify in the model case of the unit square that we are considering.

Preliminaries and the discrete commutator property
We first collect several inequalities that will be used repeatedly.We recall the standard discrete inverse inequality see e.g.[MS99], and the discrete trace inequality We will use standard estimates for the L 2 -projection Scaling the result in [BHL18, Lemma 2] we recall the Poincaré-type inequality (15) (µ 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 ∈ T h -a function ϕ is Lipschitz continuous if and only if ϕ ∈ W 1,∞ (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 ∈ K, where |ϕ| W 1,∞ (K) := ∇ϕ ∞,K and the constant C ex > 0 is the norm of the extension operator used.
Lemma 1.For all v h ∈ V h and K ∈ T h , the following inequalities hold Using the triangle inequality we may write and by (11) together with the assumption 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].
Proof.We will first show the L 2 -norm estimate By the triangle inequality For the first term, since and then By inserting ∇ϕ and ϕ, respectively, and using interpolation estimates in W 1,∞ (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 (20) 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

A priori local error estimates
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 ∈ H 2 (Ω) be the solution to (1) and (u h , z h ) ∈ [V h ] 2 the solution to (5), then and 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 We now introduce the stabilization norm on [V h ] 2 by combining the primal and dual stabilizers (v h , w h ) 2 s := s(v h , v h ) + s * (w h , w h ), and the "continuity norm" defined on H 3 2 +ε (Ω), for any ε > 0, From the jump inequality (16), standard approximation bounds for π h and the trace inequality (13), it follows that for u ∈ H 2 (Ω) . We also define the orthogonal space Proof.Integrating by parts in the convective term of a h and using ∇ • β = 0 leads to For the first term we recall the discrete approximation estimate that holds for any piecewise linear β, see e.g.[Bur05, Theorem 2.2], inf and use orthogonality to obtain For the remaining terms, applying the Cauchy-Schwarz inequality we see that Note that the proof of the above continuity estimate holds for any divergence-free piecewise linear velocity field β.To address the case of a general velocity field β ∈ W 1,∞ (Ω) one can use a similar argument by considering its piecewise linear approximation as in [BNO20, Lemma 5].Assuming that β is divergence-free, the constant would be proportional to hP e(h) otherwise it would be proportional to P e(h) 4.2.Convergence of regularization.We now prove optimal convergence for the stabilizing and data assimilation terms.
Proposition 2. (Convergence of regularization).Let u ∈ H 2 (Ω) be the solution of (1) and (u h , z h ) ∈ [V h ] 2 the solution to (5), then there holds Using both claims in Lemma 3 we may write We bound the other terms using the Cauchy-Schwarz inequality

Collecting the above bounds we have
1 2 δ ω ) (e h , z h ) s and the claim follows by applying the approximation inequality (21).
Remark 2. Compared to the result in the diffusion-dominated case [BNO20, Proposition 3], the sensitivity to data perturbations has increased by a factor of h −1 .This is due to the stronger penalty in the data term s ω (c.f.Remark 1).

Downstream estimates.
In this case we consider β = (β 1 , 0) with β 1 > 0 and the data set ω = (0, x) × (y − , y + ) touching part of the inflow boundary ∂Ω − .We aim to obtain control of the following weighted triple norm defined on V h , (23) ∂Ω + , where ϕ is given by (8).Since ϕ ∈ (0, 1), we will often use that We first consider v h ϕ as a test function in the weak bilinear form a h and obtain the following bound.
Lemma 5.There exist α > 0 and h 0 > 0 such that for all h < h 0 and v h ∈ V h we have Proof.We start with the convective term.Since ∇ • β = 0, the divergence theorem leads to 2(β Then combining with (10) we have Ω .We split the boundary term into inflow and outflow ∂Ω − .Splitting now the inflow boundary with respect to the closed set ω β and using the discrete trace inequality (13) in ω, and the weight decay (7) together with a standard global trace inequality for H 1 functions outside, we have that where in the last step we used the Poincaré-type inequality (15).Hence we obtain control over the convective terms in the triple weighted norm Let us consider the terms in a h (v h , v h ϕ) corresponding to the diffusion operator, starting with Ω + (µ∇v h , v h ∇ϕ) Ω , which we rearrange as Bounding ∇ϕ by (11) and using Cauchy-Schwarz together with µ ≤ |β|h we have that We split the boundary term µ∇v h • n, v h ϕ ∂Ω into inflow and outflow again.Similarly to (24) we have that For the outflow boundary term we use Cauchy-Schwarz and a trace inequality to obtain We denote the part of the boundary where β • n = 0 by ∂Ω 0 and use the assumption that ∂Ω 0 is away from ω β , meaning that the weight function ϕ is O(h 3 ) there.We use trace inequalities and (15) to bound Collecting the above bounds we obtain that We conclude by combining this with (25) and assuming that h is small enough and P e lim are large enough (thus absorbing the convective terms from the rhs into the lhs).Now we refine the control over the triple norm |||v h ||| ϕ by taking the projection π h (v h ϕ) ∈ V h as a test function.
Corollary 1.There exist α > 0 and h 0 > 0 such that for all h < h 0 and v h ∈ V h we have we must bound a h (v h , (π h − 1)(v h ϕ)) in a suitable way.Writing out the terms we have Let us consider the convection term first, and use orthogonality combined with (22) Integrating by parts and using that ∆v h = 0 on every element K we obtain by the trace inequality (13) and the assumption P e(h) > 1 that ) and the stability of the projection gives and hence ) Ω , collecting the contributions above we see that and hence The discrete commutator property Lemma 2 together with the ϕ-bounds (11) and (18) give that (28) Since ϕ ∈ (0, 1) and ϕ < ϕ 1 2 , it follows that for h small enough and λ large enough, given some α > 0 we have Collecting the estimates for a h (v h , (π h − 1)(v h ϕ)) and using Lemma 5 we see that from which we conclude by renaming α/2 as α.
Lemma 6.For all v h ∈ V h there holds Proof.First note that by triangle inequalities we have that up to a constant We bound these terms line by line.Using ( 27), ( 28), ( 11) and ( 18) we bound the first line by For the second line, using a global trace inequality and the stability of the projection we have Splitting the boundary into inflow and outflow and using (24), we have For the contribution of the jump term, we insert i h and bound We first observe that using ( 14) and (26), we can bound the first term by where for the last two inequalities we used the discrete commutator property Lemma 2 together with the ϕ-bounds (11) and (18).Since ϕ is Lipschitz continuous on K, ϕ| F is also Lipschitz continuous, and so ϕ| F ∈ W 1,∞ (F ).The restriction of the nodal interpolant on K onto F gives the nodal interpolant on F , hence applying Lemma 2 to F instead of K we have the discrete commutator estimate where in the last step we used (11) and (18) together with a discrete trace inequality.After summation we have that Finally we use the trivial bound (since |ϕ| < 1) We conclude the proof by summing up the above contributions.
We can now prove the following error estimate showing that, in the zone ωβ where we have stability, the convergence in the L 2 -norm is of order O(h 3 2 ) on unstructured meshes, which is known to be optimal.
Theorem 1.Let u ∈ H 2 (Ω) be the solution of (1) and (u h , z h ) ∈ [V h ] 2 the solution to (5).Then there exists h 0 > 0 such that for all h < h 0 with P e(h) 1 there holds By Cauchy-Schwarz combined with Lemma 6 and Young's inequality Applying the first equality of the consistency Lemma 3 we obtain From Lemma 6 and Young's inequality we thus have that for some ε 2 > 0, . Taking ε 2 < α/4 and combining the above bound with (31) we see that Since ε 1,2 are independent of h we can absorb them in the generic constant C and using the approximation inequality (21) together with Proposition 2, we conclude that where we used that P e(h) > 1.
4.4.Upstream estimates.In this case we consider β = (β 1 , 0) with β 1 < 0 and the data set ω = (0, x) × (y − , y + ) touching part of the outflow boundary ∂Ω + .We must choose the weight function differently and this time we take a negative ϕ given by ( 9) It seems that in this case we can not simultaneously get control of the L 2 -norm and the weighted H 1 -norm and we have to sacrifice the latter since it is not uniform in µ.We now take the weighted triple norm to be (32) ∂Ω − , and rederive the results obtained in Section 4.3, aiming for a local error estimate.Since ϕ ∈ (−1, 0), we will use that We start with an analogue of Lemma 5 by taking v h ϕ as a test function in the weak bilinear form a h and notice that since ϕ < 0 we now have that Arguing as previously in ( 24) but now for the outflow boundary, we obtain the bound and thus For the diffusion term we no longer have any positive contribution due to the change in sign of the weight function ϕ, since now Ω + (µ∇v h , v h ∇ϕ) Ω .We must therefore control this entirely using the stabilization.Integrating by parts and using the weighted trace inequality (19) To bound this by the triple norm we can simply use that |ϕ| < 1 and µ ≤ |β|h, giving that µ Hence we have that for some ε > 0, However, when P e(h)h > 1 one can obtain a better estimate due to µ ϕ .Summing these contributions we obtain the following result corresponding to Lemma 5. Lemma 7.There exists α > 0 such that for all v h ∈ V h we have when P e(h) > h −1 .Again, we can refine the control over the triple norm |||v h ||| ϕ by taking the projection π h (v h ϕ) ∈ V h as a test function and we obtain corresponding results.
Corollary 2. There exists α > 0 such that for all v h ∈ V h we have The argument in the proof of Corollary 1 remains valid with the remark that we now use the inequality |ϕ| < |ϕ| 1 2 .Lemma 8.For all v h ∈ V h there holds , when P e(h) > h −1 .Proof.We follow the proof of Lemma 6 and we focus on the bounds that are now different.As before, by the triangle inequality we have that up to a constant The first two terms can be bounded by C|||v h ||| ϕ as previously.For the third one, we can use the inverse inequality (12) and (18) to obtain Arguing as previously, we can bound the second line by C|||v h ||| ϕ using (33) instead of (24).We conclude the proof by recalling the estimate (30) for the jump term and the subsequent bounds.
We now prove the weighted error estimate in the upstream case showing that in the stability region ωβ we have quasi-optimal convergence for high Péclet numbers and a reduction of the convergence order by O(h 1 2 ) in an intermediate regime.
Theorem 2. Let u ∈ H 2 (Ω) be the solution of (1) and (u h , z h ) ∈ [V h ] 2 the solution to (5), then there holds Proof.We combine Lemma 7, Corollary 2 and Lemma 8 as in the proof of Theorem 1 and note that the argument holds verbatim when P e(h) > h −1 .Observe that when 1 P e(h) < h −1 we similarly obtain for some α > 0 and 0 < ε 1 < α/2, From Lemma 8 and Young's inequality we thus have that for some ε 2 > 0, . Taking ε 2 < α/4 and combining the above bound with (35) we see that ).Since ε 1,2 are independent of h we can absorb them the constant C and conclude the proof by using the approximation inequality (21) and Proposition 2 to obtain that The error bounds in this section contain a global Sobolev norm.This may be large in the presence of layers and it would be optimal to replace it with a local regularity measure.However, it is not clear how to reconcile this goal with the ill-posed character of the problem.Inserting everywhere the weight function as in the well-posed case [BGL09], would perturb the stability of the optimality system and, due to the lack of physical coercivity, the residual terms would then be hard to control.One could also consider a reduced transport equation (with zero diffusivity) and define the far field, where the solution is unknown, by a smooth extension, but is not obvious how this could be constructed in our context, and without requiring regularity for the right-hand side.Nonetheless, in numerical experiments we observe a good regularity behaviour, i.e. layers do not pollute the solution in the stability zone, as shown in Figures 15 and 16     We let Ω be the unit square and illustrate the performance of the numerical method (5) for different locations of the data domain ω and different regions of interest where we measure the approximation error.The computational domains are given in Figure 4 and the implementation is done using FEniCS [ABH + 15].In all the examples below we have used uniform triangulations with alternating left/right diagonals.In the definition of s Ω and s * we have taken the stabilization parameters γ = 10 −5 and γ * = 1, and ζ = 2 for s ω .The effect of different combinations of γ and γ * on the L 2 -errors is shown in Figure 5 and Figure 6 when data is given in a centered disk.Similar results are obtained when the data set is near the inflow/outflow boundary.Notice that our choice is empirically close to being optimal both locally and globally.We first show convergence plots both downstream and upstream from the data set when varying the diffusion coefficient µ and keeping the convection field β fixed.As Figure 10.Absolute L 2 -errors against mesh size h, computational domains in Figure 4a.Varying the diffusion coefficient µ for fixed β = (1, 0), exact solution u = 2 sin(5πx) sin(5πy).5.2.Interior data set.Next we consider the setting of the example discussed in the Introduction (Figure 2), where data is given in the centre of the domain.We give the convergence of the L 2 -errors in Figure 10 with the caveat that this location of the data set ω is not rigorously covered by the theoretical analysis of the previous sections.Nonetheless, the experiments are in agreement with the proven results.Notice that the L 2 -convergence is faster as µ decreases and for high Péclet numbers (above 10) one has optimal quadratic convergence both downstream and upstream, with the distinction that in the upstream case the convergence order is reduced in an intermediate regime, in agreement with the theoretical results in Section 4. Also, as expected from the error estimates proven in the first part [BNO20], when diffusion is moderately small one can see the transition towards the diffusion-dominated regime as the mesh gets refined -the convergence changes from almost quadratic to sublinear as the Péclet number decreases below 1. Figure 11 shows almost linear convergence in the H 1 -seminorm.We think this is observed due to the small contribution of the gradient term in the triple norm (23).We also remark almost no distinction between upstream and downstream for this approximation around the layer strongly deteriorates due to the crosswind position relative to the data sets.In this example the mesh is unstructured with 512 elements on a side and h ≈ 0.0025.