2.1 Case of known scalar diffusivity

In the case where the value of $\unicode[STIX]{x1D707}$ is a priori known, equations (2.3) can be recast in the form (2.5) with the help of

(2.7a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D63C}(s;t)=\left[\begin{array}{@{}cc@{}}{\displaystyle \frac{\unicode[STIX]{x2202}^{2}c}{\unicode[STIX]{x2202}x\unicode[STIX]{x2202}y}} & {\displaystyle \frac{\unicode[STIX]{x2202}^{2}c}{\unicode[STIX]{x2202}y^{2}}}\\ -{\displaystyle \frac{\unicode[STIX]{x2202}^{2}c}{\unicode[STIX]{x2202}x^{2}}} & -{\displaystyle \frac{\unicode[STIX]{x2202}^{2}c}{\unicode[STIX]{x2202}x\unicode[STIX]{x2202}y}}\end{array}\right]_{\boldsymbol{x}=\boldsymbol{x}(s;t)},\quad \boldsymbol{b}(s;t)=\left[\begin{array}{@{}c@{}}{\displaystyle \frac{\unicode[STIX]{x2202}^{2}c}{\unicode[STIX]{x2202}t\unicode[STIX]{x2202}y}}-\unicode[STIX]{x1D707}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FB}^{2}c}{\unicode[STIX]{x2202}y}}\\ -{\displaystyle \frac{\unicode[STIX]{x2202}^{2}c}{\unicode[STIX]{x2202}t\unicode[STIX]{x2202}x}}+\unicode[STIX]{x1D707}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FB}^{2}c}{\unicode[STIX]{x2202}x}}\end{array}\right]_{\boldsymbol{x}=\boldsymbol{x}(s;t)}.\end{eqnarray}$$

2.2 Case of unknown scalar diffusivity

If the value of $\unicode[STIX]{x1D707}$ is a priori unknown, we can still express $\unicode[STIX]{x1D707}$ from the original transport equation (1.1) as

(2.8) $$\begin{eqnarray}\unicode[STIX]{x1D707}={\displaystyle \frac{1}{\unicode[STIX]{x1D6FB}^{2}c}}\left({\displaystyle \frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}t}}+u{\displaystyle \frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}x}}+v{\displaystyle \frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}y}}\right).\end{eqnarray}$$

Substitution of this expression into (2.3) enables us to rewrite (2.3) in the form (2.5) with

(2.9) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}rcl@{}}\unicode[STIX]{x1D63C}(s;t)\,\, & =\,\, & \left[\begin{array}{@{}cc@{}}{\displaystyle \frac{\unicode[STIX]{x2202}^{2}c}{\unicode[STIX]{x2202}x\unicode[STIX]{x2202}y}}-{\displaystyle \frac{1}{\unicode[STIX]{x1D6FB}^{2}c}}{\displaystyle \frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}x}}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FB}^{2}c}{\unicode[STIX]{x2202}y}} & {\displaystyle \frac{\unicode[STIX]{x2202}^{2}c}{\unicode[STIX]{x2202}y^{2}}}-{\displaystyle \frac{1}{\unicode[STIX]{x1D6FB}^{2}c}}{\displaystyle \frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}y}}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FB}^{2}c}{\unicode[STIX]{x2202}y}}\\ -{\displaystyle \frac{\unicode[STIX]{x2202}^{2}c}{\unicode[STIX]{x2202}x^{2}}}+{\displaystyle \frac{1}{\unicode[STIX]{x1D6FB}^{2}c}}{\displaystyle \frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}x}}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FB}^{2}c}{\unicode[STIX]{x2202}x}} & -{\displaystyle \frac{\unicode[STIX]{x2202}^{2}c}{\unicode[STIX]{x2202}x\unicode[STIX]{x2202}y}}+{\displaystyle \frac{1}{\unicode[STIX]{x1D6FB}^{2}c}}{\displaystyle \frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}y}}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FB}^{2}c}{\unicode[STIX]{x2202}x}}\end{array}\right]_{\boldsymbol{x}=\boldsymbol{x}(s;t)},\\ \boldsymbol{b}(s;t)\,\, & =\,\, & \left[\begin{array}{@{}c@{}}{\displaystyle \frac{\unicode[STIX]{x2202}^{2}c}{\unicode[STIX]{x2202}t\unicode[STIX]{x2202}y}}-{\displaystyle \frac{1}{\unicode[STIX]{x1D6FB}^{2}c}}{\displaystyle \frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}t}}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FB}^{2}c}{\unicode[STIX]{x2202}y}}\\ -{\displaystyle \frac{\unicode[STIX]{x2202}^{2}c}{\unicode[STIX]{x2202}t\unicode[STIX]{x2202}x}}+{\displaystyle \frac{1}{\unicode[STIX]{x1D6FB}^{2}c}}{\displaystyle \frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}t}}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FB}^{2}c}{\unicode[STIX]{x2202}x}}\end{array}\right]_{\boldsymbol{x}=\boldsymbol{x}(s;t)}.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

The formulae in (2.9), however, are only well-defined at points where $\unicode[STIX]{x1D6FB}^{2}c\neq 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 $\unicode[STIX]{x1D707}$ . A number of diffusivity estimates are available based on data from surface drifters and tracers (Okubo Reference Okubo1971; LaCasce & Bower Reference LaCasce and Bower2000; Lumpkin, Treguier & Speer Reference Lumpkin, Treguier and Speer2002; McClean et al. Reference McClean, Poulain, Pelton and Maltrud2002; Zhurbas & Oh Reference Zhurbas and Oh2004; LaCasce Reference LaCasce2008; Rypina et al. Reference Rypina, Kamenkovich, Berloff and Pratt2012), satellite-observed velocity fields (Marshall et al. Reference Marshall, Shuckburgh, Jones and Hill2006; Rypina et al. Reference Rypina, Kamenkovich, Berloff and Pratt2012; Abernathey & Marshall Reference Abernathey and Marshall2013; Klocker & Abernathey Reference Klocker and Abernathey2014) and Argo float observations (Cole et al. Reference Cole, Wortham, Kunze and Owens2015). 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 versus energetic oceanic eddy regimes. Eddy-resolving simulations also demonstrate highly non-uniform spatial distribution of the eddy-induced transport (Gille & Davis Reference Gille and Davis1999; Nakamura & Chao Reference Nakamura and Chao2000; Roberts & Marshall Reference Roberts and Marshall2000). Thus, any method that involves a diffusivity assumption in a global setting, 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

(2.10) $$\begin{eqnarray}{\displaystyle \frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}t}}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}c=\unicode[STIX]{x1D735}(\unicode[STIX]{x1D707}\unicode[STIX]{x1D735}c).\end{eqnarray}$$

Taking the gradient of this equation and following the same steps as in § 2, one obtains that the velocity $\tilde{\boldsymbol{u}}$ satisfies equation (2.5) along the iso-curves of tracer concentration $c$ , with $\unicode[STIX]{x1D63C}$ from (2.7) and

(2.11) $$\begin{eqnarray}\boldsymbol{b}(s;t)=\left[\begin{array}{@{}c@{}}{\displaystyle \frac{\unicode[STIX]{x2202}^{2}c}{\unicode[STIX]{x2202}t\unicode[STIX]{x2202}y}}-{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D735}(\unicode[STIX]{x1D707}\unicode[STIX]{x1D735}c)}{\unicode[STIX]{x2202}y}}\\ -{\displaystyle \frac{\unicode[STIX]{x2202}^{2}c}{\unicode[STIX]{x2202}t\unicode[STIX]{x2202}x}}+{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D735}(\unicode[STIX]{x1D707}\unicode[STIX]{x1D735}c)}{\unicode[STIX]{x2202}x}}\end{array}\right]_{\boldsymbol{x}=\boldsymbol{x}(s;t)}.\end{eqnarray}$$

2.4 Shallow-water approximation

Due to negligible vertical velocities, equation (2.2) is a reasonable approximation for oceanic surface flows on spatial scales larger than the mesoscale (Salmon Reference Salmon1998; Bennett Reference Bennett2005; Pedlosky Reference Pedlosky2013). When vertical velocities become non-negligible, one may still find another variable other than the velocity that satisfies (2.2) and re-derive the inversion for that variable. An example is the inviscid shallow-water model with a rigid lid (Batchelor Reference Batchelor1967; Pedlosky Reference Pedlosky2013), 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 (Luettich Jr & Westerink Reference Luettich and Westerink1991; Orescanin, Elgar & Raubenheimer Reference Orescanin, Elgar and Raubenheimer2016; Slivinski et al. Reference Slivinski, Pratt, Rypina, Orescanin, Raubenheimer, MacMahan and Elgar2017).

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 $\boldsymbol{u}$ is initially independent of $z$ then it remains so for all times, i.e. we may write

(2.12) $$\begin{eqnarray}\boldsymbol{u}=\boldsymbol{u}(x,y,t).\end{eqnarray}$$

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

(2.13) $$\begin{eqnarray}c=c(x,y,t)\end{eqnarray}$$

for all times. We then invoke a rigid-lid approximation, assuming that the free-surface elevation $\unicode[STIX]{x1D702}$ – the deviation of the free surface from its mean value – is much smaller than the water depth $h(x,y)$ , i.e. $\unicode[STIX]{x1D702}\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\boldsymbol{u}$ , to be divergence free:

(2.14) $$\begin{eqnarray}\unicode[STIX]{x2202}(hu)/\unicode[STIX]{x2202}x+\unicode[STIX]{x2202}(hv)/\unicode[STIX]{x2202}y=0.\end{eqnarray}$$

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

(2.15) $$\begin{eqnarray}h{\displaystyle \frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}t}}+h\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}c=h\unicode[STIX]{x1D735}(\unicode[STIX]{x1D707}\unicode[STIX]{x1D735}c).\end{eqnarray}$$

We take the gradient of this equation, then follow the steps of § 2 to obtain an analogue of (2.5) for the transport velocity, given by

(2.16) $$\begin{eqnarray}\text{d}(h\tilde{\boldsymbol{u}})/\text{d}s=\unicode[STIX]{x1D63C}(h\tilde{\boldsymbol{u}})+\boldsymbol{b},\end{eqnarray}$$

with $\unicode[STIX]{x1D63C}$ from (2.7) and with

(2.17) $$\begin{eqnarray}\boldsymbol{b}(s;t)=\left[\begin{array}{@{}c@{}}{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}y}}\left(h{\displaystyle \frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}y}}\right)-{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}y}}(h\unicode[STIX]{x1D735}(\unicode[STIX]{x1D707}\unicode[STIX]{x1D735}c))\\ -{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}}\left(h{\displaystyle \frac{\unicode[STIX]{x2202}c}{\unicode[STIX]{x2202}y}}\right)+{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}}(h\unicode[STIX]{x1D735}(\unicode[STIX]{x1D707}\unicode[STIX]{x1D735}c))\end{array}\right]_{\boldsymbol{x}=\boldsymbol{x}(s;t)}.\end{eqnarray}$$

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

2.5 Velocity reconstruction

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

To reconstruct the velocity field $\boldsymbol{u}(\boldsymbol{x},t)$ over the largest possible domain, the initial velocities must therefore be known along a curve $\unicode[STIX]{x1D6E4}(t)$ transverse to the level curves of $c(\boldsymbol{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 $\boldsymbol{u}(\boldsymbol{x},t)$ along a non-closed level set of $c(\boldsymbol{x},t)$ , one simply has to integrate equation (2.4) both in the forward and the backward $s$ -direction until the domain boundary is reached.

The formulae (2.7)–(2.17) contain the first-order temporal derivative of the observed scalar $c(\boldsymbol{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 on $c(\boldsymbol{x},t)$ is available. We envision, however, applications of this velocity reconstruction procedure in situations where high-resolution imaging of $c(\boldsymbol{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 ((2.7), (2.9), (2.11) and (2.17)), errors and uncertainties arising from computing the necessary third derivatives are suppressed by the small $\unicode[STIX]{x1D707}$ values multiplying these terms. We illustrate this effect in § 4.2 with respect to substantial under-sampling and moderate uncertainties.

3.1 Example 1: Steady, analytic, inviscid velocity field

We consider the Mallier–Maslowe vortices (Mallier & Maslowe Reference Mallier and Maslowe1993; Gurarie & Chow Reference Gurarie and Chow2004) which are exact solutions of the two-dimensional Euler equations. Modelling a periodic row of vortices with alternating rotation directions, the corresponding velocity field derives from the streamfunction

(3.1) $$\begin{eqnarray}\unicode[STIX]{x1D713}(\boldsymbol{x})=\text{arctanh}\left[{\displaystyle \frac{\unicode[STIX]{x1D6FC}\cos \sqrt{1+\unicode[STIX]{x1D6FC}^{2}}x}{\sqrt{1+\unicode[STIX]{x1D6FC}^{2}}\cosh \unicode[STIX]{x1D6FC}y}}\right],\end{eqnarray}$$

which we plot for reference in figure 1(a) for $\unicode[STIX]{x1D6FC}=1$ .

Figure 1.

Contour plot of the steady-state streamfunction for Mallier–Maslowe vortices with $\unicode[STIX]{x1D6FC}=1$ in (3.1) (a), and the evolution of the initial concentration field (b), along these vortices into the one at $t=4.1$ (c), for the scalar diffusivity value $\unicode[STIX]{x1D707}=0.01$ .

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

Figure 2.

Reconstructed, exact and exact-interpolated velocities for Example 1 at $t=4.1$ , with the diffusivity assumed unknown: (a) $u(\boldsymbol{x},4.1)$ , (b) $v(\boldsymbol{x},4.1)$ . The reconstructed velocities have been interpolated from the characteristics (level sets in figure 1 c) to a uniform grid. The exact-interpolated velocities are obtained by evaluating equation (3.1) over the characteristics at $t=4.1$ , then interpolating to the same uniform grid used for the reconstructed velocities. This procedure is designed to factor out interpolation errors that would be present in a direct comparison with the exact velocity field evaluated directly over the numerical grid.

Next, we reconstruct the velocity field using both the known- and the unknown-diffusivity formulation. While integrating equation (1.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 (3.1). For the $u$ velocity component, this error norm, $E_{u}$ , can be written as

(3.2) $$\begin{eqnarray}E_{u}={\displaystyle \frac{||u_{reconst.}-u_{exact}||_{2}}{||u_{exact}||_{2}}},\end{eqnarray}$$

where $u_{reconst.}$ and $u_{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 figure 2(a), 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 $\unicode[STIX]{x1D6E4}(t)$ (solid red line in figure 1 c), 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(\boldsymbol{x},t)$ in this steady flow. The initial condition curve in figure 1(c) only intersects the scalar contours forming the three middle cells in figure 1(a). Hence, only these cells appear in the reconstructed velocities on the left in figure 2(a,b).

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, McWilliams & Zabusky (Reference Melander, McWilliams and Zabusky1987), over the doubly periodic square domain $(x,y)\in [0,2\unicode[STIX]{x03C0}]\times [0,2\unicode[STIX]{x03C0}]$ with the non-dimensional kinematic viscosity $\unicode[STIX]{x1D708}=10^{-5}$ . Hence, the scalar diffusion characterized by $\unicode[STIX]{x1D707}=0.01$ is faster than the vorticity diffusion, in contrast to the requirements of the SIV method of Dahm et al. (Reference Dahm, Su and Southerland1992). We initialize the numerical simulation with the velocity field

(3.3) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle u(x,y,0)=\sin (x)\cos (y)+\sin (3x)\cos (3y),\\ \displaystyle v(x,y,0)=-\cos (x)\sin (y)-\cos (3x)\sin (3y).\end{array}\right\}\end{eqnarray}$$

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 nonlinearity (Matthaeus et al. Reference Matthaeus, Stribling, Martinez, Oughton and Montgomery1991).

Figure 3.

Evolution of vorticity field for Example 2: (a) $\unicode[STIX]{x1D714}(\boldsymbol{x},490)$ and (b) $\unicode[STIX]{x1D714}(\boldsymbol{x},499.5)$ .

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

(3.4) $$\begin{eqnarray}c(\boldsymbol{x},t_{0})=4\exp [-0.5((x-\unicode[STIX]{x03C0})^{2}+(y-\unicode[STIX]{x03C0})^{2})],\end{eqnarray}$$

and with diffusivity $\unicode[STIX]{x1D707}=0.01$ (cf. figure 4). The unsteady velocity field visibly leads to intense mixing compared to Example 1.

Figure 4.

Evolution of the diffusive scalar $c(\boldsymbol{x},t)$ under the unforced, two-dimensional Navier–Stokes flow initialized by the velocity field (3.3) of Example 2: (a) $c(\boldsymbol{x},490)$ , (b) $c(\boldsymbol{x},499.5)$ . The solid red curve and the dashed black curve shown over the diffused concentration field are two different initial (non-characteristic) curves for the velocity reconstruction.

Figure 5.

Reconstructed, exact and exact-reconstructed velocities for Example 2 (with the diffusivity, $\unicode[STIX]{x1D707}$ , treated as unknown) at $t=499.5$ , using the red solid line from figure 4(b) as the initial curve: (a) $u(\boldsymbol{x},499.5)$ , (b) $v(\boldsymbol{x},499.5)$ .

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.2 Dashed – black curve

For this case (with the scalar diffusivity, $\unicode[STIX]{x1D707}$ , assumed unknown), the reconstructed velocity field recovers the positive high-vorticity region (bottom right of figure 3 b), fully capturing even the finer features (cf. the reconstructed versus the exact velocity field in figure 6). However, the negative vorticity region in the top left of figure 3(b) 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}$ .

Figure 6.

Velocities for Example 2 (diffusivity, $\unicode[STIX]{x1D707}$ , assumed unknown) at $t=499.5$ using the black dashed line from figure 4(b) as the initial curve: (a) $u(\boldsymbol{x},499.5)$ , (b)  $v(\boldsymbol{x},499.5)$ . The plotting is done identically to figure 5.

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. Reference Marshall, Adcroft, Hill, Perelman and Heisey1997). 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 Reference Munk1950), 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 s and a momentum diffusion coefficient of $40~\text{m}^{2}~\text{s}^{-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 $\unicode[STIX]{x1D707}=100~\text{m}^{2}~\text{s}^{-1}$ . Three consecutive tracer snapshots 200 s apart were saved on day 180 and were used to numerically estimate matrix $\unicode[STIX]{x1D63C}$ and vector $\boldsymbol{b}$ in (2.7), for use with the velocity inversion equation (2.5). 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 (2.5)–(2.7) 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 (2.5) – 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-hand $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-hand $1/3$ of the domain, the reconstructed velocity diverges from the true model velocity due to numerical errors arising in the computation of $\unicode[STIX]{x1D63C}$ and $\boldsymbol{b}$ and in the integration of the matrix equation (2.5) 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 $6\times 6$ square kernel (equivalent to a two-dimensional running average over ${\sim}21~\text{km}^{2}$ ), allows for reasonably accurate velocity reconstruction further into the left-hand $1/3$ of the domain (figure 8).

Figure 7.

(b,d,f) Comparison between true and reconstructed $u$ and $v$ velocity components along the three representative tracer iso-contours in an idealized model of the North Atlantic Ocean. (a,c,e) Tracer distribution with the corresponding iso-contours.

Figure 8.

Same as in figure 7 but using a spatially averaged tracer field.

Figure 9.

Velocity reconstruction along the contour shown in figure 7(c,d) (same legend for curves), but with sparse/noisy tracer data. (a) 20 times sparse tracer data (resolution changed from 3.5 to 70 km) in each direction; (b) zero-mean $\tanh \;(\text{Gaussian})$ noise with Gaussian standard deviation 0.01 and amplitude 10 % of the local temporal mean of the tracer added to the tracer data; (c) as (b) but with 20 % noise amplitude; (d) combination of (a) and (c).

Figure 10.

Contour-averaged error as a percentage of the velocity range along the contour shown in figure 7(c,d) as a function of spatial resolution. The square markers represent the case in figure 9(a).

4.2 The effect of resolution and uncertainty in the data

While our analytic velocity reconstruction is, in principle, exact under conditions described in § 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 9(a) 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 figure 7(c,d). An example of the velocity reconstruction using an under-resolved tracer field is shown in figure 9(a).

Our analysis suggests 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 the 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 to be dependent on the kinetic energy spectrum and the dominant length scales of the flow.

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

(4.1) $$\begin{eqnarray}N=A\tanh ({\mathcal{N}}(0,\unicode[STIX]{x1D70E})),\end{eqnarray}$$

where $A$ is the noise amplitude and ${\mathcal{N}}(0,\unicode[STIX]{x1D70E})$ is a zero-mean Gaussian random variable with standard deviation $\unicode[STIX]{x1D70E}$ (generated using MATLAB’s ‘randn’ function), to the model tracer field. Unlike a pure Gaussian signal, (4.1) 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 h) centred 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. (2.7)). 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 optimized for each case.

Figure 9(b) shows the velocity reconstruction along the tracer contour in figure 7(c,d) 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 h period, which is effectively the observation time for this example) and setting a standard deviation $\unicode[STIX]{x1D70E}$ equal to $0.01$ in the definition (4.1). 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, Meissner & Wentz Reference Gentemann, Meissner and Wentz2010). Aerial measurements of dye from multispectral and hyperspectral cameras have similar noise levels (about 20 %) for Rhodamine WT concentrations of about 20 p.p.b. (Clark et al. Reference Clark, Lenain, Feddersen, Boss and Guza2014; Hally-Rosendahl et al. Reference Hally-Rosendahl, Feddersen, Clark and Guza2015). 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 (Clark et al. Reference Clark, Lenain, Feddersen, Boss and Guza2014; Hally-Rosendahl et al. Reference Hally-Rosendahl, Feddersen, Clark and Guza2015) and in the coastal ocean (Rypina et al. Reference Rypina, Kirincich, Lentz and Sundermeyer2016). Figure 9(c) shows the reconstructed velocity with 20 % noise amplitude and the rest of the parameters being the same as figure 9(b). 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 9(d). 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 (2.7) would usually make our reconstruction method susceptible to noise in the data. In the example presented above, however, we find that $\mathit{O}((\unicode[STIX]{x2202}^{2}c/\unicode[STIX]{x2202}t\unicode[STIX]{x2202}x),(\unicode[STIX]{x2202}^{2}c/\unicode[STIX]{x2202}t\unicode[STIX]{x2202}y))\gg \mathit{O}(\unicode[STIX]{x1D707}(\unicode[STIX]{x2202}\unicode[STIX]{x1D6FB}^{2}c/\unicode[STIX]{x2202}x),\unicode[STIX]{x1D707}(\unicode[STIX]{x2202}\unicode[STIX]{x1D6FB}^{2}c/\unicode[STIX]{x2202}y))$ holds throughout the domain, so the right-hand side of (2.7) 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.1) is defined as

(4.2) $$\begin{eqnarray}Pe={\displaystyle \frac{UL}{\unicode[STIX]{x1D707}}}.\end{eqnarray}$$

Using the empirical relation between length scale $L$ [cm] and diffusivity $\unicode[STIX]{x1D707}$ [ $\text{cm}^{2}~\text{s}^{-1}$ ] obtained through synthesis from various dye release experiments in the ocean (Okubo Reference Okubo1971), we can let $\unicode[STIX]{x1D707}=0.0103L^{1.15}$ . Therefore, for typical velocity values, $U=10$ – $100~\text{cm}~\text{s}^{-1}$ , the Péclet number estimated for length scales, $L=0.5$ –10 $^{4}$  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 (2.7) small. Thus, in our velocity reconstruction procedure, the sensitivity to noise would enter largely through the second derivatives in (2.7) or (2.11). 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 figure 9(b–d), the major source of error is the first-order temporal scalar gradient, $\unicode[STIX]{x2202}c/\unicode[STIX]{x2202}t$ . Using noisy scalar data to evaluate purely spatial derivative terms of all orders in (2.7) 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 $\unicode[STIX]{x1D70E}$ in (4.1)). 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 (4.1). 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 Oceans, 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.