Stability estimate for scalar image velocimetry

In this paper we analyse the stability of the system of partial differential equations modelling scalar image velocimetry. We first revisit a successful numerical technique to reconstruct velocity vectors ${u}$ from images of a passive scalar field $\psi$ by minimising a cost functional, that penalises the difference between the reconstructed scalar field $\phi$ and the measured scalar field $\psi$, under the constraint that $\phi$ is advected by the reconstructed velocity field ${u}$, which again is governed by the Navier-Stokes equations. We investigate the stability of the reconstruction by applying this method to synthetic scalar fields in two-dimensional turbulence, that are generated by numerical simulation. Then we present a mathematical analysis of the nonlinear coupled problem and prove that, in the two dimensional case, smooth solutions of the Navier-Stokes equations are uniquely determined by the measured scalar field. We also prove a conditional stability estimate showing that the map from the measured scalar field $\psi$ to the reconstructed velocity field $u$, on any interior subset, is H\"older continuous.


Introduction
Scalar image velocimetry (SIV) is the reconstruction of the fluid velocity field u from measurements of a scalar field ψ, that is advected by u.This idea, which dates back to the works [10,23,24], is applied in weather forecasting models using e.g.satellite images of clouds or ocean temperature.For a background of the technique see for instance [14] and references therein.
SIV also finds applications in medical flow imaging and in experimental fluid mechanics.For instance in the laser induced fluorescence (LIF) technique, a fluid is seeded with fluorescent molecules and laser light is focused into a thin sheet, where it is absorbed by the fluorescent molecules followed by spontaneous emission of light which is then recorded by a camera (see e.g.[21] and references therein).Assuming that the recorded light intensity is proportional to the fluorescence concentration ψ, the velocity field u := u 1 e x + u 2 e y , with e x , e y the Cartesian unit vectors, can be reconstructed by invoking the scalar transport equation, where λ is the scalar diffusivity.A direct inversion of the scalar transport equation only provides the component u ⊥ of u that is normal to ψ isolines, Finding all components of u requires additional constraints, for example, conservation of hydrodynamic variables.
Inspired by recent developments on computational methods for SIV [13], in this work we study mathematically the stability of the map from ψ to u.We show, in particular, that under suitable a priori assumptions, the velocity field is uniquely determined by the measured scalar field, thereby giving a partial theoretical justification to the computational approaches.The analysis uses conditional stability estimates as a workhorse for understanding of inverse problems in a spirit similar to that of the seminal work by Bukhgeȋm and Klibanov [4].
To fix a configuration for the stability analysis we consider a soap film experiment from [13].In that work, a soap film with a thickness of ∼ 10µ m is formed in the 10 cm gap between two vertical parallel nylon wires.Turbulence in the film is generated by piercing the film with an array of cylindrical obstacles, resulting in chaotically interacting wake vortices.These perturbations are accompanied by film thickness fluctuations which behave similar as those of a passive scalar.The soap film is illuminated and light reflections are recorded at high speed.The recorded interference pattern depends on the film thickness fluctuations and therefore behaves similar as a passive scalar.
Motivated by this experiment, we consider a particular Cartesian geometry, Ω = (−a, a) × (−b, b), with a, b > 0, and write Q = Ω × (0, T ) for T > 0. Our proofs generalize to many other two dimensional settings but the above choice allows for convenient notations.We assume that the flow satisfies the slip or no-slip conditions on the vertical boundaries, x = ±a, and that the scalar field satisfies the no flux condition there, but the boundary conditions for both the velocity u and the scalar field ψ on the top and bottom boundaries y = ±b are unknown.
Determination of u given ψ is possible only if ψ satisfies some non-degeneracy condition.Indeed, a constant ψ satisfies (1) but gives no information on u.To state our main result in a simplified form (see Theorem 2 below for the precise formulation), we assume that the spatial derivative of ψ does not vanish identically on the left boundary x = −a at any time.
The simplified formulation is as follows.Let u and ũ be smooth velocity fields satisfying the Navier-Stokes equations, together with, say, the no-slip boundary condition on x = ±a, and let ψ and ψ be smooth scalar fields satisfying (1) with u and ũ, respectively, together with the no flux boundary condition on x = ±a.Suppose, moreover, that for all t ∈ (0, T ) there is y ∈ (−b, b) such that ∂ y ψ(−a, y, t) = 0. Then for every space time domain ω, such that ω ⊂ Q (here and below X denotes the closure of the set X), there exists α ∈ (0, 1) and C > 0 such that the following stability estimate holds whenever u and ũ satisfy an a priori bound.Observe that this estimate in particular shows that the scalar field ψ uniquely determines the velocity field u.It also gives an upper bound similar to those applied in [2,7] to obtain sharp error estimates for finite element data assimilation methods.Although the estimate (2) can not be directly applied in that context, it is a first step towards an analysis of the error propagation in computational velocity reconstruction based on SIV.
Our approach combines a stream function formulation for the two dimensional Navier-Stokes' equations, with the classical pressure velocity formulation.First, using a stream function Θ to represent the velocity field u, the convection-diffusion equation for the scalar field φ defines a transport equation for Θ.Provided the scalar field satisfies the non-degeneracy condition ∂ y ψ = 0 on the left boundary, this transport equation can be solved in a neighbourhood of the left boundary.This way we show that u can be reconstructed locally given the scalar field ψ.Then to extend this local reconstruction to an arbitrary subdomain in the interior of the space time domain, we apply a unique continuation result for the transient linearized Navier-Stokes' equations [16].Observe that the local reconstruction step is possible whenever the velocity u is known along a curve segment transverse to the level curves of ψ.In our configuration this holds on the lateral boundaries thanks to the chosen boundary conditions.A variant of this local reconstruction approach, using classical techniques, has recently been applied to SIV in [20].
The outline of the paper is as follows.First in section 2 we recall a recent technique for computational velocity reconstruction, and study its stability computationally in Section 3. Then in Section 4 we prove the conditional stability estimate.The paper ends with concluding remarks in Section 5.

An example of a computational method for SIV
In this section we describe a variational method for the reconstruction of the space and time dependent velocity field u, pressure field p and scalar field φ from a measured space and time dependent scalar field ψ, also referred to as the true scalar field.We assume perfect knowledge of ψ over a space time volume of Ω × (0, T ), where Ω and (0, T ) represent the corresponding spatial and temporal dimensions.The reconstruction method divides (0, T ) into segments.The the i-th segment starts and ends at t 0 = (i − 1)τ and t 1 = t 0 + τ , respectively, where τ is referred to as the segment time.The reconstruction scheme solves a sequence of optimisation problems for the unknown state variable w = (u, p, φ) at the start of each segment, i.e. at t = t 0 .We use a subscript on a field variable to indicate a time instance, e.g.w 0 = w(t 0 ).Finding w 0 in each segment involves an iterative scheme, and the initial guess for the iteration is taken from the reconstructed field w 1 at t = t 1 obtained in the preceding segment.
2.1.The cost function and its minimisation.Defining the state variable as: w = (u, p, φ), the velocity reconstruction method finds the initial conditions for the state variable w 0 in each time segment, by minimising the deviation between the reconstructed scalar field φ and the measured scalar field ψ integrated over Ω × (0, T ).The corresponding cost functional for the method reads: where the norm • is based on the L 2 (Ω) inner product, denoted by •, • .Equation ( 3) is minimised under the constraint that w 1 is related to w 0 via the conservation equations of fluid momentum, fluid mass and scalar field: Here ν is the fluid kinematic viscosity and λ is scalar diffusivity.Adding constraint (4) to equation (3) using the Lagrange multiplier ŵ = (û, p, φ) results in the following constrained cost functional, which is referred to as the Lagrangian L : Minimising L of (5) w.r.t.w 0 involves computing the gradient of L w.r.r w 0 , i.e. δL /δw 0 , following [13] we obtain: This expression for δL /δw 0 contains the Lagrange multiplier ŵ = (û, p, φ), whose evolution equation and initial and final conditions read: (7a) Equation (7a) governs the evolution of the Lagrange multiplier ŵ = (û, p, φ), showing that û is incompressible and is forced along gradients of φ, and that both û and φ are advected by u and are subjected to diffusion.The diffusion coefficients −ν and −λ of these transport equations are negative, and therefore these equations are integrated backward in time from t = t 1 to t = t 0 .The 'starting' conditions ŵ1 at t = t 1 are given by equation (7b), while the 'final' conditions ŵ0 at t = t 0 define the optimisation update direction of w 0 via (6).This direction approaches zero, when w 0 reaches an extremum of L , which corresponds to the condition of (7c).
To find w 0 we use the Polak-Rebiere variant of the conjugate gradient method [18], which updates w 0 along a search direction h related to δL /δw 0 .The initial guess for w 0 is w 1 from the previous time segment, and the step length along h is varied using Brent's line minimisation algorithm [3], until the minimum of the functional J , from (3), in this direction is found.The conjugate gradient algorithm is continued until the relative change in J between two consecutive iterations drops below 0.01.A reconstruction typically require ∼ 10 2 conjugate gradient steps and ∼ 10 Brent minimisation steps per conjugate gradient step.Therefore the computational effort of both methods is equivalent to that of ∼ 10 3 computational fluid dynamics simulations.
It was shown in [13] that the reconstruction method method produces unstable results when the segment time τ is larger than the flow correlation time T , which for the system described above corresponds to τ 0.1.The instability is related to the ill posed-ness of the initial value problem for chaotic systems, which corresponds to the cost functional developing multiple minima, when τ exceeds T , see for instance [17].In order to stabilise the method for these cases regularisation terms may be added to the functional of (3).In this work we restrict ourselves to τ 0.1, which does not require the use of these regularisation terms.

Numerical investigation of the stability of SIV
3.1.Setup.We apply the computational velocity reconstruction to a two dimensional, incompressible, decaying turbulent flow field in a space time domain Ω × (0, T ).The spatial domain Ω = (−π, π) 2 is a bi-periodic square with side length L = 2π and the temporal domain (0, T ) has a size of T = 8, which is referred to as the reconstruction time.The objective is to reconstruct the velocity field u from the measured space and time dependent scalar field ψ.To distinguish between the reconstructed velocity field and the reference velocity field, giving the advection of ψ, we denote the former by u and the latter by v.It is re-emphasised that we have access to perfect information of ψ on Ω × (0, T ).Both ψ and v start from random initial conditions and the initial conditions are normalised, such that U = v = 1 and ψ = 1 at t = 0.The diffusivity is λ = 2 × 10 −3 and the viscosity is ν = 1 × 10 −3 , which corresponds to a Reynolds number of Re = U L/ν = 6.3 × 10 3 based on the initial velocity scale U and a Schmidt number of Sc = ν/λ = 1/2.We furthermore use a segment time of τ = 8 × 10 −2 .
The true velocity field v and the true (measured) scalar field ψ are generated on a 256 2 grid by numerically integrating (4).Subsequently the scalar measurement is interpolated to a 128 2 grid and the reconstruction fields u and φ are obtained on this 128 2 grid by iteratively integrating equations (4,6,7).On both the 256 2 and 128 2 grids, equations (4, 7a) are advanced in time using a computational time step of ∆t = 10 −3 .Spatial derivatives in these equations are computed using the Fourier basis functions.Time integration is performed using the second-order explicit Adams-Bashforth scheme for the advection terms and the second-order implicit Crank-Nicolson scheme for the diffusion terms.For further details on the numerical methods see [11,12,13].
3.2.Results. Figure 1(left) shows the reconstruction error ǫ = u − v|| 2 as a function of time, where it is recalled, that u is the reconstructed velocity and v is the reference velocity field.It is seen, that with time, ǫ first decreases and after t > 4, the error levels off.
Recalling the time segmentation described in the first paragraph if Section 2, the time dependent behaviour in Figure 1(left) indicates that the reconstruction depends on the quality of the initial guess for the initial condition at the start of each segment.In the first segment the initial guess is zero, while in each consecutive segment the initial guess becomes closer to the reference solution, explaining the observed decrease in ǫ with time.Also recall that the stability of the unique continuation problem for parabolic equations degenerates towards t = 0, so a similar effect is expected also for globally coupled solutions.
We now study the convergence of the reconstructed u as a function of the true (measured) ψ.To this end we generate multiple simulations of true scalar fields ψ that are advected by different 'true' velocity fields v.For these simulations the conditions at t = 0 for ψ are identical but the conditions for v differ.Differences in a quantity q between two simulations are denoted by q − q.The relative difference in the initial conditions is varied between 10 −5 and 10 −1 .To study the convergence of u as a function of ψ, we plot in Figure 1(right) u − ũ 2 as a function of ψ − ψ 2 .It is noted that we average this norm over a time interval of 6 < t < 8.During this time interval the computational method has converged roughly to a steady state, as show in Figure 1(left).
Figure 1(right) shows that there are two regimes in the dependence of the variation of the reconstructed velocity u − ũ 2 on the variation of the true scalar These regimes are indicated by the dashed lines in Figure 1(right) which have a slope of zero and one, respectively.In these regimes, the variation of the true scalar ψ − ψ 2 (and of the true velocity v − ṽ 2 ) is smaller and larger, respectively, than the error of the reconstructed velocity u − v 2 .Consequently, the variation of the reconstructed velocity u − ũ 2 in these two regimes is set by u − v 2 (constant) and by ψ − ψ 2 , respectively.The level of stagnation depends on the time and space resolution of the Navier-Stokes' solver, which limits the accuracy of the reconstructed velocity.The second regime of Eq. ( 8) corresponds to Lipschitz stability of the velocities.In view of the classical stability result [8] for data assimilation problems subject to the heat equation, Lipschitz stability is feasible thanks to the fact that only the initial data of u are unknown in this computational example.We also show that the local reconstruction step is Lipschitz stable, see Remark 1 below.We leave a precise numerical analysis of the computational error as a topic of a future work, see [7] for such an analysis of data assimilation problems subject to the heat equation.
For illustration purposes, we show in Figure 2(top) two true scalar fields, ψ and ψ, at t = 6, that are advected by two 'true' velocity fields v and ṽ.At t = 0 it was assumed that ψ = ψ and that ṽ = v + ∆v where ∆v / v = 10 −1 .Figure 2(bottom) shows the curl of the corresponding reconstructed velocity fields, u and ũ, at t = 6.

Stability analysis
In this section we will derive a stability estimate that gives a mathematical interpretation of the stability illustrated by the right plot of Figure 1.Compared to the numerical example we consider a more challenging setting where not only the initial data but also the boundary data on parts of the boundary are unknown.
To fix the ideas we consider a simplified geometric configuration similar to that of a two-dimensional soap film experiment, described above.Let u be a solution to the Navier-Stokes equations in Q: On Σ we assume that u satisfies either homogeneous Dirichlet conditions, u = 0, or slip conditions, i.e. u • n = 0, where n is the outward pointing normal on Σ.
On the top and bottom boundaries the boundary data are unknown.The viscosity coeffcient ν > 0 is known.We suppose that a passive scalar ψ satisfies in Q, together with the no flux boundary condition Observe that since u 1 = 0 on Σ this reduces to a Neumann condition in practice.The diffusivity coefficient λ > 0 is a known constant.
Recall that the inverse problem to find u given ψ is clearly unsolvable if ψ is a constant function, so we make the standing assumption that this is not the case.More precisely we will assume that ψ is non-constant in space on Σ.
It is well-known that the Navier-Stokes equations admit smooth solutions in the two dimensional case, see e.g.[22,Remark 3.7,p. 303].For simplicity, we assume that both u and ψ are smooth in Ω × (0, T ).
Note that smoothness of ψ in the interior of Q actually follows from the interior Schauder estimates for (10) and smoothness of u, see e.g.[15,Theorem 8.12.1,p. 131].However, as we do not assume any boundary conditions on the horizontal boundaries y = ±b, we do not enter into discussion of smoothness properties of ψ.

Stream function.
As Ω is a two-dimensional simply connected domain, the vanishing of divergence ∇ • u = 0 implies that there is Θ such that Observe that (10) implies on Ω, we can write equivalently viewing the vector field X as the differential operator defined by (13).
As ∂ y Θ = u 1 = 0 on Σ, due to the boundary condition on u, and as the stream functions Θ and Θ + c, with c ∈ R a constant, give the same u, we may assume that We can view ( 14)-( 15) as a transport equation for Θ.Observe that the vector field X and the right-hand side ζ are known as ψ is known.

4.2.
On stability of linear transport equations.In this section we study linear transport equations in an abstract setting and the meaning of the variables here are different from those above.
Consider the equation with the initial condition u| t=0 = 0. We are interested in the continuity properties of the map (f, β) → u.However, let us first recall the standard continuity result for the map f → u with β fixed.The exposition below was inspired by [9,Theorem 7,p. 131] where more complicated Hamilton-Jacobi equations are considered.Contrary to this reference, we need to keep track of the dependence on β of the constants in the estimates. Let Furthermore, define r(t) = R(T − t), let B(r) ⊂ R n be the open ball of radius r, with a fixed centre and outward pointing normal n B , and consider the energy Then for a solution u of ( 17), where Indeed, using we get where we used |n B • β| ≤ R = −r ′ .Now (18) follows from the Cauchy-Schwarz inequality.We write Grönwall's inequality implies the following lemma.
for all u ∈ C 1 (C) satisfying (17) and u| t=0 = 0.The constant C depends only on T and the norm of β.
We need also higher order energy estimates.For notational simplicity, let us now consider only the case n Applying (18) to v gives where the constant depends only on the L ∞ (0, T ; W k,∞ (R)) norm of β.The sum of the estimates (20) for k = 0, . . ., K gives Grönwall's inequality gives the following lemma.
for all u ∈ C k+1 (C) satisfying ( 17) and u| t=0 = 0.The constant C depends only on T and the norm of b.
We can apply a similar argument also to the time derivative.We will need only the first derivative in time v = ∂ t u.Then Writing Ẽ = B(r) v 2 dx and applying (18) to v gives where the constant depends only on W 1,∞ ((0, T ) × R) norm of β.The sum of this and ( 21) with K = 1, together with Grönwall's inequality, gives the following lemma.
Lemma 3. Let β ∈ W 1,∞ ((0, T ) × R).Then there is C > 0 such that for all u ∈ C 2 (C) satisfying (17) and u| t=0 = 0.The constant C depends only on T and the norm of β. Then where the constant depends only on T and R 1 .
Proof.The function w = u 1 − u 2 satisfies and from Lemma 3, Using f 2 H 1 (C) ≤ R 1 , the claim follows from Lemma 3.

4.4.
Global recovery.We recall the three cylinders inequality for the linearized Navier-Stokes equation from [16], see also [1] for an earlier related result.Let u ∈ H 1 (0, T ; H 2 loc (Ω)) be a nontrivial solution of (9) with associated pressure p ∈ L 2 (0, T ; H 1 (Ω)).Let B(x, R) ⊂ R 2 denote the ball of radius R with the centre x.

Concluding remarks
We have shown that for the SIV problem, the velocity field u is uniquely determined by the measured scalar field ψ, and u depends continuously on ψ.The stability is of Hölder type for the interior estimates considered here.Due to the nonlinearity of the map ψ → u, the scalar field in the right hand side of the stability estimate of Theorem 2 is measured in the H 4 -norm.This is much stronger than the L 2 -norm of the velocities in the left hand side, but it seems unlikely that it can be improved by much in the framework exposed here.The consequence of this lack of balance in the estimate is that computationally we must expect the error in the velocity to be larger than that in the scalar field, even if α ≈ 1.
An outstanding challenge is to further develop the analysis so that it allows for error estimates for a computational method.Several building blocks for such a development are available, for convection-diffusion equations and transport in [6,5], for parabolic problems in [7] and for the stationary linearized Navier-Stokes' equation in [2].In those references finite element methods are considered, but the arguments can be reinterpreted in the context of spectral or Fourier methods that we considered here.

Figure 1 .
Figure 1.(Left) Reconstruction error versus time.(Right) Variation in reconstructed velocity versus variation in measured scalar.
. Remark 1.Using notation from the above proposition, it is immediate from equation (12) that