Analytic Reconstruction of a Two-Dimensional Velocity Field from an Observed Diffusive Scalar

TL;DR

This paper presents a method to uniquely reconstruct two-dimensional velocity fields from scalar measurements in incompressible flows, using high-resolution tracer data and PDE solving, with applications in oceanography.

Contribution

It introduces a novel approach to invert scalar fields for velocity reconstruction in 2D flows, extending to shallow-water equations and spatially varying diffusivity.

Findings

Unique velocity reconstruction possible with high-res scalar and transverse velocity measurements.

Method applicable to Navier-Stokes flows and ocean circulation models.

Potential for extending local velocity data to larger areas using scalar fields.

Abstract

Inverting an evolving diffusive scalar field to reconstruct the underlying velocity field is an underdetermined problem. Here we show, however, that for two-dimensional incompressible flows, this inverse problem can still be uniquely solved if high-resolution tracer measurements, as well as velocity measurements along a curve transverse to the instantaneous scalar contours, are available. Such measurements enable solving a system of partial differential equations for the velocity components by the method of characteristics. If the value of the scalar diffusivity is known, then knowledge of just one velocity component along a transverse initial curve is sufficient. These conclusions extend to the shallow-water equations and to flows with spatially dependent diffusivity. We illustrate our results on velocity reconstruction from tracer fields for planar Navier- Stokes flows and for a barotropic ocean circulation model. We also discuss the use of the proposed velocity reconstruction in oceanographic applications to extend localised velocity measurements to larger spatial domains with the help of remotely sensed scalar fields.

Equations

$$\frac{\partial c}{\partial t}+\mathbf{u}\cdot \nabla c=\mu {\nabla }^{2}c,$$

$$\frac{{\partial }^{2}c}{\partial t\partial x}+u\frac{{\partial }^{2}c}{\partial x^2}+\frac{\partial u}{\partial x}\frac{\partial c}{\partial x}+v\frac{{\partial }^{2}c}{\partial x\partial y}+\frac{\partial v}{\partial x}\frac{\partial c}{\partial y}=\mu \frac{\partial {\nabla }^{2}c}{\partial x},$$

$$\frac{{\partial }^{2}c}{\partial t\partial y}+u\frac{{\partial }^{2}c}{\partial x\partial y}+\frac{\partial u}{\partial y}\frac{\partial c}{\partial x}+v\frac{{\partial }^{2}c}{\partial y^2}+\frac{\partial v}{\partial y}\frac{\partial c}{\partial y}=\mu \frac{\partial {\nabla }^{2}c}{\partial y}.$$

$$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0$$

$$\frac{\partial v}{\partial x}\frac{\partial c}{\partial y}-\frac{\partial v}{\partial y}\frac{\partial c}{\partial x}=-u\frac{{\partial }^{2}c}{\partial x^2}-v\frac{{\partial }^{2}c}{\partial x\partial y}-\frac{{\partial }^{2}c}{\partial t\partial x}+\mu \frac{\partial {\nabla }^{2}c}{\partial x},$$

$$\frac{\partial u}{\partial x}\frac{\partial c}{\partial y}-\frac{\partial u}{\partial y}\frac{\partial c}{\partial x}=u\frac{{\partial }^{2}c}{\partial x\partial y}+v\frac{{\partial }^{2}c}{\partial y^2}+\frac{{\partial }^{2}c}{\partial t\partial y}-\mu \frac{\partial {\nabla }^{2}c}{\partial y}.$$

$$\frac{dx(s;t)}{ds}=\frac{\partial c(\mathbf{x}(s;t),t)}{\partial y},\frac{dy(s;t)}{ds}=-\frac{\partial c(\mathbf{x}(s;t),t)}{\partial x},$$

$$\frac{d\tilde{\mathbf{u}}(s;t)}{ds}=\mathbf{A}(s;t)\tilde{\mathbf{u}}(s;t)+\mathbf{b}(s;t),\tilde{\mathbf{u}}(s;t)\in {\mathbb{R}}^2,s\in \mathbb{R},$$

$$\tilde{\mathbf{u}}(s;t)=\mathbf{\Phi }(s;t)\tilde{\mathbf{u}}(0;t)+{\int }_0^s\mathbf{\Phi }(s-\sigma ;t)\mathbf{b}(\sigma ;t)d\sigma$$

$$\mathbf{A}(s;t)={\left[\begin{array}{cc}\frac{{\partial }^{2}c}{\partial x\partial y}& \frac{{\partial }^{2}c}{\partial y^2}\\ -\frac{{\partial }^{2}c}{\partial x^2}& -\frac{{\partial }^{2}c}{\partial x\partial y}\end{array}\right]}_{\mathbf{x}=\mathbf{x}(s;t)},\mathbf{b}(s;t)={\left[\begin{array}{c}\frac{{\partial }^{2}c}{\partial t\partial y}-\mu \frac{\partial {\nabla }^{2}c}{\partial y}\\ -\frac{{\partial }^{2}c}{\partial t\partial x}+\mu \frac{\partial {\nabla }^{2}c}{\partial x}\end{array}\right]}_{\mathbf{x}=\mathbf{x}(s;t)}.$$

$$\mu=\frac{1}{\nabla^{2}{c}}\Big{(}\frac{\partial c}{\partial t}+u\frac{\partial c}{\partial x}+v\frac{\partial c}{\partial y}\Big{)}.$$

$$\begin{array}{rl}\mathbf{A}(s;t)& ={\left[\begin{array}{cc}\frac{{\partial }^{2}c}{\partial x\partial y}-\frac{1}{{\nabla }^{2}c}\frac{\partial c}{\partial x}\frac{\partial {\nabla }^{2}c}{\partial y}& \frac{{\partial }^{2}c}{\partial y^2}-\frac{1}{{\nabla }^{2}c}\frac{\partial c}{\partial y}\frac{\partial {\nabla }^{2}c}{\partial y}\\ -\frac{{\partial }^{2}c}{\partial x^2}+\frac{1}{{\nabla }^{2}c}\frac{\partial c}{\partial x}\frac{\partial {\nabla }^{2}c}{\partial x}& -\frac{{\partial }^{2}c}{\partial x\partial y}+\frac{1}{{\nabla }^{2}c}\frac{\partial c}{\partial y}\frac{\partial {\nabla }^{2}c}{\partial x}\end{array}\right]}_{\mathbf{x}=\mathbf{x}(s;t)},\\ \mathbf{b}(s;t)& ={\left[\begin{array}{c}\frac{{\partial }^{2}c}{\partial t\partial y}-\frac{1}{{\nabla }^{2}c}\frac{\partial c}{\partial t}\frac{\partial {\nabla }^{2}c}{\partial y}\\ -\frac{{\partial }^{2}c}{\partial t\partial x}+\frac{1}{{\nabla }^{2}c}\frac{\partial c}{\partial t}\frac{\partial {\nabla }^{2}c}{\partial x}\end{array}\right]}_{\mathbf{x}=\mathbf{x}(s;t)}.\end{array}$$

$$\frac{\partial c}{\partial t}+\mathbf{u}\cdot \nabla c=\nabla (\mu \nabla c).$$

$$\mathbf{b}(s;t)={\left[\begin{array}{c}\frac{{\partial }^{2}c}{\partial t\partial y}-\frac{\partial \nabla (\mu \nabla c)}{\partial y}\\ -\frac{{\partial }^{2}c}{\partial t\partial x}+\frac{\partial \nabla (\mu \nabla c)}{\partial x}\end{array}\right]}_{\mathbf{x}=\mathbf{x}(s;t)}.$$

$$\mathbf{u}=\mathbf{u}(x,y,t).$$

$$c=c(x,y,t)$$

$$\partial (hu)/\partial x+\partial (hv)/\partial y=0.$$

$$h\frac{\partial c}{\partial t}+h\mathbf{u}\cdot \nabla c=h\nabla (\mu \nabla c).$$

$$d(h\tilde{\mathbf{u}})/ds=\mathbf{A}(h\tilde{\mathbf{u}})+\mathbf{b},$$

$$\mathbf{b}(s;t)={\left[\begin{array}{c}\frac{\partial }{\partial y}(h\frac{\partial c}{\partial y})-\frac{\partial }{\partial y}(h\nabla (\mu \nabla c))\\ -\frac{\partial }{\partial x}(h\frac{\partial c}{\partial y})+\frac{\partial }{\partial x}(h\nabla (\mu \nabla c))\end{array}\right]}_{\mathbf{x}=\mathbf{x}(s;t)}.$$

$$\psi(\mathbf{x})=\text{arctanh}\Big{[}\frac{\alpha\cos\sqrt{1+\alpha^{2}}x}{\sqrt{1+\alpha^{2}}\cosh{\alpha y}}\Big{]},$$

$$E_u=\frac{||u_{\text{reconst.}}-u_{\text{exact}}|{|}_2}{||u_{\text{exact}}|{|}_2},$$

$$\begin{array}{rl}u(x,y,0)=& \sin (x)\cos (y)+\sin (3x)\cos (3y),\\ v(x,y,0)=& -\cos (x)\sin (y)-\cos (3x)\sin (3y).\end{array}$$

$$c(\mathbf{x},t_0)=4\exp [-0.5((x-\pi {)}^2+(y-\pi {)}^2)],$$

$$N=A\tanh (\mathcal{N}(0,\sigma )),$$

$$Pe=\frac{UL}{\mu }.$$

Full text

1 Introduction and related work

The measurement of an evolving scalar tracer is sometimes considered easier than that of the underlying velocity (Schlichtholz, 1991), both in nature and in a laboratory setting. Determining the flow velocity more accurately from scalar field observations is, therefore, an exciting perspective in a plethora of applications including oceanography, meteorology, medicine and chemical engineering (Wunsch, 1987).

The transport equation for a passive scalar with instantaneous concentration $c(\mathbf{x},t)$ and constant diffusivity $\mu\geq 0$ under an incompressible, planar velocity field $\mathbf{u}(\mathbf{x},t)$ is given by

[TABLE]

with the spatial coordinates $\mathbf{x}=(x,y)$ varying over some domain $U\subset\mathbb{R}^{2}$ and $t$ denoting the time. Note that the velocity field $\mathbf{u}$ only influences the advection of $c$ via its inner product $\mathbf{u}\cdot\mathbf{\nabla}c$ with the instantaneous concentration gradient. This leads to the classic observation that the iso-scalar component of the velocity is inconsequential for the scalar field evolution (Dahm et al., 1992; Kelly, 1989). In other words, the iso-scalar velocity component falls in the null space of the velocity-reconstruction problem, necessitating further conditions for a unique solution to this inverse problem (Kelly, 1989).

Seminal work by Fiadeiro & Veronis (1984) and Wunsch (1985) considered velocity reconstruction from discretized observations of a steady scalar distribution in a two-dimensional, steady channel flow with uniform velocity. After spatial discretization, these authors obtained an under-determined algebraic system of equations for the values of the streamfunction at the grid points. The minimum-norm solution obtained to this problem by Fiadeiro & Veronis (1984) and Wunsch (1985) was found to be qualitatively inaccurate by Veronis (1986) even for the simple uniform flow considered. Providing additional boundary conditions (Fiadeiro & Veronis, 1984) or information on the flow anisotropy (Wunsch, 1985) leads to improved solutions. Alternatively, if the flow is dispensed with two scalars of different diffusivity, the system becomes over-determined because the same streamfunction must satisfy two linearly independent scalar equations. This reformulation, however, allows for a physically realistic solution that minimizes the least-squares of the solution error, as obtained by Fiadeiro & Veronis (1984). As expected, the reconstructed velocity departs from the true velocity significantly in regions where the two scalar fields are similarly distributed.

An alternative approach to inverting the scalar transport equation (1) has led to the emergence of Scalar Imaging Velocimetry (SIV) (see Wallace & Vukoslavčević 2010 for a review). When the diffusion of the scalar is slower than that of the carrier fluid (i.e., the ratio of fluid viscosity and scalar diffusivities is larger than one), the fluid velocity at nearby points varies slower than the scalar intensity (Dahm et al., 1992). As a consequence, the velocity component normal to the iso-scalar varies faster than the overall velocity. For an initial estimate, the overall velocity can therefore be assumed locally spatially constant and recovered from a projection onto neighboring normals to iso-scalar curves. As the calculation of iso-scalar-normal velocity requires division by $|\nabla c|$, the velocity reconstructed in this fashion is limited to regions where $|\nabla c|$ is well separated from zero (Dahm et al., 1992).

In an another approach relaxing the requirement of slow scalar diffusion, Su & Dahm (1996) formulated the velocity reconstruction as a minimization problem for an integral composed of the transport equation, the divergence and a regularizer term depending on the velocity gradient norm. Kelly (1989) also implemented a similar formulation for ocean flows in which the functional to be minimized involves the transport equation, the horizontal divergence, the relative vorticity and the kinetic energy. These minimization techniques can return an accurate (albeit never exact) velocity field if the weight factors for each component of the functional are tuned appropriately.

Here, we point out that for two-dimensional incompressible flows, one can use the method of characteristics to analytically reformulate the partial differential equation (PDE) (1) into a system of ordinary differential equations (ODEs) for the spatial evolution of the velocity components along an instantaneous level curve of $c(\mathbf{x},t)$. Integration of these ODEs enables the reconstruction of the velocity field over the domain of scalar observations. This procedure assumes that the velocity is available along a curve $\Gamma(t)$ transverse to the level curves of $c(\mathbf{x},t)$. If the value of the diffusivity $\mu$ is known, knowledge of one velocity component along $\Gamma(t)$ turns out to be sufficient. We give the derivation of these results in section 2 and test the approach on two numerical examples in section 3.

Our technique can be particularly advantageous in oceanographic applications, in which often only localised velocity field measurements are available. In such situations, independent observations of larger-scale physical, chemical or biological tracer fields offer an opportunity to extend the local velocity data to larger domains, as we illustrate in section 4. We also discuss oceanography-inspired extensions of our main results to diffusivities with spatial dependence and to shallow-water equations in section 2.

2 Formulation

Taking the spatial gradient of the scalar transport equation (1) leads to the equations

[TABLE]

There is no general recipe for solving such a coupled system of linear PDEs. Remarkably, however, the incompressibility condition

[TABLE]

enables us to rewrite this system as

[TABLE]

Both equations in system (4) share the same characteristics $\mathbf{x}(s;t)=(x(s;t),y(s;t))$, which satisfy the system of ODEs

[TABLE]

with the instantaneous time $t$ playing the role of a parameter. These characteristics are pointwise normal to the gradient $\nabla c$ and hence coincide with instantaneous level curves of $c(\mathbf{x},t)$ at any time $t$. Along the characteristics (5), the PDEs in (4) can be re-written as a linear, two-dimensional, non-autonomous system of ODEs for the velocity $\mathbf{\tilde{u}}(s;t)\mathrel{\mathop{:}}=\mathbf{u}(\mathbf{x}(s;t),t)$. This ODE is of the form

[TABLE]

with the form of the matrix $\mathbf{A}(s;t)\in\mathbb{R}^{2\times 2}$ and the vector $\mathbf{b}(s;t)\in\mathbb{R}^{2}$ depending upon wether the constant scalar diffusivity, $\mu$, is known or unknown, as we discuss below separately. In either case, once the normalized fundamental matrix solution $\boldsymbol{\Phi}(s;t)$ of the homogenous part $\frac{d\mathbf{\tilde{u}}}{ds}=\mathbf{A}(s;t)\mathbf{\tilde{u}}$, satisfying $\boldsymbol{\Phi}(0;t)=\mathbf{I}$, has been numerically determined, we can integrate (6) directly to obtain the velocity

[TABLE]

at any location $\mathbf{x}(s;t)$ along the level curve of $c(\mathbf{x},t)$ containing the point $\mathbf{x}(0;t)$.

2.1 Case of known scalar diffusivity

In the case where the value of $\mu$ is a priori known, equations (4) can be recast in the form (6) with the help of

[TABLE]

2.2 Case of unknown scalar diffusivity

If the value of $\mu$ is a priori unknown, we can still express $\mu$ from the original transport equation (1) as

[TABLE]

Substitution of this expression into (4) enables us to rewrite (4) in the form (6) with

[TABLE]

The formulae in equation (10), however, are only well-defined at points where $\nabla^{2}c\neq\mathbf{0}$.

2.3 Spatially-dependent diffusivity

Most non-eddy-resolving ocean circulation models are diffusion-based, i.e., assume that the small-scale, un- and under-resolved eddy motions lead to a diffusive tracer transfer—a process that can be fully characterized by just one parameter, the lateral eddy diffusivity $\mu$. A number of diffusivity estimates are available based on data from surface drifters and tracers (Okubo, 1971; Zhurbas & Oh, 2004; LaCasce, 2008; Rypina et al., 2012; LaCasce & Bower, 2000; Lumpkin et al., 2002; McClean et al., 2002) satellite-observed velocity fields (Marshall et al., 2006; Abernathey & Marshall, 2013; Klocker & Abernathey, 2014; Rypina et al., 2012) and Argo float observations (Cole et al., 2015). All these observations suggest that the ocean eddy diffusivity is strongly spatially varying with up to two orders of magnitude difference between regions with quiescent vs. energetic oceanic eddy regimes. Eddy-resolving simulations also demonstrate highly nonuniform spatial distribution of the eddy-induced transport (Gille & Davis, 1999; Nakamura & Chao, 2000; Roberts & Marshall, 2000). Thus, any method that involves diffusivity assumption in a global settings, including the velocity reconstruction method described here, would need to take into account spatial variation in diffusivity.

In this case, the appropriate advection-diffusion equation takes the form

[TABLE]

Taking the gradient of this equation and following the same steps as in Section 2, one obtains that the velocity $\tilde{\mathbf{u}}$ satisfies equation (6) along the iso-curves of tracer concentration $c$, with $\mathbf{A}$ from eq. (8) and

[TABLE]

2.4 Shallow-water approximation

Due to negligible vertical velocities, eq. (3) is a reasonable approximation for oceanic surface flows on spatial scales larger than the mesoscale (Pedlosky, 2013; Salmon, 1998; Bennett, 2005). When vertical velocities become non-negligible, one may still find another variable other than the velocity that satisfies eq. (3) and re-derive the inversion for that variable. An example is the inviscid shallow-water model with a rigid lid (Pedlosky, 2013; Batchelor, 1967), as we shall discuss next. This shallow-water approximation is often employed in physical oceanography to describe ocean flows in shallow tidally-driven bays, such as, for example, the Katama Bay located on the island of Martha’s Vineyard, MA (Slivinski et al., 2017; Orescanin et al., 2016; Luettich Jr & Westerink, 1991).

Consider a shallow sheet of inviscid unstratified fluid overlaying spatially-dependent bathymetry. Fluid parcels are allowed to have three non-zero components of velocity, $\{u,v,w\}\neq 0$, but if the horizontal velocity $\mathbf{u}$ is initially independent of $z$ then it remains so for all times, i.e., we may write

[TABLE]

If the initial tracer distribution is independent of $z$, i.e., $c_{0}=c_{0}(x,y)$, then we will also have

[TABLE]

for all times. We then invoke a rigid-lid approximation, assuming that the free surface elevation $\eta$—the deviation of the free surface from its mean value—is much smaller than the water depth $h(x,y)$, i.e., $\eta\ll h\ll L$, where $L$ is the characteristic horizontal length scale of motion. Then the water depth is time-independent and can be considered known because the bathymetry of a shallow bay can be mapped out and the sea surface height is assumed to be constant due to the rigid lid approximation.

In this case, mass conservation requires the transport velocity, $h\mathbf{u}$, to be divergence free:

[TABLE]

To proceed, we multiply the advection- diffusion equation (11) by $h$ to obtain

[TABLE]

We take the gradient of this equation, then follow the steps of section 2 to obtain an analogue of eq. (6) for the transport velocity, given by

[TABLE]

with $\mathbf{A}$ from eq. (8) and with

[TABLE]

As before, equation (17) should be solved along the iso-contours of $c$, which are the characteristics of the PDE (16).

2.5 Velocity reconstruction

For an experimentally or numerically obtained scalar field, $c(\mathbf{x},t)$, the characteristic $\mathbf{x}(s;t)$ emanating from any point $\mathbf{x}_{0}(t)\mathrel{\mathop{:}}=\mathbf{x}(0;t)$ can be computed numerically from equation (5). To recover the velocity $\mathbf{\tilde{u}}(s;t)$ at a point $\mathbf{x}(s;t)$ of such a characteristic via (7), we need to know the velocity $\mathbf{\tilde{u}}(0;t)$ at the initial location $\mathbf{x}_{0}(t)$. The PDE (4) is well-posed only if such a set of initial velocity vectors is available along a non-characteristic curve $\Gamma(t)$, i.e., along a curve segment that is transverse to the characteristics at each of its points. If, however, the diffusivity $\mu$ is a priori available, then only one component of $\mathbf{\tilde{u}}(0;t)=\mathbf{{u}}(\mathbf{x}_{0},t)$ needs to be known along the non-characteristic curve $\Gamma(t)$, because the other component can be determined from (1). For example, if only the $u(\mathbf{x},t)$ velocity is known at a point of $\mathbf{x}\in\Gamma(t)$, then $v(\mathbf{x},t)$ can be uniquely determined from (1) as long as the contour of $c(\mathbf{x},t)$ is not locally parallel to the $x$ direction, i.e., as long as ${\partial c}(\mathbf{x},t)/{\partial y}\neq 0$.

To reconstruct the velocity field $\mathbf{u}(\mathbf{x},t)$ over the largest possible domain, the initial velocities must therefore be known along a curve $\Gamma(t)$ transverse to the level curves of $c(\mathbf{x},t)$ over the largest possible domain. Most of the characteristics (level curves) in our examples will be closed, but this is not a requirement for our method to be applicable. Indeed, to obtain $\mathbf{u}(\mathbf{x},t)$ along a non-closed level set of $c(\mathbf{x},t)$, one simply has to integrate equation (5) both in the forward and the backward $s$-direction until the domain boundary is reached.

The formulae (8)-(18) contain the first-order temporal derivative of the observed scalar $c(\mathbf{x},t)$, as well as up to second and up to third-order spatial derivatives, respectively. This presents challenges for velocity reconstruction over domains where only scarce information is available on $c(\mathbf{x},t)$. We envision, however, applications of this velocity reconstruction procedure in situations where high-resolution imaging of $c(\mathbf{x},t)$ is available, enabling an accurate computation of the necessary derivatives.

In geophysical flows, the ratio of advective transport rate to the diffusive rate (i.e. the dimensionless Péclet number) tends to be large. For such large Péclet numbers, in the formulae relevant for velocity reconstruction under a known diffusivity (eqs. (8), (10), (12) and (18)), errors and uncertainties arising from computing the necessary third derivatives are suppressed by the small $\mu$ values multiplying these terms. We illustrate this effect in section 4.2 with respect to substantial under-sampling and moderate uncertainties.

3 Numerical examples

We now carry out our proposed velocity reconstruction scheme in two examples. In both, we generate $c(\mathbf{x},t)$ numerically by evolving an initial scalar distribution $c_{0}(\mathbf{x})\mathrel{\mathop{:}}=c(\mathbf{x},t_{0})$ under the PDE (1) with a known velocity field $\mathbf{u}(\mathbf{x},t)$, then reconstruct $\mathbf{u}(\mathbf{x},t)$ solely from observations of $c(\mathbf{x},t)$. We solve for $c(\mathbf{x},t)$ by discretizing (1) using the pseudospectral method. The necessary spatial derivatives of $c(\mathbf{x},t)$ in the reconstruction formulae are evaluated in the frequency domain using MATLAB’s FFT algorithm; the time integration is performed in real space using MATLAB’s variable-order adaptive time step solver, ODE45. The standard 2/3 dealiasing is applied where the energy in the last one third of the wavenumbers is artificially removed. In both examples, the scalar diffusivity is taken to be $\mu=0.01$ and a 512 $\times$ 512 grid is used. We perform the integration of the characteristic ODEs (equations (5) and (6)) using MATLAB’s ODE113 to achieve high accuracy.

3.1 Example 1: Steady, analytic, inviscid velocity field

We consider the Mallier–Maslowe vortices (Mallier & Maslowe, 1993; Gurarie & Chow, 2004) which are exact solutions of the 2D Euler equations. Modeling a periodic row of vortices with alternating rotation directions, the corresponding velocity field derives from the streamfunction

[TABLE]

which we plot for reference in the left panel of figure 1 for $\alpha=1$. \floatsetup[figure]style=plain,subcapbesideposition=top

Over the square domain $\mathbf{x}\in[0,2\sqrt{2}\pi]\times[-\sqrt{2}\pi,\sqrt{2}\pi]$, we generate the diffusive scalar distribution $c(\mathbf{x},t)$ from the initial concentration field $c(\mathbf{x},0)=4\exp[-0.2((x-\sqrt{2}\pi)^{2}+y^{2})]$ by solving the transport equation (1). No advective transport is possible among the cells defined by the streamfunction in the left panel of figure 1, but diffusion enables the transport of the scalar $c$ across cell boundaries. We will reconstruct the underlying velocity field generated by (19) at the time $t=4.1$. The initial and the evolved scalar fields are shown in the middle and right panels of figure 1. The non-characteristic initial curve $\Gamma(t)$ for the characteristics is formed by a set of 100 points on the solid red curve in the right panel of figure 1, along which we consider the velocities known.

Next, we reconstruct the velocity field using both the known- and the unknown-diffusivity formulation. While integrating equation (1), we record velocities along each characteristic at 3000 locations, irrespective of their lengths. We define the relative normed error over the computational domain as the ratio of the $L_{2}$ norm of the difference between the reconstructed and the original velocity field to the $L_{2}$ norm of the original velocity field defined in equation (19). For the ${u}$ velocity component, this error norm, $E_{u}$, can be written as

[TABLE]

where $u_{\text{reconst.}}$ and $u_{\text{exact}}$ are the reconstructed and exact $u$-velocities, respectively. The error, $E_{v}$, for the $v$ velocity component is defined similarly. We have found $E_{u}$ and $E_{v}$ for the known-diffusion case (not shown) to be $1.6\times 10^{-3}$ and $1.2\times 10^{-3}$, respectively, and for the unknown-diffusion case to be $1.7\times 10^{-2}$ and $1.1\times 10^{-2}$, respectively. For reference, we also show exact-interpolated velocities in the right subplots in figure 2 and 2, i.e., velocities evaluated exactly along the 3000 locations along the characteristics, then interpolated onto the regular grid used for visualization.

Figure 2 shows that throughout the domain spanned by the characteristics emanating from the inital-condition curve $\Gamma(t)$ (solid red line in the right panel of figure 1), there is generally close agreement between the reconstructed and the exact velocity profiles. The domain of the velocity reconstruction is limited by the weak mixing of the scalar field $c(\mathbf{x},t)$ in this steady flow. The initial condition curve in the right panel of figure 1 only intersects the scalar contours forming the three middle cells in the left panel of figure 1. Hence, only these cells appear in the reconstructed velocities on the left in figures 2 and 2.

3.2 Example 2: Unsteady Navier–Stokes velocity field

We now apply our velocity reconstruction procedure to an unsteady flow obtained from the direct numerical solution of the unforced, two-dimensional Navier–Stokes equations. We solve the vorticity-streamfunction formulation of the equations, as described, e.g., in Melander et al. (1987), over the doubly periodic square domain $\left(x,y\right)\in[0,2\pi]\times[0,2\pi]$ with the non-dimensional kinematic viscosity $\nu=10^{-5}$. Hence, the scalar-diffusion characterized by $\mu=0.01$ is faster than the vorticity-diffusion, in contrast to the requirements of the SIV method of Dahm et al. (1992). We initialize the numerical simulation with the velocity field

[TABLE]

Eventually, two distinct regions of opposite vorticity with different strengths emerge, subsequently evolving periodically within the domain, as shown in figure 3. For longer simulation times ($t\approx 500$ here), these final structures exhibit slow decay due to dissipation and weak non-linearity (Matthaeus et al., 1991).

Using the velocity field from $t_{0}=490$ onwards, we advect a scalar field $c(\mathbf{x},t)$ with initial distribution

[TABLE]

and with diffusivity $\mu=0.01$ (cf. figure 4). The unsteady velocity field visibly leads to intense mixing compared to Example 1.

We have reconstructed the velocity field for various times between $t_{0}=490$ and $t=499.5$ with similar outcomes in all cases. Here we only show the result for $t=499.5$ in figures 5 and 6 for the solid red and dashed black lines of figure 4 as initial curves, respectively.

3.2.1 Solid red initial curve

If the initial curve $\Gamma(t)$ for the characteristic ODEs (5)-(6) is the solid red curve in figure 4, then the reconstructed velocity field captures most features of the true velocity field. For the known diffusivity formulation, the reconstruction errors $E_{u}$ and $E_{v}$ are $1.9\times 10^{-3}$ and $2.5\times 10^{-3}$, respectively. For the unknown diffusivity formulation, we have obtained $E_{u}=5.2\times 10^{-2}$ and $E_{v}=5.9\times 10^{-2}$. For the more erroneous latter case, the overall qualitative accuracy is still compelling, as we show in figure 5. Most of the local inaccuracies arise from interpolation errors, due to a lack of penetration of the external scalar level curves into the cores of high-vorticity regions (cf. the lower-right region of figure 4). This is confirmed by a comparison of the exact-interpolated and the exact plots in figures 5 and 5. While the spatial domain of the reconstructed velocity is limited by the presence of scalar transport barriers, the numerical calculation is overall accurate (cf. the reconstructed and exact-interpolated plots of figure 5).

3.2.2 Dashed–black curve

For this case (with the scalar diffusivity, $\mu$, assumed unknown), the reconstructed velocity field recovers the positive high-vorticity region (bottom right of figure 3), fully capturing even the finer features (cf. the reconstructed vs. the exact velocity field in figure 6). However, the negative vorticity region in the top left of figure 3 is not fully captured in the reconstructed velocity plots of figure 6. The error from this calculation is nevertheless still low: $E_{u}=6.92\times 10^{-2}$ and $E_{v}=6.23\times 10^{-2}$.

4 Oceanographic applications

Oceanic velocities are key variables in understanding a number of environmental processes that range from climate change and global warming to transport of pollutants and various bio/geo/chemical tracers. Obtaining reliable estimates for ocean currents with high accuracy and resolution, however, remains a challenging task.

Oceanic velocities can be measured directly by current meters or Acoustic Doppler Current Profiler (ADCP) instruments attached to ships, stationary moorings, or towed or autonomous sampling platforms. Such estimates provide vertical profiles or sections of velocities at select locations or along select ship tracks but only at limited times and locations. Ship-based ADCP sections, however, can provide the required initial conditions for our proposed velocity reconstruction method along all tracer contours intersecting the ship track.

Surface drifters also provide estimates of near-surface ocean currents. These semi-Lagrangian instruments consist of a surface-drifting buoy equipped with a GPS tracker and a subsurface drogue, or sail, underneath. Finite differencing consecutive GPS fixes with respect to the time interval between neighboring GPS transmissions provides estimates for the underlying ocean currents at the depths of the drogue, typically 20m. The Global Surface Drifting Buoy Array maintained by NOAA is part of their Global Ocean Observing System (http://www.aoml.noaa.gov/phod/gdp/objectives.php). Currently, this drifter array has 1445 operational drifters that are distributed throughout the World Oceans with at least 1 drifter per each 5 x 5 deg bin at any given time. Although the 5 degree resolution of the resulting velocity field obtained from these drifters via finite differencing appears coarse, it would yield a much higher-resolution global estimate of the near-surface oceanic velocities when used as a boundary condition for the velocity inversion proposed here along all observed tracer contours passing through drifters positions.

A third alternative for collecting ocean current data is indirect remote sensing from either land-based or satellite-based radars. For example, satellite-based estimates of sea surface height provide nearly global geostrophic velocities at the ocean surface (cf. https://www.aviso.altimetry.fr). However, these velocity estimates only provide the geostrophic component of the total velocity field and only on spatial scales larger than meso-scale. Land-based high-frequency radars, on the other hand, can measure full surface velocity at a higher spatio-temporal resolution, but only over limited spatial domains.

Satellite-based measurements of diffusive tracers–such as sea surface temperature, sea surface salinity, or ocean color–are also available. Applying the velocity reconstruction technique described in this paper to these tracer distributions, along with the available information about the localized velocity values from current meters, ADCPs, drifters, moorings or HF radars, would potentially allow for the estimation of the oceanic velocities in a much larger (possibly global) domain.

Clearly, the opportunities we have described above for ocean current reconstruction come with their own challenges, including non-diffusive decay of tracers, sporadic availability of the tracers in space and time, presence of noise in data, lack of exact two-dimensionality and incompressibility, and spatial dependence of the diffusivity. Some of these challenges we have already addressed in sections 2.3 and 2.4, sensitivity to noise and undersampling will be explored in section 4.2, and the rest we leave to future work. As a proof of concept, we consider below a numerical example in which we apply our proposed velocity inversion to a conservative tracer in a barotropic ocean circulation model.

4.1 Application to an ocean circulation model

We seek to reconstruct velocities from a tracer distribution generated numerically by the ocean circulation model MITgcm (Marshall et al., 1997). We keep the complexity of the model to the minimum and use the simplest form of the advection-diffusion equation with a constant diffusivity coefficient. In its linear limit, the model is configured to reproduce the solution described by the Munk model of the wind-driven ocean circulation (Munk, 1950), though the numerical model allows for nonlinear features to develop, such as eddies and jets.

The simulation has been configured on a beta-plane with a rigid lid, 3.5 km spatial resolution, a temporal step of 200 seconds and a momentum diffusion coefficient of 40 m2s*-1*. The fluid is unstratified and initially is at rest, but is subject to a temporally constant stress at the surface boundary that drives the gyre flows that develop. Momentum input by the wind at the surface is removed by boundary friction at the lateral boundaries.

Consistent with the Munk model, we mimic the time-mean-observed winds with an idealized structure that is zonally uniform and exerts stress only in the zonal direction. The winds have a cosine structure in the meridional direction, with westerly winds at mid latitudes and easterly winds to the north and south. The net vorticity input into the domain is zero, with the negative vorticity input into the northern half of the domain balanced by the positive vorticity in the south of the domain. A 2-gyre ocean circulation develops, with intense western boundary currents mimicking the Gulf Stream and Labrador currents of the Atlantic Ocean, and the Kuroshio and Oyashio of the Pacific. The two gyres are separated by a meandering jet –an analogue of the Gulf Stream or Kuroshio Extension currents.

After an initial spin up of 300 days, a tracer with a constant meridional gradient was released into the model, and its evolution over the subsequent 180 days was simulated numerically using an advection-diffusion scheme with constant diffusivity coefficient $\mu=100$ m2s*-1*. Three consecutive tracer snapshots 200 s apart were saved on day 180 and were used to numerically estimate matrices $\mathbf{A}$ and $\mathbf{b}$ in equation (8), for use with the velocity inversion equation (6). All tracer derivatives were computed using finite differences, in a manner similar to what is feasible in ocean applications involving remotely sensed velocity and tracer data.

We have applied the velocity inversion method with known-scalar diffusivity (eqs. (6)-(8)) to reconstruct velocities along several tracer iso-contours covering different parts of the model domain. We show a comparison between the reconstructed velocity and the ground truth (actual model velocity on day 480) in figure 7 for three representative tracer contours. In this calculation, we assume that the true velocity—our initial condition for integrating eq. (6)—is known at the northernmost point of the contour, then perform the integration going southward along the entire length of the contour. The reconstructed velocities provide useful information for oceanographic studies in roughly the right 2/3 of the model domain, far enough from the western boundary where the eddy activity is most vigorous and the tracer contours are most filamented. Close to the western boundary, in the left 1/3 of the domain, the reconstructed velocity diverges from the true model velocity due to numerical errors arising in the computation of $\mathbf{A}$ and $\mathbf{b}$ and in the integration of the matrix equation (6) along the tracer contours. The reconstructed velocity is still useful along the initial segment of the iso-contour before the error accumulation causes significant mismatch. The reconstruction performs better if we integrate from the equator going poleward to the north and south (thus decreasing the length of the tracer contour two times), or if the true velocity is known at more than one point along the tracer contour. In that case, we can apply the velocity inversion to the segments of the contour instead of the entire contour.

In the case when the spatio-temporal resolution of the scalar field is significantly higher than the dominant length and/or time-scales of the flow, the reconstruction can be improved by smoothing the tracer field in space and/or time to suppress noise arising during the calculation of tracer derivatives. In our example, applying a spatial smoothing, specifically, the two-dimensional convolution filter ‘conv2’ in MATLAB with a constant-coefficient 6x6 square kernel (equivalent to a two dimensional running average over $\sim$21 km2), allows for reasonably accurate velocity reconstruction further into the left 1/3 of the domain (figure 8).

4.2 The effect of resolution and uncertainty in the data

While our analytic velocity reconstruction is, in principle, exact under conditions described in section 2, its practical implementations will have to address sparsity and uncertainty in the available observational data. A general analysis of these issues with real ocean data is beyond the scope of this paper, but we provide an initial illustration of the sensitivities to noise and resolution for the oceanographic model we have studied.

We illustrate the outcome of velocity reconstruction from under-resolved tracer fields in figures 9a and 10, where we sub-sampled the full model tracer field over an increasingly sparse grid. The sub-sampled field was then interpolated to its original grid using MATLAB’s build-in cubic interpolator, and the velocity inversion was applied along the same level contour of the tracer, as in the middle panel of figure 7. An example of the velocity reconstruction using an under-resolved tracer field is shown in figure 9a.

Our analysis suggest that the inversion is reasonably insensitive to under-sampling as the contour-averaged reconstructed velocity error (expressed as a percentage of the range of respective velocity along the iso-contour) in both zonal and meridional direction (figure 10) is less than $5\%$ for grids with spacing up to about 70 km which is 20 times the original grid spacing. The critical grid spacing is likely dependent on the kinetic energy spectrum and the dominant length scales of the flow.

Since real data is noisy, we have also investigated the applicability of the velocity inversion to model fields with increasing noise levels. We add noise pointwise,

[TABLE]

where $A$ is the noise amplitude and $\mathcal{N}(0,\sigma)$ is a zero mean Gaussian random variable with standard deviation $\sigma$ (generated using MATLAB’s ‘randn’ function), to the model tracer field. Unlike a pure Gaussian signal, (23) allows for defining a maximal noise amplitude normally guaranteed for a remote-sensing instrument. As a first step in our reconstruction algorithm, we use a linear regression model using the 37 time slices (about 2 hours) centered around the current time to temporally filter the input scalar-field data. We then smooth the noisy tracer concentration fields in space using MATLAB’s ‘conv2’ filter with constant coefficient square kernels, both before and after calculating a scalar derivative required in the reconstruction (cf. equation (8)). For purely spatial derivatives of any order, a 17$\times$17 kernel is used in all the examples presented below, which correspond to about 58 km. The size of spatial filter for terms involving temporal derivatives is optimised for each case.

Figure 9b shows the velocity reconstruction along the tracer contour in the middle panel of 7 but using the noisy scalar data with noise amplitude $A$ equal to $10\%$ of the local temporal mean concentration across the 37 time slices (equivalent to a 2-hour period, which is effectively the observation time for this example) and setting a standard deviation $\sigma$ equal to $0.01$ in the definition (23). The shape of the iso-contour remains almost the same as that without the noise. The kernel size for spatial filtering of terms involving the temporal derivatives is 30$\times$30 or, equivalently, 103 km. For these parameters, we have found the contour-averaged errors relative to the model velocities in the zonal and meridional direction to be 5.6% and 7.8%.

In the remote sensing of scalars, such as sea-surface temperature through satellite microwave radiometers, noise levels up to 18% may be present (Gentemann et al., 2010). Aerial measurements of dye from multispectral and hyperspectral cameras have similar noise levels (about 20%) for Rhodamine WT concentrations of about 20ppb (Hally-Rosendahl et al., 2015; Clark et al., 2014). Dye releases of this magnitude on spatio-temporal scales of several km and several hours have been successfully carried out both in the surf-zone (Hally-Rosendahl et al., 2015; Clark et al., 2014) and in the coastal ocean (Rypina et al., 2016). Figure 9c shows the reconstructed velocity with 20% noise amplitude and the rest of the parameters being the same as figure 9b. The kernel size for spatial filtering of terms involving the temporal derivatives is 45$\times$45 (or 155 km). The contour-averaged errors in zonal and meridional velocity for this case are 9% and 10.8%. The joint effect of 20% amplitude noise and under-sampling (20 times less resolution) is shown in figure 9d. The contour-averaged errors in the zonal and meridional velocity are 10% and 10.5% for this case.

The presence of second and third order derivatives in equation (8) would usually make our reconstruction method susceptible to noise in the data. In the example presented above, however, we find that $\mathcal{O}(\frac{\partial^{2}c}{\partial t\partial x},\frac{\partial^{2}c}{\partial t\partial y})\gg\mathcal{O}(\mu\frac{\partial{\nabla^{2}c}}{\partial x},\mu\frac{\partial{\nabla^{2}c}}{\partial y})$ holds throughout the domain, so the right-hand side of equation (8) is dominated by the second—rather than the third—derivative terms. Recall that for a characteristic velocity, $U$, and characteristic length scale, $L$, the Péclet number for the scalar transport equation (1) is defined as

[TABLE]

Using the empirical relation between length scale $L$ [cm] and diffusivity $\mu$ [cm2/s] obtained through synthesis from various dye release experiments in the ocean (Okubo, 1971), we can let $\mu=0.0103L^{1.15}$. Therefore, for typical velocity values, $U=10$–100 cm/s, the Péclet number estimated for length scales, $L=0.5$–104 km ranges from 61 to 2706. Evolution of tracers in flows with large Péclet numbers is advection-dominated, rendering the terms involving the third derivatives in (8) small. Thus, in our velocity reconstruction procedure, the sensitivity to noise would enter largely through the second derivatives in eqs. (8) or (12). Therefore, the most effort for noise suppression in the purely spatial derivatives via smoothing or filtering should be directed towards these terms.

In the noisy scalar examples presented in figures 9b to 9d, the major source of error is the first-order temporal scalar gradient, $\partial c/\partial t$. Using noisy scalar data to evaluate purely spatial derivative terms of all orders in equation (8) and using model scalar data for the temporal derivatives, we have found that an accurate velocity reconstruction (where contour average velocity errors remains less than about 10%) can be obtained for much higher noise levels (both $A$ and $\sigma$ in equation (23)). Hence, the higher-order derivatives in the formulation of our velocity reconstruction method do not necessarily pose a challenge for accurate velocity reconstruction from noisy scalar data.

In summary, our numerical results indicate that the velocity reconstruction proposed here is of value to physical oceanography when applied to tracer fields with spatial resolution of several km, temporal resolution of minutes to hours, and noise levels of up to 20% of the local temporal mean and a standard deviation of 0.01 for the noise model represented by equation (23). Such resolution can be achieved for airplane-mounted radars but is currently beyond the capabilities of the existing satellites that have resolution from tens to hundreds of km spatially and from days to weeks temporally.

Even though our idealized model was able to reproduce several of the major qualitative features of the real ocean circulation in the Atlantic and Pacific ocean, a more realistic baroclinic free-surface model configuration is required in order to better simulate real oceanic conditions. We also note that because sea surface temperature and salinity are influenced by vertical mixing, air-sea interaction, and precipitation/evaporation/river run off, these two tracer fields are only purely diffusive on short temporal scales. Additional challenges will arise in real ocean applications that are associated with small but non-zero horizontal divergence of surface velocity, which was exactly zero in our numerical model due to the rigid lid approximation, the space-dependent diffusivity, and the presence of higher-level non-Gaussian noise in data.

5 Conclusion

It is well established that unique velocity reconstruction from an observed scalar field is an underdetermined problem (Dahm et al., 1992; Kelly, 1989). We have pointed out here, however, that if the velocity at a point is known, then the velocity field all along the scalar level curve emanating from that point can be uniquely expressed using sufficiently well-resolved observations of the scalar. In a well-mixing flow, therefore, velocities can in principle be reconstructed on global domains from highly localized velocity measurements combined with global imaging of a scalar field.

While our procedure also works with unknown scalar diffusivities, we have found in our examples the reconstruction errors $E_{u}$ and $E_{v}$ to be about an order of magnitude smaller when the diffusivity is a priori known. This is due to the numerical errors in the evaluation of higher scalar derivatives involved in the unknown diffusivity formulation (cf. equation (10) vs. (8)). To further improve accuracy for the latter formulation, once the velocity field is recovered, a pointwise estimate for the scalar diffusivity could be obtained at each point in the domain spanned by the scalar level curves emanating from an initial curve. A subsequent statistical analysis could then establish an acceptable estimate for $\mu$ from these pointwise values, which would in turn enable use of the known-diffusivity formulation of our approach. Depending on accuracy requirements, this modification could even be turned into an iterative procedure, involving a gradual, simultaneous refinement of both the velocity field and the diffusivity estimate.

Extensions of our results to three-dimensional flows may be possible but remain mathematically challenging because of the lack of a general recipe for solving multi-dimensional systems of linear partial differential equations. An extension to general compressible flows appears beyond reach because incompressibility or specific form of shallow water equations secures the common set of characteristic curves critical to our analytic solution strategy. A case with temporally constant but spatially varying density of the fluid in 2D will have a similar formulation as the shallow water case.

The approach proposed here has the potential to extend spatially localized ocean velocity measurements along contours of satellite-observed tracer fields, such as temperature, salinity and ocean color (phytoplankton) measurements. In Section 4, we have provided a first proof of concept by reconstructing larger-scale velocities generated by a global circulation model, with the tracer-field output of the same model serving as observational input to our procedure. We have also demonstrated the robustness of this reconstruction under sparsification of the observed scalar field and the addition of noise that models observational inaccuracies. More work is required to investigate applicability of this method to real ocean data. Further work can build on the extension of our theory we have given in sections 2.3-2.4 for spatially dependent diffusivities and 2D-divergent but mass-conserving shallow water velocities.

We acknowledge useful conversations with Markus Holzner and Larry Pratt, as well as partial funding from the Turbulent Superstructures priority program of the German National Science Foundation (DFG), and from the NASA grant #NNX14AH29G to IR.