# Iterative Methods for Navier–Stokes Inverse Problems

Liam O’Connor1,2, Daniel Lecoanet1,2, Evan H. Anders3, Kyle C. Augustson1,2,
Keaton J. Burns4, Geoffrey M. Vasil5, Jeffrey S. Oishi6, Benjamin P. Brown7
1Department of Engineering Sciences and Applied Mathematics, Northwestern
University, Evanston, IL, 60208 USA  2Center for Interdisciplinary Exploration
and Research in Astrophysics, Northwestern University, Evanston, IL, 60201 USA
3Kavli Institute for Theoretical Physics, University of California Santa
Barbara, Santa Barbara, CA, 93106 USA  4Department of Mathematics,
Massachusetts Institute of Technology, Cambridge, MA, 02142 USA  5School of
Mathematics, Edinburgh University, EH9 3FD, UK  6Department of Physics and
Astronomy, Bates College, Lewiston, ME, 04240 USA  7Department of
Astrophysical and Planetary Sciences, University of Colorado Boulder, Boulder,
CO, 80309 USA

###### Abstract

Even when the partial differential equation underlying a physical process can
be evolved forward in time, the retrospective (backward in time) inverse
problem often has its own challenges and applications. Direct Adjoint Looping
(DAL) is the defacto approach for solving retrospective inverse problems, but
it has not been applied to deterministic retrospective Navier–Stokes inverse
problems in 2D or 3D. In this paper, we demonstrate that DAL is ill-suited for
solving retrospective 2D Navier–Stokes inverse problems. Alongside DAL, we
study two other iterative methods: Simple Backward Integration (SBI) and the
Quasi-Reversible Method (QRM). Our iterative SBI approach is novel while
iterative QRM has previously been used. Using these three iterative methods,
we solve two retrospective inverse problems: 1D Korteweg–de Vries–Burgers
(decaying nonlinear wave) and 2D Navier–Stokes (unstratified Kelvin–Helmholtz
vortex). In both cases, SBI and QRM reproduce the target final states more
accurately and in fewer iterations than DAL. We attribute this performance gap
to additional terms present in SBI and QRM’s respective backward integrations
which are absent in DAL.

††preprint: APS/123-QED

##  I Introduction

Consider a retrospective inverse problem where the variable of interest
$\overline{X}$ is constrained by a partial differential equation (PDE)
involving space $x$ and time $t$:

$\displaystyle\mathcal{F}[\overline{X}(x,t)]=0.$ (1)

We aim to determine the unique initial condition $\overline{X}(x,0)$
corresponding to a given final state $\overline{X}(x,t_{f})$ where the final
time $t_{f}$ is known. These “retrospective inverse problems” are often ill-
conditioned or ill-posed, even when the forward problem is readily solvable.

Despite potential difficulties, retrospective inverse problems arise in a wide
variety of contexts. Scientists and engineers have documented numerous cases
where retrospective inversion is feasible and practical. [1] inverted chemical
reaction fronts governed by a 1D reaction-diffusion equation. [2] inverts
tumor growth using biophysically-motivated constraints on the initial
condition. [3] inverted supersonic supernova expansions governed by the 1D
compressible Navier–Stokes equation. [4, 5] examined the history of buoyancy-
driven flows in Earth’s mantle by approximating palaeo temperature
distributions as well as variable parameters such as the thermal diffusivity.

Although there are problem-specific details in their methodologies, [1, 2, 3,
4, 5] use the same optimal control technique known as Direct Adjoint Looping
(DAL). Using DAL, they minimize the cost functional

$\displaystyle\mathcal{J}_{f}^{X}$
$\displaystyle\equiv\frac{1}{2}\Big{\langle}\Big{|}X^{\prime}(x,t_{f})\Big{|}^{2}\Big{\rangle}=\frac{1}{2}\Big{\langle}\Big{|}X(x,t_{f})-\overline{X}(x,t_{f})\Big{|}^{2}\Big{\rangle},$
(2)

where the angled brackets $\langle\cdot\rangle$ denote a spatial integral over
the problem domain and the trial solution $X$ has some deviation
$X^{\prime}\equiv X-\overline{X}$. They construct sequences of trial
solutions, such as $\\{X_{n}(x,t)\\}$, where each element $X_{n}$ is computed
from a corresponding trial initial condition $X_{n}(x,0)$. Convergence to the
target initial condition $\overline{X}(x,0)$ is seldom guaranteed because the
user provides an initial guess $X_{0}(x,0)$.

Related 2D Navier–Stokes inverse problems have been studied analytically from
several perspectives. [6] uses boundary conditions as control inputs (with a
given initial condition) and verifies controllability of certain 2D fluid
systems. Later, the existence of solutions to a more demanding retrospective
Navier–Stokes inverse problem with final overdetermination was established by
[7]. In [7], the authors do not attempt to solve their inverse problem.
Instead they provide an analysis of the problem where a separable source term
(as well as the initial condition) are left to be determined. Following this,
the data-completion problem of determining an unknown viscosity was studied
analytically by [8].

More recently, researchers extended these analytical studies of Navier–Stokes
inverse problems by implementing data assimilation algorithms. [9] studies an
inverse problem governed by the parabolized Navier–Stokes (PNS) equation.
Following this, [10] studies data assimilation algorithms for non-
deterministic Navier–Stokes inverse problems. They demonstrate several
advantages of using filters for numerical weather prediction. Yet more
recently, [11, 12, 13, 14] used neural networks to solve analogous
Navier–Stokes data-completion problems.  In particular, [14] compares DAL with
a physics-informed neural network approach while solving inverse problems
governed by the viscous Burgers equation as well as 2D incompressible
Navier–Stokes.

In this paper, we will also compare the performance of several iterative
methods (such as DAL) by applying them to a retrospective Navier–Stokes
inverse problem. Various Navier–Stokes inverse problems have been studied
analytically and numerically using optimal control. However, we will
demonstrate that conventional optimal control methods are ill-suited for this
class of inverse problems, even when a numerically-consistent dataset is
available.

With DAL, we can minimize an arbitrary cost functional such as
$\mathcal{J}_{f}^{X}$. First, we must solve the forward problem on $t:0\to
t_{f}$. This forward problem is constrained by some PDE (such as
Navier–Stokes) which we express generally as eqn 1. At this point, we can
evaluate $\mathcal{J}_{f}^{X}$ and compare the trial ($X(x,t_{f})$) and target
($\overline{X}(x,t_{f})$) final states. Information from time $t=t_{f}$ is
then propagated backward, to time $t=0$, by solving the linear adjoint problem
$\mathcal{F}^{\dagger}[\chi,X]=0$ on $t:t_{f}\to 0$, where $\chi(x,t)$ denotes
the adjoint variable. Backward integration of the adjoint equations allows us
to compute the gradient or functional derivative of $\mathcal{J}_{f}^{X}$ with
respect to $X(x,0)$. Using this gradient, we adjust our trial initial
condition by some amplitude $\gamma$:
$X_{n+1}(x,0)=X_{n}(x,0)+\gamma\chi_{n}(x,0)$.

In practice, researchers often customize the DAL algorithm to improve
performance in specific cases. For example, [15] used the conjugate-gradient
method with a discrete adjoint formulation to optimize initial conditions
subject to norm constraints. For the PNS inverse problem, [9] combined
Principle Orthogonal Decomposition (POD) with four-dimensional variational
assimilation (4DVAR), an optimal control technique analogous to DAL. In both
examples, researchers reduced the computational cost of each DAL iteration by
reformulating the adjoint while simultaneously reducing the required number of
DAL iterations by implementing quasi-Newton minimization routines. For this
investigation, we take a different approach. Instead of reformulating the
adjoint or changing the way in which gradient information is used, we will
demonstrate that DAL’s shortcomings in the context of Navier–Stokes inversion
can be mitigated by appending the adjoint system with additional advective
terms.

The additional advective terms not present in DAL are present in two existing
nonlinear backward integration methods: Simple Backward Integration (SBI) and
the Quasi-Reversible Method (QRM). SBI and QRM are designed to invert the
target final state by approximating the (ill-posed) integration of the
original equations backward in time. SBI accomplishes this by reversing the
sign of the problematic diffusion term, i.e.
$a\partial_{x}^{2}X\to-a\partial^{2}_{x}X$. In contrast, QRM is a
regularization framework which preserves the ill-posed diffusion term [16].
With QRM, we append a small hyperdiffusion term to the otherwise ill-posed
constraint equation, i.e. $a\partial_{x}^{2}X\to
a(\partial^{2}_{x}+\varepsilon\partial^{4}_{x})X$ where $0<\varepsilon\ll 1$.

We introduce these concepts as backward integration methods rather than
iterative methods because their previous implementations were rarely iterated.
SBI (which has never been iterated) was introduced by [17] to infer the
history of Earth’s convective mantle in the viscinity of the southern African
superswell. Following this, [4] and [5] used SBI to construct initial guesses
$X_{0}(x,0)$, which they refined iteratively using DAL. SBI was also examined
by [18] to dismiss a time-reversibility hypothesis on Subgrid-Scale (SGS)
Large-Eddy Simulations (LES).

[19] applied iterative QRM to ill-posed linear inverse problems and derives
rigorous statements of convergence. More recently, [20] used iterative QRM to
compute the unknown source in a nonlinear hyperbolic equation. Using Carleman
estimates, they demonstrate QRM’s convergence to the correct source at an
exponential rate. While an analysis involving Carleman estimates is outside
the scope of this paper, we refer to [21].

We compare DAL, SBI, and QRM by solving two retrospective inverse problems,
respectively constrained by Korteweg–de Vries–Burgers (KdVB) and
incompressible Navier–Stokes. The novelty of our work lies in the iterative
application of SBI and QRM. At the start of iteration $n$, we evolve some
initial trial state $X_{n}(x,0)$ to its final state $X_{n}(x,t_{f})$. Then, we
initialize SBI and QRM using the final deviation $X^{\prime}_{n}(x,t_{f})$ in
order to approximate the initial deviation $X^{\prime}_{n}(x,0)$. Finally,
assuming we know $X^{\prime}_{n}(x,0)$, we refine the trial state as
$X_{n+1}(x,0):=X_{n}(x,0)-X^{\prime}_{n}(x,0)$ to conclude the iteration.
Equivalently, iterative SBI and QRM minimize the inaccessible objective
functional

$\displaystyle\mathcal{J}_{0}^{X}$
$\displaystyle\equiv\frac{1}{2}\Big{\langle}\Big{|}X^{\prime}(x,0)\Big{|}^{2}\Big{\rangle}=\frac{1}{2}\Big{\langle}\Big{|}X(x,0)-\overline{X}(x,0)\Big{|}^{2}\Big{\rangle}$
(3)

by repeatedly approximating its gradient $X^{\prime}(x,0)$.
$\mathcal{J}_{0}^{X}$ measures the trial state’s error at the initial time,
whereas $\mathcal{J}_{f}^{X}$ does so at the final time.

The remainder of this paper is organized as follows: In Section II we provide
a conceptual description of our iterative methods (DAL, SBI, and QRM). In
Section III we illustrate these methods by applying them to a retrospective
KdVB inverse problem. In Section IV we generalize these methods for multiple
spatial dimensions, and compare their implementations using a 2D
incompressible Navier–Stokes inverse problem. In Section V, we demonstrate
that the adjoint equation which yields the gradient of $\mathcal{J}_{0}^{X}$
is ill-posed, nonlinear, and equivalent to the perturbed constraint eqn 1. In
Section VI we summarize our results and discuss future applications for
iterative SBI and QRM.

##  II Overview of Methods

The Korteweg–de Vries–Burgers (KdVB) equation

$\displaystyle\mathcal{F}[\overline{u}(x,t)]\equiv\partial_{t}\overline{u}+\overline{u}\partial_{x}\overline{u}-a\partial_{x}^{2}\overline{u}+b\partial_{x}^{3}\overline{u}=0$
(4)

provides a concise context for demonstrating each method (DAL, SBI, QRM). Let
$\overline{u}(x,t)$ denote to the target solution where $0\leq t\leq t_{f}$ in
a 1D periodic domain $0\leq x<L_{x}$. The retrospective inverse problem is
solved by recovering the target initial state $\overline{u}(x,0)$ where the
target final state $\overline{u}(x,t_{f})$ is given. The general strategy has
us construct trial solutions $u(x,t)$ satisfying eqn 4, which we compute from
trial initial states $u(x,0)$. We measure a trial state’s proximity from the
target using the cost functionals $\mathcal{J}^{u}_{0}$ and
$\mathcal{J}^{u}_{f}$, which are defined for an arbitrary variable $X$ in eqns
2 and 3. Both $\mathcal{J}^{u}_{0}$ and $\mathcal{J}^{u}_{f}$ vary with
respect to the trial initial state $u(x,0)$. However, given that our knowledge
of the target solution is limited to the final time $t_{f}$, it follows that
we can only evaluate $\mathcal{J}^{u}_{0}$ if the target solution is known
apriori.

Figure 1: Conceptual illustration of iterative methods. Purple: the target
solution $\overline{u}(x,t)$ is known only at the final time $t_{f}$. Blue: we
evolve a trial state $u(x,0)$ to $u(x,t_{f})$ where we then initialize one of
three backward integration systems using a single compatibility condition.
Each backward integration gives a solution at $t=0$ which we use to update our
trial initial condition. We seek the trial solution’s deviation
$u^{\prime}=u-\overline{u}$ at the initial time $t=0$. Yellow: Direct Adjoint
Looping (DAL) is a linear backward integration. For advective systems such as
KdVB and Navier–Stokes, DAL advects the final deviation $u^{\prime}(x,t_{f})$
by the trial solution $u$. Orange: Simple Backward Integration (SBI) is a
hybrid method, which supplements the linear DAL system with two additional
advective terms. Red: the Quasi-Reversible Method (QRM) aims to compute the
deviation by introducing a small hyperdiffusion term with coefficient
$\varepsilon$. The QRM backward integration is ill-posed when $\varepsilon=0$
and numerically unstable as $\varepsilon\to 0$.

###  II.1 Direct Adjoint Looping (DAL)

In general, we can cast an inverse problem into the form of an optimization
problem by defining and minimizing an appropriate cost functional [22]. DAL
provides an efficient avenue for computing a cost functional’s gradient. For
KdVB, we minimize $\mathcal{J}^{u}_{f}$ by computing its gradient with respect
to $u(x,0)$. Once this gradient is known, we refine the trial initial state
$u(x,0)$ by following the direction of steepest descent.

Using DAL, the gradient is computed as follows. We define the Lagrangian
$\mathcal{L}$

$\displaystyle\mathcal{L}$
$\displaystyle\equiv\int_{0}^{t_{f}}\Big{\langle}\mu(x,t)\cdot\mathcal{F}[u(x,t)]\Big{\rangle}dt+\mathcal{J}^{u}_{f}$
(5)

where $\mu$ is the adjoint variable corresponding to $u$. $\mathcal{L}$’s
variations with respect to $u$ and $\mu$ must disappear at the optimum
$\overline{u}$, yielding

$\displaystyle 0<t<t_{f}:\quad$ $\displaystyle
0=\partial_{t}u+u\partial_{x}u-a\partial_{x}^{2}u+b\partial_{x}^{3}u$ (6)
$\displaystyle t=t_{f}:\quad$ $\displaystyle\mu=-u^{\prime}$ (7)
$\displaystyle 0<t<t_{f}:\quad$ $\displaystyle
0=\partial_{t}\mu+u\partial_{x}\mu+a\partial_{x}^{2}\mu+b\partial_{x}^{3}\mu.$
(8)

The gradient of $\mathcal{J}^{u}_{f}$ resides in the manifold of solutions
satisfying eqns 6, 7, and 8. Figure 1 illustrates how these equations are
solved sequentially. Given a trial initial condition, we start by solving the
forward problem eqn 6 from $t:0\to t_{f}$ (shown in blue). Next, we initialize
the adjoint variable using the compatibility condition eqn 7. Finally, we
solve the adjoint eqn 8 backward in time $t:t_{f}\to 0$ (shown in yellow),
yielding the desired gradient $-\mu(x,0)$.  The adjoint involves the trial
solution $u$, so we must store data from the forward solve in memory. In
practice this leads to redundant computations (due to memory constraints),
offering another area for tangible improvement.

Once $\mathcal{J}^{u}_{f}$’s gradient is known, we use gradient descent with
the popular Barzilai–Borwein method [23] to select step sizes. In addition to
gradient descent, we also minimize $\mathcal{J}^{u}_{f}$ using the second-
order sparse L-BFGS algorithm [24] which utilizes the gradient’s history over
several iterations.  We select these two optimization routines for comparison
after observing that they converge faster than alternatives. 222In addition to
the results shown in section III, we perform numerical experiments on the KdVB
inverse problem using conjugate-gradient and line search routines (provided by
scipy version 1.12.0). In both cases we use the default parameters. For the
line search, we use the negative gradient for the descent direction.

###  II.2 SBI and QRM Setup

The KdVB retrospective inverse problem can be solved by finding the initial
deviation $u^{\prime}(x,0)$ of some trial state $u(x,0)$.
$u^{\prime}=u-\overline{u}$ obeys a perturbation equation

$\displaystyle
0=\partial_{t}u^{\prime}-a\partial_{x}^{2}u^{\prime}+b\partial_{x}^{3}u^{\prime}+u\partial_{x}u^{\prime}+u^{\prime}\partial_{x}u-u^{\prime}\partial_{x}u^{\prime}$
(9)

which is ill-posed over $t:t_{f}\to 0$. We cannot compute
$u^{\prime}(x,t<t_{f})$ numerically due to the intrinsic loss of information
in diffusive systems. DAL circumvents this obstacle by computing a gradient,
whereas SBI and QRM approximate the deviation $u^{\prime}$. Let $\tilde{\mu}$
denote our SBI and QRM approximations for $-u^{\prime}$.

###  II.3 Simple Backward Integration (SBI)

SBI is one method for approximating the ill-posed backward integration, in
which we reverse the sign of the problematic diffusion term, i.e.
$a\partial_{x}^{2}\tilde{\mu}\to-a\partial_{x}^{2}\tilde{\mu}$ such that

$0=\partial_{t}\tilde{\mu}+a\partial_{x}^{2}\tilde{\mu}+b\partial_{x}^{3}\tilde{\mu}+u\partial_{x}\tilde{\mu}+\tilde{\mu}\partial_{x}u+\tilde{\mu}\partial_{x}\tilde{\mu}.$
(10)

This backward integration is well-posed. Eqn 10 contains every term in eqn 8
(DAL) with two additional advective terms: $\tilde{\mu}\partial_{x}u$ and
$\tilde{\mu}\partial_{x}\tilde{\mu}$.

For the KdVB inverse problem, we first initialize eqn 10 with
$\tilde{\mu}(x,t_{f})=-u^{\prime}(x,t_{f})$, just as with with DAL. We then
solve eqn 10 on $t:t_{f}\to 0$. Finally, SBI yields an approximation
$\tilde{\mu}(x,0)\approx-u^{\prime}(x,0)$. The SBI backward integration loop
is illustrated in orange in Figure 1, where the additional SBI terms cause
$\tilde{\mu}$ to deviate from the linear DAL calculation of $\mu$. After
completing an SBI loop, we update the trial state
$u_{n+1}(x,0)=u_{n}(x,0)+\tilde{\mu}_{n}(x,0)$. This is analogous to
performing gradient descent with a fixed step size of unity.333SBI and QRM do
not generally return any sensible gradient, so the specified updates are not
technically gradient descent. Attempting a constant step size of
$\gamma\sim\mathcal{O}(1)$ with DAL is unstable for the direct minimization of
$\mathcal{J}^{u}_{f}$. SBI can often be accelerated using a larger step size.

In previous investigations, SBI was applied once to the target solution. Our
iterative implementation differs, because we apply SBI to the deviation
equation. Despite the inaccuracies introduced by reversing the sign of the
diffusive term, we will demonstrate that iterative SBI is an effective
algorithm. SBI has already been implemented to develop initial guesses for
retrospective inverse problems. We extend this work by comparing SBI with DAL
rather than using them separately for different purposes.

###  II.4 Quasi-Reversible Method (QRM)

QRM deals with the ill-posed diffusion term by introducing a small, well-posed
hyperdiffusion term: $a\partial_{x}^{2}\tilde{\mu}\to
a(\partial_{x}^{2}\tilde{\mu}+\varepsilon\partial_{x}^{4}\tilde{\mu})$. This
method replaces eqn 9 by

$0=\partial_{t}\tilde{\mu}-a(\partial_{x}^{2}\tilde{\mu}+\varepsilon\partial_{x}^{4}\tilde{\mu})+b\partial_{x}^{3}\tilde{\mu}+u\partial_{x}\tilde{\mu}+\tilde{\mu}\partial_{x}u+\tilde{\mu}\partial_{x}\tilde{\mu}.$
(11)

$0<\varepsilon\ll 1$ is a free parameter chosen according to the numerical
resolution as well as the solutions’ spatial scales. Notice how the ill-posed
term remains, but the equation becomes well-posed because small-scale modes do
not grow arbitrarily fast.

The QRM loop is carried out using the same procedure as SBI. We initialize eqn
11 with $\tilde{\mu}(x,t_{f})=-u^{\prime}(x,t_{f})$ and update the trial state
$u_{n+1}(x,0)=u_{n}(x,0)+\tilde{\mu}_{n}(x,0)$. In Figure 1, we illustrate how
QRM (red) recovers the unknown deviation precisely as $\varepsilon\to 0$.
However this limit is numerically unstable.

##  III KdVB Numerical Experiments

We implement each iterative method (DAL, SBI, QRM) to solve an inverse problem
constrained by KdVB (eqn 4). The target solution $\overline{u}(x,t)$ (shown
left in Figure 2) consists of a dissipating nonlinear wave where $a=0.02$ and
$b=0.04$. The spatial domain is 1D periodic $x\in[0,2\pi]$ and we solve eqn 4
over $t:0\to 3\pi$. The target initial condition

$\displaystyle\overline{u}(x,0)$
$\displaystyle=3\cosh^{-2}\left(\frac{x-\pi}{2\sqrt{b}}\right),$ (12)

would propagate rightwards as a stable soliton with a constant speed of one if
diffusion were not present (i.e. if $a=0$ in eqn 4). This retrospective
problem is nontrivial, as the final state cannot be evolved backwards in time
via conventional timestepping algorithms.

Figure 2: Target KdVB solution $\overline{u}(x,t)$ (left) plotted as a
function of space $x$ and time $t$. We aim to recover the initial condition of
a nonlinear wave which propagates right at an initial speed of one. The
problem’s spatial domain is 1D periodic. We include diffusion ($a=0.04$) such
that its speed decreases as the wave propagates. Beginning at $t=3\pi$, SBI
(center) evolves $\overline{u}(x,t_{f})$ backward to $t=0$ by reversing the
sign of the diffusive term. The wave diffuses during the forward and backward
integrations, such that its amplitude and speed of propagation decrease in
both directions. We also perform the QRM backward integration (right) on
$t:3\pi\to 0$. Unlike SBI, QRM repopulates small-scale modes which are lost
during forward integration. As the wave’s peak regains its former amplitude,
its speed of propagation increases mimicking $\overline{u}(x,t)$.

The target and trial solutions along with each backward integration are
carried out using the Dedalus open-source pseudospectral python framework
[27]. We represent $u,\;\mu,$ and $\tilde{\mu}$ as vectors of 128 real Fourier
modes. While computing $u(x,t)$, we store the solution vector at each
timestep. Linear (diffusive and dispersive) terms are treated implicitly
whereas the nonlinear advective term is evaluated explicitly by multiplying
$u$’s grid data on a 3/2 dealiased grid. The same holds when we perform
backward integrations, where terms such as $\mu\partial_{x}u$ must be treated
explicitly to avoid dense matrix operations at every timestep.  Compared to
DAL, the SBI and QRM backward integration equations contain more of these
terms. Consequently, each SBI and QRM iteration requires approximately $30\%$
more computation time. For every solve we use a 2nd-order Runge–Kutta scheme
with a fixed timestep $\Delta t=0.01$.  We verify the adjoint solver’s
accuracy by comparing its output to a gradient obtained via finite-difference.
Our code for this paper is available at https://github.com/liamoconnor9/adjop

We initialize each iterative method using an initial guess $u^{0}(x,0)=0$.
(See Appendix A for the same comparison beginning with an SBI initial guess.)
The first SBI iteration is shown in the center panel of Figure 2. Diffusion of
the nonlinear wave causes its amplitude and speed of propagation to decrease
such that the wave’s peak does not return to its original position. The first
QRM iteration (shown right in Figure 2) reconstructs the wave by inverting the
advection alongside diffusion. In this case, the wave’s peak returns to its
approximate initial position. For QRM we use $\varepsilon=0.01$ whereas SBI
involves no free parameter.

Figure 3: Evolution of the KdVB trial initial conditions $u_{n}(x,0)$ using an
initial guess $u_{0}(x,0)=0$ (shown in green). We compare two DAL runs (top
row) with SBI (bottom left) and QRM (bottom right). (Top left) gradient
descent approaches the target initial state gradually. (Top right) DAL
minimization of $\mathcal{J}^{u}_{f}$ using L-BFGS approaches the target state
more rapidly, but this method introduces undesired medium-scale features. SBI
captures the target state faster and more directly than DAL. QRM resembles the
target state immediately $(n=1)$. However, this method introduces low-
amplitude small-scale errors which do not subside at $n=200$.

The top row of Figure 3 illustrates the evolution of DAL trial initial states
$u_{n}(x,0)$. We minimize $\mathcal{J}^{u}_{f}$ using gradient descent (top
left) and L-BFGS (top right). For advective systems such as KdVB, the
amplitude of each feature is proportional to its speed of propagation.
Consequently, gradient descent gives a sequence of initial states which
increase in amplitude while shifting leftwards. $u_{10}(x,0)$ contains a
shallow peak centered near $x=\pi/2$. $u_{50}(x,0)$ is centered near $x=0$
with approximately twice the amplitude of $u_{10}(x,0)$. Similarly,
$u_{200}(x,0)$ is centered near $x=3\pi/2$ with approximately triple the
amplitude of $u_{10}(x,0)$. Their respective final states ($u_{10}(x,3\pi)$,
$u_{50}(x,3\pi)$, and $u_{200}(x,3\pi)$) all contain peaks which approximately
align with the target’s ($\overline{u}(x,3\pi)$).

Rather than following the direction of steepest descent, L-BFGS anticipates
changes in the gradient by constructing a sparse Hessian representation.
Although this quasi-Newton method recovers the target initial state better
than gradient descent at $n=200$, L-BFGS still does not yield a quality
approximation of $\overline{u}(x,0)$. This case introduces additional wave
features in $u_{10}(x,0)$, $u_{50}(x,0)$, and $u_{200}(x,0)$ which do not
appear in the target.

To locate the target state, we introduce additional advective terms in the
backward integration. The bottom row of Figure 3 illustrates the evolution of
trial initial states $u_{n}(x,0)$ refined via SBI (bottom left) and QRM
(bottom right) over 200 iterations. SBI does not effectively approximate the
initial deviation $u^{\prime}_{n}(x,0)$ (as shown for $n=0$ in Figure 2), nor
does it follow the direction of steepest descent of $\mathcal{J}^{u}_{f}$.
This algorithm smoothly guides the initial states’ peaks directly toward the
target. In this case, $u_{50}(x,0)$ is indistinguishable from
$\overline{u}(x,0)$. QRM loosely captures the target’s structure in a single
iteration by approximating the initial deviation $u^{\prime}_{n}(x,0)$. In
this case, $u_{10}(x,0)$ is nearly indistinguishable from $\overline{u}(x,0)$.

After one iteration, the SBI and QRM initial states’ peaks are centered
rightwards of the target feature (see $u_{1}(x,0)$ in Figure 3, bottom left
and right panels). These nonlinear backward integrations couple the
deviation’s amplitude with the resulting direction of our update for
$u_{n}(x,0)$. SBI refines each trial initial state by advecting the final
deviation backward in time as it continues to diffuse. This decreases the
speed of backward advection such that each update does not return to the
position of the initial deviation. This is illustrated by $u_{10}(x,0)$ and
$u_{50}(x,0)$ in the bottom left panel of Figure 3. Their peaks shift
leftward, approaching the target with each iteration. When the trial final
state intercepts the target, the deviation shrinks and SBI’s backward
diffusion has less influence on the direction of each subsequent update. QRM
refines each trial state by advecting the final deviation while simultaneously
inverting diffusion. The hyperdiffusion term slows backward advection to a
lesser extent than SBI, yielding smaller rightward misalignments. These
updates approximate the initial deviation as it shrinks with each iteration,
such that iterative QRM behaves like a one-shot method (e.g. [28]).

For QRM and SBI, we plot the deviation
$u^{\prime}(x,0)=u_{200}(x,0)-\overline{u}(x,0)$ at iteration $n=200$ in
Figure 4. Both curves have relatively large features coinciding with the
target initial condition’s peak at $x=\pi$. The SBI deviation has a larger
maximum magnitude while the QRM deviation oscillates at a particular
wavenumber. When we increase $\varepsilon$, these oscillations subside, but it
requires more iterations to obtain a comparable result. We select
$\varepsilon=0.01$ to illustrate the presence of these oscillations, while
still achieving the desired outcome in a small number of iterations.

Figure 4: Deviation of SBI and QRM trial initial conditions, plotted as a
function of $x$ at iteration $n=200$. Both curves have prominent peaks which
overlap near $x=\pi$. The SBI deviation has a slightly larger magnitude while
the QRM deviation oscillates at a particular wavenumber. Figure 5: Initial
errors ($\mathcal{J}^{u}_{0}$, top) and final errors ($\mathcal{J}^{u}_{f}$,
bottom) for the 1D KdVB inverse problem plotted as a function of iteration
$n$. We compare DAL, SBI, and QRM. DAL performs slightly better when paired
with the second-order sparse minimization routine L-BFGS (LB
$\mathcal{J}^{u}_{f}$), whereas gradient descent (GD $\mathcal{J}^{u}_{f}$)
performs the worst. Our iterative applications of SBI and QRM minimize
$\mathcal{J}^{u}_{0}$ and $\mathcal{J}^{u}_{f}$ more effectively than DAL. QRM
minimizes both $\mathcal{J}^{u}_{0}$ and $\mathcal{J}^{u}_{f}$ more rapidly
than SBI.

Errors from each optimization study are plotted in Figure 5. DAL does not
effectively minimize $\mathcal{J}^{u}_{0}$ or $\mathcal{J}^{u}_{f}$. In
similar investigations, L-BFGS has received significant attention thanks to
its DAL-compatible features [5]. It converges faster than gradient descent
which is also reflected in our experiment. However, the SBI and QRM iterative
methods both outperform DAL by several orders of magnitude. QRM minimizes both
$\mathcal{J}^{u}_{0}$ and $\mathcal{J}^{u}_{f}$ more rapidly than SBI, as
shown in Figure 5. At $n\approx 100$, SBI stalls and QRM approaches floating-
point precision in terms of $\mathcal{J}^{u}_{f}$. We attribute this gap to an
energy discrepancy between the SBI and QRM updates. For QRM applied to a
highly diffusive system, the energy $\langle|\tilde{\mu}(x,t)|^{2}\rangle$
increases during backward integration (provided $\varepsilon$ is sufficiently
small) such that
$\langle|\tilde{\mu}(x,0)|^{2}\rangle\sim\langle|u^{\prime}(x,0)|^{2}\rangle$.
This is illustrated in the first QRM iteration (Figure 2, right). In contrast,
SBI diffuses during backward integration (Figure 2, center) such that
$\langle|\tilde{\mu}(x,0)|^{2}\rangle<\langle|u^{\prime}(x,0)|^{2}\rangle$. We
can eliminate this energy discrepancy and often accelerate SBI by increasing
the step sizes of the updates (e.g.
$u_{n+1}(x,0)=u_{n}(x,0)+1.1\mu_{n}(x,0)$).

##  IV 2D Navier–Stokes Numerical Experiments

Consider 2D flow with spatial coordinate $\bm{x}=x\bm{\hat{x}}+y\bm{\hat{y}}$.
The velocity $\bm{\overline{u}}(\bm{x},t)$ obeys the incompressible
Navier–Stokes equation

$\displaystyle\partial_{t}\bm{\overline{u}}+\bm{\overline{u}}\cdot\nabla\bm{\overline{u}}+\nabla\overline{p}=\nu\nabla^{2}\bm{\overline{u}}\quad\textrm{and}\quad\nabla\cdot\bm{\overline{u}}=0$
(13)

where $\overline{p}$ is the pressure and $\rm{Re}\equiv\nu^{-1}=50,000$. The
domain is doubly-periodic with $x\in[0,1]$ and $y\in[-1,1]$, such that its
boundary conditions are given by
$\bm{\overline{u}}|_{x=0}=\bm{\overline{u}}|_{x=1}$ and
$\bm{\overline{u}}|_{y=-1}=\bm{\overline{u}}|_{y=1}$.

For our inverse problem, we use divergence-cleaning to construct the target
initial condition $\bm{\overline{u}}(\bm{x},0)$. First, consider the following
velocity field which has a nonzero divergence

$\displaystyle\bm{\overline{v}}(\bm{x})\cdot\bm{\hat{x}}$
$\displaystyle=\frac{1}{2}\bigg{[}\tanh\bigg{(}10\Big{(}y-\frac{1}{2}\Big{)}\bigg{)}-\tanh\bigg{(}10\Big{(}y+\frac{1}{2}\Big{)}\bigg{)}\bigg{]}\;;$
$\displaystyle\bm{\overline{v}}(\bm{x})\cdot\bm{\hat{y}}$
$\displaystyle=\frac{1}{10}\sin\big{(}2\pi
x\big{)}\Bigg{[}\exp\bigg{(}-100\Big{(}y-\frac{1}{2}\Big{)}^{2}\bigg{)}$ (14)
$\displaystyle\qquad\qquad\qquad\qquad-\exp\bigg{(}-100\Big{(}y+\frac{1}{2}\Big{)}^{2}\bigg{)}\Bigg{]}.$

Next, we perform divergence cleaning on $\bm{\overline{v}}$ by solving the
boundary value problem $\nabla^{2}p_{v}+\nabla\cdot\bm{\overline{v}}=0$. The
divergence-free target initial condition is then given by
$\bm{\overline{u}}(\bm{x},0)=\bm{\overline{v}}(\bm{x})+\nabla p_{v}(\bm{x})$.

For our inverse problem, we run a target simulation by evolving
$\bm{\overline{u}}(\bm{x},0)$ according to eqn 13 until $t_{f}=20$. The target
initial condition has two thin layers of positive and negative vorticity
concentrated near $y=1/2$ and $y=-1/2$. Advection causes these vortex sheets
to wind counter-clockwise and clockwise while viscosity smooths its small-
scales.

Given the final velocity field $\bm{\overline{u}}(\bm{x},t_{f})$, we develop
trial solutions $\bm{u}(\bm{x},t)$ (also obeying eqn 13) which are meant to
approach the target solution $\bm{\overline{u}}(\bm{x},t)$. All PDE solves
(except the GD $\mathcal{J}^{\bm{u}}_{f}$ double-resolution study) are
performed on a 128 by 256 grid of real Fourier modes using the Dedalus
pseudospectral framework [27]. We timestep with $\Delta t=0.002$ using a
second-order Runge–Kutta scheme. Nonlinear operations are evaluated on a 3/2
dealiased grid. We plot the upper half the domain (unit square) because the
solutions have an approximate symmetry about $y=0$. The target solution’s
vorticity evolution is illustrated in the top row of Figure 6.

###  IV.1 DAL

Let $\bm{\mu}$ and $\Pi$ be the adjoint variables respectively corresponding
to $\bm{u}$ and $p$. The Navier–Stokes adjoint equations are then given by

$\displaystyle 0$ $\displaystyle=\nabla\cdot\bm{\mu}$ (15)
$\displaystyle\bm{0}$
$\displaystyle=\partial_{t}\bm{\mu}+\bm{u}\cdot\nabla\bm{\mu}\hbox{\pagecolor{white}$\displaystyle-\bm{\mu}\cdot(\nabla\bm{u})^{T}$}+\nabla\Pi+\nu\nabla^{2}\bm{\mu}.$

[29] provides a detailed derivation of eqn 15. Using DAL, we minimize two
objective functionals with associated compatibility conditions

$\displaystyle\mathcal{J}^{\bm{u}}_{f}$
$\displaystyle\equiv\frac{1}{2}\left\langle||\bm{u^{\prime}}(\bm{x},t_{f})||^{2}\right\rangle\rightarrow\bm{\mu}(\bm{x},t_{f})=-\bm{u^{\prime}}(\bm{x},t_{f});$
(16) $\displaystyle\mathcal{J}^{\omega}_{f}$
$\displaystyle\equiv\frac{1}{2}\left\langle\left(\omega^{\prime}(\bm{x},t_{f})\right)^{2}\right\rangle\rightarrow\bm{\mu}(\bm{x},t_{f})=\nabla^{\perp}\omega^{\prime}(\bm{x},t_{f}).$
(17)

The skew-gradient $\nabla^{\perp}\equiv(-\partial_{y},\;\partial_{x})$ and the
vorticity $\omega\equiv\nabla^{\perp}\cdot\bm{u}$. As with KdVB, a prime
denotes the deviation of a trial variable from its target, e.g.
$\omega^{\prime}\equiv\omega-\overline{\omega}$ and so on.
$\mathcal{J}^{\bm{u}}_{f}$ quantifies the trial states’ proximity to the
target at the final time $t_{f}$ in terms of velocity whereas
$\mathcal{J}^{\omega}_{f}$ does so using vorticity. We implement the gradient
descent (GD) and L-BFGS optimization routines to minimize these two cost
functionals.  We verify the accuracy of the adjoint solver by using the
finite-difference approximation to reproduce the gradient at low resolution.

###  IV.2 SBI and QRM

Our iterative procedures for SBI and QRM are analogous to those described in
the previous section. The ill-posed backward integration of $\bm{u^{\prime}}$
obeys the nonlinear perturbation equation

$\displaystyle 0$ $\displaystyle=\nabla\cdot\bm{u^{\prime}}$ (18)
$\displaystyle\bm{0}$
$\displaystyle=\partial_{t}\bm{u^{\prime}}+\bm{u}\cdot\nabla\bm{u^{\prime}}+\bm{u^{\prime}}\cdot\nabla\bm{u}-\bm{u^{\prime}}\cdot\nabla\bm{u^{\prime}}+\nabla
p^{\prime}-\nu\nabla^{2}\bm{u^{\prime}}.$ Next we substitute the
approximations $-\bm{\tilde{\mu}}$ for $\bm{u^{\prime}}$ and $-\tilde{\Pi}$
for $p^{\prime}$. For SBI, we reverse the sign of the viscous term
$\nu\nabla^{2}\bm{\tilde{\mu}}$. Its corresponding approximation for eqn 18 is
given by $\displaystyle 0$ $\displaystyle=\nabla\cdot\bm{\tilde{\mu}}$ (19)
$\displaystyle\bm{0}$
$\displaystyle=\partial_{t}\bm{\tilde{\mu}}+\bm{u}\cdot\nabla\bm{\tilde{\mu}}+\bm{\tilde{\mu}}\cdot\nabla\bm{u}+\bm{\tilde{\mu}}\cdot\nabla\bm{\tilde{\mu}}+\nabla\tilde{\Pi}+\nu\nabla^{2}\bm{\tilde{\mu}}$
For QRM, we introduce a small hyperdiffusion term, approximating eqn 18 as
$\displaystyle 0$ $\displaystyle=\nabla\cdot\bm{\tilde{\mu}}$ (20)
$\displaystyle\bm{0}$
$\displaystyle=\partial_{t}\bm{\tilde{\mu}}+\bm{u}\cdot\nabla\bm{\tilde{\mu}}+\bm{\tilde{\mu}}\cdot\nabla\bm{u}+\bm{\tilde{\mu}}\cdot\nabla\bm{\tilde{\mu}}+\nabla\tilde{\Pi}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad-\nu(\nabla^{2}+\varepsilon\nabla^{4})\bm{\tilde{\mu}}.$

For this problem, we set $\varepsilon=0.001$, because the oscillatory errors
described in Section III trigger a numerical instability when
$\varepsilon\lesssim 0.001$.

Figure 6: Target (top row) and trial vorticity snapshots in the upper half of
the domain at simulation times $t=0,5,10,20$. The second, third, and bottom
row trial solutions are obtained by applying DAL over 200 iterations. GD
$\mathcal{J}^{\bm{u}}_{f}$ (second row) uses gradient descent to minimize
$\mathcal{J}^{\bm{u}}_{f}$. For this case, we use double resolution in space
and time (256 by 512 modes, $\Delta t=0.001$) to ensure that each simulation
is resolved. Despite this extra resolution, we observe symmetry breaking when
$t\geq 5$ due to numerical errors. LB $\mathcal{J}^{\bm{u}}_{f}$ (third row)
uses L-BFGS to minimize $\mathcal{J}^{\bm{u}}_{f}$ with the default
resolution. Although the initial condition is under-resolved in this case, we
still recover a final state which resembles the target. GD
$\mathcal{J}^{\omega}_{f}$ (bottom row) uses gradient descent to minimize
$\mathcal{J}^{\omega}_{f}$. The DAL trial initial conditions (bottom three
rows) have small-scale banded features with high concentrations of positive
vorticity. Their final states have varying degrees of resemblance to the
target at $t=20$.

###  IV.3 Results

We initialize each iterative method using an initial guess
$\bm{u}^{0}(\bm{x},0)=\bm{0}$. For DAL, we minimize the velocity error using
gradient descent (GD $\mathcal{J}^{\bm{u}}_{f}$) and L-BFGS (LB
$\mathcal{J}^{\bm{u}}_{f}$). We run GD $\mathcal{J}^{\bm{u}}_{f}$ at double
resolution to ensure that the trial state’s evolution is well-resolved. We
also minimize the vorticity error via gradient descent (GD
$\mathcal{J}^{\omega}_{f}$). L-BFGS minimization of $\mathcal{J}^{\omega}_{f}$
frequently gives unstable trial initial conditions which cannot be evolved
without additional resolution. Figure 6 consists of vorticity snapshots
belonging to the target solution and trial solutions at $t=0,5,10,20$. The
trial solutions in these cases are obtained using DAL until iteration $n=200$.
GD $\mathcal{J}^{\bm{u}}_{f}$ (second row) contains highly concentrated bands
of positive vorticity encompassing a larger, lower-amplitude vortex. Although
this case is resolved, we observe symmetry breaking when $t\gg 0$. This trial
solution adequately captures the target at $t=20$, but not as well as LB
$\mathcal{J}^{\bm{u}}_{f}$. LB $\mathcal{J}^{\bm{u}}_{f}$ (third row) is
clearly under-resolved at $t=0$, with the same highly-concentrated features as
GD $\mathcal{J}^{\bm{u}}_{f}$ along with small-scale ringing throughout the
domain. These ringing features subside after a brief interval ($t\geq 5$),
after which, the trial solution approaches the target as $t\to 20$. GD
$\mathcal{J}^{\omega}_{f}$ (bottom row) is well-resolved, but its vortex’s
amplitude is appreciably lower than that of the target. In this case, we still
observe band-like regions of concentrated vorticity encompassing a larger low-
amplitude vortex. The DAL trial initial conditions (illustrated in the left
column, lower three rows) do not share similar structures with the target (top
left). Despite this, advection sweeps the concentrated vorticity bands into a
single coherent vortex after 20 advective times. This process illustrates ill-
conditioning of the Navier–Stokes inverse problem.

Figure 7: Same as Figure 6 using the SBI and QRM iterative methods. Despite
the ill-conditioning which hinders DAL in Figure 6, SBI and QRM both recover
the target initial condition (top left) to a greater extent. SBI1 and SBI200
(middle two rows) denote SBI trial states after 1 and 200 iterations
respectively. The QRM initial state (bottom left) loosely resembles the target
initial condition at $n=200$.

Figure 7 illustrates the analogous trial vorticity snapshots using SBI and
QRM. In the top row, we reproduce the target. The second row (SBI1) denotes
the trial state after a single SBI iteration. Its initial state resembles the
target’s better than any of the DAL runs in Figure 6. Its final state is
nearly indistinguishable from the target. The third row (SBI200) denotes the
SBI trial solution after 200 iterations. The bottom row (QRM) denotes the
trial solution after 200 QRM iterations. This case shares the most resemblance
with the target throughout the time domain. At $t=0$, QRM has curved arm-like
features which mimic the target’s initially confined shear layer. At $t\geq
5$, this trial solution resembles the target. In general, SBI and QRM capture
the target state better than DAL without introducing the undesired features
shown in Figure 6.

Figure 8: Velocity error
$\mathcal{J}^{\bm{u}}(t)\equiv\frac{1}{2}\big{\langle}|\bm{u}(\bm{x},t)-\bm{\overline{u}}(\bm{x},t)|^{2}\big{\rangle}$
vs. time of Navier–Stokes trial solutions at iteration $n=200$. GD
$\mathcal{J}^{\bm{u}}_{f}$ and LB $\mathcal{J}^{\bm{u}}_{f}$ denote DAL
minimization of $\mathcal{J}^{\bm{u}}_{f}$ via gradient descent and L-BFGS
respectively. SBI200 and QRM denote the trial solutions refined by SBI and QRM
respectively. Every trial solution’s error descreases nonmonotonically with
respect to time. QRM has the smallest error at the initial and final times,
though SBI200 has less error over large intervals in the time domain.

Figure 8 shows the time-dependent velocity error
$\mathcal{J}^{\bm{u}}(t)\equiv\frac{1}{2}\langle|\bm{u}(\bm{x},t)-\bm{\overline{u}}(\bm{x},t)|^{2}\rangle$
for the trial solutions GD $\mathcal{J}^{\bm{u}}_{f}$, LB
$\mathcal{J}^{\bm{u}}_{f}$, SBI, and QRM at iteration $n=200$. Due to the
inverse-problem’s ill-conditioning, all trial states begin with substantially
larger errors than their respective final states. These errors decrease
nonmonotonically on $t:0\to 20$. The DAL trial solutions’ (GD
$\mathcal{J}^{\bm{u}}_{f}$ and LB $\mathcal{J}^{\bm{u}}_{f}$) errors increase
on $0<t<2$ as their small-scale vortex bands (illustrated in the left column
of Figure 6) rapidly disappear. LB $\mathcal{J}^{\bm{u}}_{f}$ performs
significantly better than GD $\mathcal{J}^{\bm{u}}_{f}$ throughout most of the
time domain, especially near the final state. In contrast, the velocity errors
of SBI and QRM decrease immediately after initialization. Throughout the time
domain, SBI and QRM perform better than either DAL case.

Figure 9: Initial (top, $\mathcal{J}^{\bm{u}}_{0}$) and final (bottom,
$\mathcal{J}^{\bm{u}}_{f}$) velocity errors for the 2D Navier–Stokes inverse
problem over 20 iterations starting with an initial guess
$\bm{u}^{0}(\bm{x},0)=\bm{0}$. We compare how Direct Adjoint Looping (DAL),
Simple Backward Integration (SBI), and the Quasi-Reversible Method (QRM)
minimize these errors. For DAL we compute the gradients of
$\mathcal{J}^{\bm{u}}_{f}$ and $\mathcal{J}^{\omega}_{f}$ (shown black and
grey). Using the gradient of $\mathcal{J}^{\bm{u}}_{f}$, we implement gradient
descent and L-BFGS (solid black and solid grey respectively).
$\mathcal{J}^{\bm{u}}_{f}$ decreases rapidly when SBI and QRM are used. Their
respective performances are comparable in the first 20 iterations.

Figure 9 plots the initial ($\mathcal{J}^{\bm{u}}_{0}$, top) and final
($\mathcal{J}^{\bm{u}}_{f}$, bottom) velocity errors vs. iteration $n$ for
each iterative method until $n=20$. Every method’s performance is comparable
in terms of $\mathcal{J}^{\bm{u}}_{0}$ except GD $\mathcal{J}^{\omega}_{f}$
which does not appreciably decrease this error. In terms of
$\mathcal{J}^{\bm{u}}_{f}$, SBI (light blue solid) and QRM (dark blue dotted)
both outperform every DAL case by several orders of magnitude.

Figure 10: Same as Figure 9 extended over 200 loops. DAL minimization of
$\mathcal{J}^{\omega}_{f}$ (dotted black) is ineffective in terms of both
$\mathcal{J}^{\bm{u}}_{0}$ (top) and $\mathcal{J}^{\bm{u}}_{f}$ (bottom).
After the initial 20 iterations, every other method decreases
$\mathcal{J}^{\bm{u}}_{0}$ by less than an order of magnitude. DAL
minimization of $\mathcal{J}^{\bm{u}}_{f}$ using L-BFGS (grey, LB
$\mathcal{J}^{\bm{u}}_{f}$) is the most effective DAL case in terms of both
errors. At $n=200$, LB $\mathcal{J}^{\bm{u}}_{f}$ does not minimize
$\mathcal{J}^{\bm{u}}_{f}$ to the same extent as SBI and QRM do in very few
iterations. Figures 6 and 7 illustrate their corresponding trial vorticities,
where it is clear that LB $\mathcal{J}^{\bm{u}}_{f}$ is under-resolved and
does not approximate the target initial state.

Figure 10 extends Figure 9 to iteration $n=200$. Both DAL gradient descent
algorithms (GD $\mathcal{J}^{\bm{u}}_{f}$ and GD $\mathcal{J}^{\omega}_{f}$)
stall after $n=20$, whereas LB $\mathcal{J}^{\bm{u}}_{f}$ continues to
minimize the final error $\mathcal{J}^{\bm{u}}_{f}$. The relative success of
LB $\mathcal{J}^{\bm{u}}_{f}$ indicates that the gradient of
$\mathcal{J}^{\bm{u}}_{f}$ has large curvature. SBI and QRM continue to
outperform all DAL runs, with QRM decreasing $\mathcal{J}^{\bm{u}}_{0}$
slightly more. With QRM, the errors decrease nonmonotonically whereas SBI
exhibits a smooth decline.

Whereas previous work used SBI and QRM to perform a single loop, we
demonstrate in Figures 9 and 10 that these methods can be iterated. SBI and
QRM both decrease the errors by more than an order of magnitude in a single
iteration. For $n>150$ in Figure 10, SBI and QRM continue to gradually
minimize $\mathcal{J}^{\bm{u}}_{f}$. This prolonged minimization of
$\mathcal{J}^{\bm{u}}_{f}$ using iterative SBI is a surprising result of our
study. In terms of $\mathcal{J}^{\bm{u}}_{0}$, every method fails to recover
the target state due to ill-conditioning of this inverse problem. QRM
surpasses the other methods in this respect, as illustrated in Figure 10.

##  V Derivation of Ill-Posed Nonlinear Adjoint Equation

Consider the KdVB inverse problem from Sections II and III. Throughout this
investigation, we assert that SBI and QRM perform better than DAL because they
approximate the initial deviation rather than computing a gradient. In fact,
the gradient of $\mathcal{J}^{u}_{0}$ taken with respect to $u(x,0)$ is the
unknown deviation $u^{\prime}(x,0)=u(x,0)-\overline{u}(x,0)$.
$\mathcal{J}^{u}_{0}$ measures the squared distance between the trial and
target initial states. Its levelsets form hyperspheres centered around
$\overline{u}(x,0)$ and its paths of steepest descent are radial (pointing
directly at $\overline{u}(x,0)$). The gradient of $\mathcal{J}^{u}_{0}$ embeds
global information, whereas the gradient of $\mathcal{J}^{u}_{f}$ gives local
information in the space of trial initial conditions.

Using optimal control, we will demonstrate that the ill-posed backward problem
(eqn 9) is equivalent to the adjoint which minimizes $\mathcal{J}^{u}_{0}$.
First we multiply eqn 9 by $u^{\prime}$ and integrate over the space-time
domain, yielding $\mathcal{J}^{u}_{0}$ in terms of $u^{\prime}$, $u$, and
$\mathcal{J}^{u}_{f}$.

$\displaystyle\mathcal{J}^{u}_{0}=\mathcal{J}^{u}_{f}-\int_{0}^{t_{f}}\langle
au^{\prime}\partial_{x}^{2}u^{\prime}-u^{\prime}\partial_{x}[uu^{\prime}]\rangle
dt.$ (21)

Eqn 21 amounts to an energy equation for the deviation between our target and
trial states. Notice how terms which are conservative of the $L_{2}$ norm
cancel ($u^{\prime}\partial_{x}u^{\prime}$ and $b\partial_{x}^{3}u^{\prime}$)
whereas terms with non-zero energy flux ($\partial_{x}[uu^{\prime}]$ and
$a\partial_{x}^{2}u^{\prime}$) persist. Using eqn 21 we construct a new
Lagrangian

$\displaystyle\hat{\mathcal{L}}$
$\displaystyle=\int_{0}^{t_{f}}\langle\hat{\mu}(x,t)\cdot\mathcal{F}[u(x,t)]\rangle
dt+\mathcal{J}^{u}_{0}$ (22)
$\displaystyle=\int_{0}^{t_{f}}\Big{\langle}\hat{\mu}\big{[}\partial_{t}u+u\partial_{x}u-a\partial_{x}^{2}u+b\partial_{x}^{3}u\big{]}$
$\displaystyle\qquad\qquad\qquad\qquad-
au^{\prime}\partial_{x}^{2}u^{\prime}+u^{\prime}\partial_{x}[uu^{\prime}]\Big{\rangle}dt+\mathcal{J}^{u}_{f}.$

with a corresponding adjoint variable $\hat{\mu}$. Taking the variation of
$\hat{\mathcal{L}}$ with respect to $\hat{\mu}(x,t)$ still returns the
original constraint eqn 4 (on $0<t<t_{f}$) and compatibility condition eqn 7
(at $t=t_{f}$). However, varying $\hat{\mathcal{L}}$ with respect to $u(x,t)$
along $0<t<t_{f}$ gives a new adjoint equation

$\displaystyle 0$
$\displaystyle=-\partial_{t}\hat{\mu}-u\partial_{x}\hat{\mu}-a\partial_{x}^{2}\hat{\mu}-b\partial_{x}^{3}\hat{\mu}$
$\displaystyle\qquad-2a\partial_{x}^{2}u^{\prime}+u^{\prime}\partial_{x}u-u^{\prime}\partial_{x}u^{\prime}.$
(23)

Finally, by summing eqns 9 and 23, it is evident that
$\hat{\mu}(x,t)+u^{\prime}(x,t)=0$ not just at $t=t_{f}$, but throughout the
time domain. We substitute $\hat{\mu}\to-u^{\prime}$ in $\hat{\mathcal{L}}$’s
adjoint, yielding

$\displaystyle
0=\partial_{t}\hat{\mu}-a\partial_{x}^{2}\hat{\mu}+b\partial_{x}^{3}+u\partial_{x}\hat{\mu}+\hat{\mu}\partial_{x}u-\hat{\mu}\partial_{x}\hat{\mu}.$
(24)

This new adjoint (equivalent to eqn 9) is ill-posed and effectively nonlinear.
The same demonstration can be applied for any retrospective inverse problem
with an advective constraint, including the Navier–Stokes inverse problem
outlined in Section IV. If eqn 24 could be solved, we would obtain the exact
deviation in our trial initial condition by performing a single adjoint loop.
For example, the non-diffusive Korteweg–de Vries equation and the inviscid
Euler equation admit well-posed nonlinear adjoints as these equations are
time-reversible. However when $a>0$ (diffusive forward problem), eqn 24 cannot
be solved backwards using conventional numerical solvers. The ill-posed
backward integration problem motivates the SBI and QRM modifications. As
$\varepsilon\to 0$, QRM backward integration approaches the ill-posed
nonlinear adjoint which minimizes $\mathcal{J}^{u}_{0}$.  From this
perspective, QRM is a regularized form of DAL. SBI backward integration
combines the linear adjoint which minimizes $\mathcal{J}^{u}_{f}$ with the
additional terms appearing in the ill-posed nonlinear adjoint.

Given that the SBI and QRM backward integrations return approximations for the
gradient of $\mathcal{J}^{u}_{0}$, it follows that we can repeatedly apply the
approximate gradients to refine our trial state. Taking the dynamical systems
perspective, we can treat these iterative methods as discrete maps where the
target state is a fixed point. Our results demonstrate that these fixed points
are generally attracting for advection–diffusion systems. Potential
difficulties might arize in future investigations due to limit-cycles or even
strange attractors. We have simultaneously highlighted an important limitation
of using gradients to optimize initial conditions: gradients can have large
curvatures, even when the trial state is in close proximity to the extremum
(see Appendices A and B). By including additional terms in our backward
integrations, we relax the gradient’s sensitive dependence on the trial
initial condition by replacing the gradient with a smoother discrete map.
Higher-order algorithms for finding fixed points could also be paired with SBI
and QRM to further accelerate convergence.

##  VI Discussion and Conclusions

###  VI.1 Summary of Results

In this investigation, we extend several previous studies involving Direct
Adjoint Looping (DAL), Simple Backward Integration (SBI), and the Quasi-
Reversible Method (QRM) by applying them to a pair of retrospective inverse
problems. DAL has been used previously to solve inverse problems and optimize
simulation inputs, as it allows us to compute the gradient of an arbitrary
functional. This gradient can be used to refine the trial initial condition
toward a local extremum. However, the number of iterations required to
approximate said extremum depends on the user’s initial guess.

The SBI and QRM methods are designed to approximate the trial solution’s
deviation rather than taking a gradient. Using the approximate deviation at
the initial state, we update the trial initial condition using a “step size”
of unity. In practice, we cannot compute the initial deviation because the
necessary backward integration is ill-posed. SBI and QRM provide well-posed
approximations for the ill-posed backward integration.

SBI has been used to make coarse approximations and generate initial guesses
for further refinement via DAL. We extend this use by developing an iterative
algorithm which acts on the trial state’s deviation. With this method, we
evolve the trial state’s final deviation $u^{\prime}(x,t_{f})$ backwards in
time using a modified constraint equation. This modified constraint equation
is well-posed because we reverse the sign of the diffusive term.

Iterative QRM has previously been implented. With this backward integration,
we preserve the ill-posed diffusion term by introducing a small hyperdiffusion
term to relax the numerical instability. We must also tune this hyperdiffusion
term’s coefficient $\varepsilon$ according to the problem’s dynamics and
numerical resolution.

Our first inverse problem is constrained by the Korteweg–de Vries–Burgers
(KdVB) equation. The target solution consists of a nonlinear wave whose
advective velocity decays due to diffusion. We demonstrate that the gradient
of $\mathcal{J}^{u}_{f}$, obtained via DAL or otherwise, is not an efficient
means of recovering the target initial state. DAL performs poorly when coupled
with gradient descent as well as the L-BFGS optimization routine. When
following the direction of steepest descent in the space of trial initial
conditions, we accumulate undesired wave features. This is due to high
curvature in the path of steepest descent. QRM captures the target final state
almost immediately while SBI does so more gradually.

The second inverse problem is constrained by the 2D incompressible
Navier–Stokes equation. The target solution is a symmetric Kelvin–Helmholtz
vortex with $\rm{Re}\equiv\nu^{-1}=50,000$. We evolve the target state over 20
advective time units, then use the final state to recover its associated
initial condition. The problem domain is doubly-periodic which likely reduces
the ill-conditioning of our retrospective inverse problem. We apply each
iterative method until iteration $n=200$. Here we utilize the flexibility of
DAL by minimizing the velocity error $\mathcal{J}^{\bm{u}}_{f}$ as well the
vorticity error $\mathcal{J}^{\omega}_{f}$. DAL minimization of
$\mathcal{J}^{\omega}_{f}$ performs poorly as illustrated in Figures 6, 9 and
10. Using L-BFGS to minimize $\mathcal{J}^{\bm{u}}_{f}$, we produce an initial
condition whose final state resembles the target after 200 iterations (Figure
6). This L-BFGS case successfully minimizes $\mathcal{J}^{\bm{u}}_{f}$, as
illustrated in Figure 10. However, the trial initial condition is under-
resolved because L-BFGS accumulates small-scale features. This highlights
severe ill-conditioning in the Navier–Stokes inverse problem, as we are able
to minimize the final error without approximating the target initial
condition. We implement gradient descent with double resolution to demonstrate
this more clearly. This gradient descent case does not approximate the target
final state as well as L-BFGS, even though the simulations are well-resolved.

Just as with KdVB, SBI and QRM outperform DAL by every metric in the
Navier–Stokes inverse problem. Although every method decreases
$\mathcal{J}^{\bm{u}}_{0}$ with a similar trend (Figures 9 and 10), SBI and
QRM (Figure 7) approximate the target initial state better than DAL (Figure 6)
at $n=200$. The methods’ respective behaviors differ in close proximity to
their target, as illustrated in Figure 10. Specifically, QRM exhibits
nonmonotonicity in both $\mathcal{J}^{\bm{u}}_{0}$ and
$\mathcal{J}^{\bm{u}}_{f}$ whereas SBI decreases these quantities with each
iteration.

###  VI.2 Future Work

Optimal control via DAL is the most common approach for retrospective inverse
problems. However, SBI and QRM are more effective than DAL when applied to the
2D Navier–Stokes inverse problem.  In future work, we will implement SBI and
QRM for 3D Navier–Stokes inverse problems. Chaos and the direct energy cascade
(as opposed to the inverse energy cascade in 2D) will likely present
additional challenges. We will also use SBI and QRM to solve other optimal
control problems, such as the mantle convection inverse problems studied by
[4] and [5]. Numerical weather prediction (NWP) could also be enhanced by
implementing these nonlinear methods. SBI and QRM are well-suited for this
application because NWP models are increasingly nonlinear. Predictions are
limited by our ability to incorporate real-time observations in a reasonable
number of iterations [30].

Another relevant study is [31], which computes time-periodic solutions of
vortex sheets with surface tension. Using DAL, they minimize a cost functional
which resembles $\mathcal{J}^{u}_{f}$. Although they are not afforded a target
final state, we can interpret their adjoint loop as a dynamic inversion where
the target state changes with each iteration. Conceptually, any properly-
constrained optimal control problem has a corresponding retrospective inverse
problem. The compatibility condition may differ (e.g. [31] initializes their
adjoint with the deviation between their trial initial and final states), but
in general the adjoint evolution equations remain the same. Given that our SBI
and QRM loops tend toward the mutual extremum of $\mathcal{J}^{u}_{0}$ and
$\mathcal{J}^{u}_{f}$, it is reasonable to suggest that these methods might
also have attracting extrema when seeking periodic solutions. Efficient
algorithms for computing time-periodic solutions would have direct
applications in studying cyclic dynamos and anticipating nonlinear resonances
in engineering systems such as turbines.

## Acknowledgments

The authors thank Emma Kaufman, Benjamin Hyatt, Whitney Powers, Ilaria
Fontana, and Adrian Fraser for their valuable suggestions. We thank the
Dedalus development team. L.O., D.L., K.A., K.J.B., J.S.O., and B.P.B. are
supported in part by NASA HTMS grant 80NSSC20K1280. D.L., K.J.B., J.S.O., and
B.P.B. are supported in part by NASA OSTFL grant 80NSSC22K1738. L.O. is
supported by the National Science Foundation Graduate Research Fellowship
under Grant No. (DGE-2234667). Computations were conducted with support by the
NASA High End Computing Program through the NASA Advanced Supercomputing (NAS)
Division at Ames Research Center on Pleiades with allocation GID s2276.

##  Appendix A KdVB Inverse Problem Objectives using SBI Initial Guess

Figure 11: Same is Figure 3 using the SBI initial guess (shown in green).

For the KdVB inverse problem described in Section III, we initialize
$u^{0}(x,0)$ using an initial guess constructed via SBI. This SBI initial
guess is computed by evolving $\overline{u}(x,t_{f})$ backwards using a
modified constraint equation (eqn 4 with the diffusive term’s sign reversed).
For DAL, this implementation mimics [5, 4] who demonstrated that their SBI
initial guess converges more rapidly compared to a neutral guess of zero.

Using the SBI initial guess does improve our DAL approximations of
$\overline{u}(x,0)$, as shown in the top row of Figure 11. L-BFGS (top right)
continues to outperform gradient descent (top left), though both cases have
significant deviations from $\overline{u}(x,0)$. SBI and QRM both outperform
the DAL cases, converging in approximately the same number of iterations as in
Section III.

Figure 12 plots the initial (top) and final (bottom) errors for the KdVB
inverse problem using the SBI initial guess. LB $\mathcal{J}^{u}_{f}$ (shown
in grey) performs appreciably better than in Figure 5. However SBI and QRM
both surpass DAL, with QRM performing best by several orders of magnitude in
terms of $\mathcal{J}^{u}_{f}$.

Figure 12: Same is Figure 5 using initial guess $u^{0}(x,0)$ generated via
SBI.

##  Appendix B Navier–Stokes Inverse Problem Objectives using SBI Initial
Guess

For the Navier–Stokes inverse problem described in Section IV, we compare each
method (DAL, SBI, and QRM) starting with an initial guess generated via SBI.
This case mimics [5, 4], who began by rewinding their target state from
$t:t_{f}\to 0$ using SBI. They subsequently refined their SBI initial guess by
using DAL to minimize an error functional (analogous to DAL
$\mathcal{J}^{\bm{u}}_{f}$, shown in black). Figures 13 and 14 plot the
initial ($\mathcal{J}^{\bm{u}}_{0}$, top) and final
($\mathcal{J}^{\bm{u}}_{f}$, bottom) velocity errors as functions of iteration
$n$. Figure 13 shows the first 20 iterations while Figure 14 terminates at
$n=200$. In the black curve (GD $\mathcal{J}^{\bm{u}}_{f}$), we use DAL with
gradient descent to minimize $\mathcal{J}^{\bm{u}}_{f}$.

Figure 13: Same as Figure 9 using initial guess $\bm{u}^{0}(\bm{x},0)$
generated via SBI. Figure 14: Same as Figure 10 using initial guess
$\bm{u}^{0}(\bm{x},0)$ generated via SBI.

First, we confirm that using an SBI initial guess with DAL does approach the
target final state faster than in Section IV. Next we implement SBI and QRM as
before. Just as with the previous initial guess $\bm{u}^{0}(\bm{x},0)=\bm{0}$,
we find that SBI and QRM outperform DAL, albeit by a narrower margin (ranging
from less than one to two orders of magnitude in terms of
$\mathcal{J}^{\bm{u}}_{f}$). Crucially, SBI and QRM both outperform DAL in
very few iterations, suggesting that the efficacy of these methods is somewhat
robust to the initial guess. Every iterative method increases the initial
error $\mathcal{J}^{\bm{u}}_{0}$ (top) compared to the SBI initial guess,
except for QRM which eventually decreases $\mathcal{J}^{\bm{u}}_{0}$ below its
original value near $n\approx 190$ (shown in dotted dark blue). Iterative SBI
appears to follow a similar trend at a slower pace, whereas GD
$\mathcal{J}^{\bm{u}}_{f}$ shows no such trend. This behavior with GD
$\mathcal{J}^{\bm{u}}_{f}$ indicates that our Navier–Stokes inverse problem is
highly non-convex, even in close proximity to the extremum. In this case,
L-BFGS minimization of $\mathcal{J}^{\bm{u}}_{f}$ gives unstable trial initial
conditions. We obtain similar results by initializing each iterative method
using the QRM trial initial state at $n=200$.

## References

## References

  * Lukyanenko  _et al._ [2021] D. Lukyanenko, T. Yeleskina, I. Prigorniy, T. Isaev, A. Borzunov, and M. Shishlenin, Inverse problem of recovering the initial condition for a nonlinear equation of the reaction–diffusion–advection type by data given on the position of a reaction front with a time delay, Mathematics 9, 1 (2021).
  * Subramanian  _et al._ [2020] S. Subramanian, K. Scheufele, M. Mehl, and G. Biros, Where did the tumor start? an inverse solver with sparse localization for tumor growth models, Inverse Problems 36, 10.1088/1361-6420/ab649c (2020).
  * Kabanikhin  _et al._ [2020] S. I. Kabanikhin, I. M. Kulikov, and M. A. Shishlenin, An algorithm for recovering the characteristics of the initial state of supernova, Computational Mathematics and Mathematical Physics 60, 10.1134/S0965542520060135 (2020).
  * Liu and Gurnis [2008] L. Liu and M. Gurnis, Simultaneous inversion of mantle properties and initial conditions using an adjoint of mantle convection, Journal of Geophysical Research: Solid Earth 113, 10.1029/2008jb005594 (2008).
  * Li  _et al._ [2017] D. Li, M. Gurnis, and G. Stadler, Towards adjoint-based inversion of time-dependent mantle convection with nonlinear viscosity, Geophysical Journal International 209, 10.1093/gji/ggw493 (2017).
  * Imanuvilov [1997] O. Y. Imanuvilov, Local exact controllability for the 2-d navier-stokes equations with the navier slip boundary conditions, Turbulence Modeling and Vortex Dynamics , 148 (1997).
  * Fan and Nakamura [2009] J. Fan and G. Nakamura, Well-posedness of an inverse problem of navier-stokes equations with the final overdetermination, Journal of Inverse and Ill-Posed Problems 17, 565 (2009).
  * Fan  _et al._ [2010] J. Fan, M. D. Cristo, Y. Jiang, and G. Nakamura, Inverse viscosity problem for the navier-stokes equation, Journal of Mathematical Analysis and Applications 365, 750 (2010).
  * Du  _et al._ [2013] J. Du, I. M. Navon, J. Zhu, F. Fang, and A. K. Alekseev, Reduced order modeling based on pod of a parabolized navier–stokes equations model ii: Trust region pod 4d var data assimilation, Computers & Mathematics with Applications 65, 380 (2013).
  * Law and Stuart [2012] K. J. Law and A. M. Stuart, Evaluating data assimilation algorithms, Monthly Weather Review 140, 3757 (2012).
  * Fan  _et al._ [2020] T. Fan, K. Xu, J. Pathak, and E. Darve, Solving inverse problems in steady-state navier-stokes equations using deep neural networks (2020), arXiv:2008.13074 [math.NA] .
  * Frerix  _et al._ [2021] T. Frerix, D. Kochkov, J. A. Smith, D. Cremers, M. P. Brenner, and S. Hoyer, Variational data assimilation with a learned inverse observation operator (2021), arXiv:2102.11192 [cs.LG] .
  * Zhao  _et al._ [2022] Q. Zhao, D. B. Lindell, and G. Wetzstein, Learning to solve pde-constrained inverse problems with graph networks (2022), arXiv:2206.00711 [cs.LG] .
  * Mowlavi and Nabi [2023] S. Mowlavi and S. Nabi, Optimal control of PDEs using physics-informed neural networks, Journal of Computational Physics 473, 111731 (2023).
  * Mannix  _et al._ [2022] P. M. Mannix, C. S. Skene, D. Auroux, and F. Marcotte, Discrete adjoint-based control: A robust gradient descent procedure for optimisation with pde and norm constraints (2022), arXiv:2210.17194 [physics.flu-dyn] .
  * [16] R. Lattès, J. Lions, and R. Bellman,  _The Method of Quasi-reversibility: Applications to Partial Differential Equations_ (American Elsevier Publishing Company).
  * Conrad and Gurnis [2003] C. P. Conrad and M. Gurnis, Seismic tomography, surface uplift, and the breakup of gondwanaland: Integrating mantle convection backwards in time, Geochemistry, Geophysics, Geosystems 4, 10.1029/2001GC000299 (2003).
  * Fang  _et al._ [2012] L. Fang, W. J. T. Bos, S. Liang, J.-P. Bertoglio, and L. Shao, Time reversibility of navier-stokes turbulence and its implication for subgrid scale models, Journal of Turbulence 13, 1 (2012).
  * Dardé [2016] J. Dardé, Iterated quasi-reversibility method applied to elliptic and parabolic data completion problems, Inverse Problems and Imaging 10, 379 (2016).
  * Nguyen and Klibanov [2022] L. H. Nguyen and M. V. Klibanov, Carleman estimates and the contraction principle for an inverse source problem for nonlinear hyperbolic equations, Inverse Problems 38, 035009 (2022).
  * Imanuvilov  _et al._ [2022] O. Y. Imanuvilov, L. Lorenzi, and M. Yamamoto, Carleman estimate for the Navier–Stokes equations and applications, Inverse Problems 38, 085006 (2022), publisher: IOP Publishing.
  * Steiner and Reichl [2012] W. Steiner and S. Reichl, The optimal control approach to dynamical inverse problems, Journal of Dynamic Systems, Measurement and Control, Transactions of the ASME 134, 10.1115/1.4005365/455882 (2012).
  * Barzilai and Borweigh [1988] J. Barzilai and J. M. Borweigh, Two-Point Step Size Gradient Methods, IMA Journal of Numerical Analysis 8, 141 (1988), https://academic.oup.com/imajna/article-pdf/8/1/141/2402762/8-1-141.pdf .
  * Byrd  _et al._ [1995] R. H. Byrd, P. Lu, J. Nocedal, and C. Zhu, A limited memory algorithm for bound constrained optimization, SIAM Journal on Scientific Computing 16, 1190 (1995), https://doi.org/10.1137/0916069 .
  * Note [2] In addition to the results shown in section III, we perform numerical experiments on the KdVB inverse problem using conjugate-gradient and line search routines (provided by scipy version 1.12.0). In both cases we use the default parameters. For the line search, we use the negative gradient for the descent direction.
  * Note [3] SBI and QRM do not generally return any sensible gradient, so the specified updates are not technically gradient descent. Attempting a constant step size of $\gamma\sim\mathcal{O}(1)$ with DAL is unstable for the direct minimization of $\mathcal{J}^{u}_{f}$. SBI can often be accelerated using a larger step size.
  * Burns  _et al._ [2020] K. J. Burns, G. M. Vasil, J. S. Oishi, D. Lecoanet, and B. P. Brown, Dedalus: A flexible framework for numerical simulations with spectral methods, Physical Review Research 2, 023068 (2020).
  * Nabi  _et al._ [2022] S. Nabi, P. Grover, and C. P. Caulfield, Robust preconditioned one-shot methods and direct-adjoint-looping for optimizing reynolds-averaged turbulent flows, Computers and Fluids 238, 10.1016/j.compfluid.2022.105390 (2022).
  * Kerswell  _et al._ [2014] R. R. Kerswell, C. C. T. Pringle, and A. P. Willis, An optimization approach for analysing nonlinear stability with transition to turbulence in fluids as an exemplar, Reports on Progress in Physics 77, 085901 (2014).
  * Bonavita  _et al._ [2018] M. Bonavita, P. Lean, and E. Holm, Nonlinear effects in 4d-var, Nonlinear Processes in Geophysics 25, 713 (2018).
  * Ambrose  _et al._ [2015] D. Ambrose, M. Kondrla, and M. Valle, Computing time-periodic solutions of a model for the vortex sheet with surface tension, Quarterly of Applied Mathematics 73, 10.1090/s0033-569x-2015-01364-8 (2015).