# A generalized framework to transport geophysical fields: a differential
geometry perspective

Yicun Zhen Corresponding author<EMAIL_ADDRESS>Department of
Oceanography, Hohai University, Nanjing, Jiangsu, China  Valentin Resseguier
LAB, SCALIAN DS, Rennes, France  INRAE, OPAALE, Rennes, France  Bertrand
Chapron  Laboratoire d’Océanographie Physique et Spatiale, Ifremer, Plouzané,
France

###### Abstract

To estimate displacements of physical fields, a general framework is proposed.
Considering that for each state variable, a tensor field can be associated,
ways these displacements act on different state variables will differ
according to the tensor field definitions. This perspective provides a
differential-geometry-based reformulation of the generalized optical flow (OF)
algorithm. Using the proposed framework, optimisation procedures can
explicitly ensure the conservation of certain physical quantities (total mass,
total vorticity, total kinetic energy, etc.). Existence and uniqueness of the
solutions to the local optimisation problem are demonstrated, leading to a new
nudging strategy using all-available observations to infer displacements of
both observed and unobserved state variables. Using the proposed nudging
method before EnKF, numerical results show that ensemble data assimilation
better preserves the intrinsic structure of underline physical processes if
the ensemble members are aligned with the observations.

##  1 Introduction

Since chaotic divergence is an intrinsic property of turbulent geophysical
systems, data assimilation is often necessary to improve model trajectories of
state variables. Numerous strategies have then been proposed, including
nudging methods, 3-Dimensional variational method [17], Kalman filter based
methods [13][14][19][15], 4-Dimensional variational methods [6], and particle
filters [7].

If $x^{b}(t)$ denotes the state vector at time $t$ estimated by the model, and
$y^{o}(t)$ the observed state variables at the same time $t$, a first step is
to quantify differences between $x^{b}(t)$ and $y^{o}(t)$. Different
interpretations of this space of errors result in different data assimilation
algorithms. For illustration, consider $x^{b}$ and $y^{o}$, Fig. 1. The
difference simply relates to shifted positions of two bright spots.
Numerically $y^{o}-x^{b}$ can possibly result in two different nudging
strategies :

Figure 1: Suppose that the difference between the model estimate $x^{b}$ and
observed state $y^{o}$ only differ by the position, then apparently a nudging
strategy that gradually moves $x^{b}$ rightward is more reasonable than the
linear nudging strategy.

  * Strategy 1 (nonlinear): $x^{b}(x)\longleftarrow x^{b}({x}-\epsilon v(x))$;

  * Strategy 2 (linear): $x^{b}(x)\longleftarrow(1-\epsilon)x^{b}(x)+\epsilon y^{o}(x)$,

where $v(x)$ is a vector field determined by $y^{o}$ and $x^{b}$, encoding the
displacement of each point in the domain. The linear strategy is the most
commonly used nudging strategy. However, a nonlinear nudging strategy seems
more natural. It can preserve the main physical feature of the physical field,
e.g. bright spots for this example, during the whole nudging process. In
general, a position error can be represented by an invertible smooth map
$T:\Omega\longrightarrow\Omega$ to represent the displacement of each point in
the domain $\Omega$.

A numerical evaluation of $T$ for two physical fields $S_{1}$ and $S_{2}$, can
then target the minimization of a cost function of the following form:

$\displaystyle\begin{cases}T=\displaystyle\operatorname*{arg\,min}_{T\in
C^{\infty}_{\text{inv}}(\Omega,\Omega)}\|S_{1}-T^{\\#}S_{2}\|^{2}+a\|T\|^{2}\\\
\text{Proper boundary conditions}\end{cases}$ (1)

in which $C^{\infty}_{\text{inv}}(\Omega,\Omega)$ refers to the set of smooth
1-1 and onto maps from $\Omega$ to itself. $T^{\\#}S_{2}$ represents how
$S_{2}$ transforms under the map $T$. $\|T\|$, a chosen norm for the map $T$,
serves as a regularization term. $a>0$ is a pre-chosen positive constant that
weights this regularization term. Directly solving the optimisation problem
(1) may technically be too difficult. Instead, it might be easier to
iteratively solve an "infinitesimal" version :

$\displaystyle\begin{cases}v_{i}=\displaystyle\operatorname*{arg\,min}_{{\bf
v}:\Omega\rightarrow\mathbb{R}^{n}}\frac{\partial\|S_{1}-(Id+sv)^{\\#}S_{2}\|^{2}}{\partial
s}\big{|}_{s=0}+a\|v\|^{2}\\\ \text{ Proper boundary condition for
$v_{i}$},\end{cases}$ (2)

where $v_{i}$ is the vector field at the $i-$th iterative step, $Id$ the
identity map on $\Omega$. For each $i$, $v_{i}$ is the "optimal" vector field
along which $S_{2}$ the most fast transforms towards $S_{1}$. At each
iterative step, $S_{2}$ is replaced with $(Id+\epsilon v_{i})^{\\#}S_{2}$ for
some pre-chosen small $\epsilon$. Iteratively solving (2) results in a
displacement flow $\Phi$:

$\displaystyle\frac{\partial\Phi(s,x)}{\partial
s}\Big{|}_{s=k\epsilon}=v_{k}(\Phi(k\epsilon,x)).$ (3)

And the displacement map $T(x)=\Phi(N\epsilon,x)\approx(Id+\epsilon
v_{N})\circ\dots\circ(Id+\epsilon v_{1})$ is our candidate solution for (1).
However, not every smooth 1-1 and onto map $T$ can be approximated by a
sequential composition of small displacements. Therefore the optimisation
problems (1) and (2) can be fundamentally different. Our present motivation is
on a reformulation and solution of (2).

A key is to design a method to include the definition of $T^{\\#}$ and the
choice of norms for the two terms in Eq.(2). The definition of $T^{\\#}$ is
essential and some relevant algorithms are reviewed or briefly discussed in
section 2 and section 3, including the optical flow (OF) algorithms [1, 29,
16], the large deformation diffeomorphisms [26, 23, 3, 4], the metamorphoses
strategy [25, 24], and some other algorithms [20, 2]. Although some of these
algorithms have been widely applied to geophysical observations, $T^{\\#}$ are
commonly used without considering the dynamics of the geophysical fields. The
displacement map $T$ may then cause severe structural errors when applied to
unobserved state variables. Fig.2 illustrates an example of directly applying
$T$ by composition. The original physical field $S$ is a rotational vector
field, rotating counter-clockwise. Suppose that a displacement map $T$ is to
be estimated from observations of other state variables. $T$ is assumed to be
a clockwise rotation by $90^{\circ}$. A direct application of $T$ then
transforms $S$ to a displacement field $S\circ T^{-1}$ displaying completely
different features.

Figure 2: Suppose that the original field $S$ is a rotating wind field and the
given displacement map $T$ is clockwise rotation by $90^{\circ}$. Then the
direct composition of $S$ and $T$ results in a wind field of completely
different feature.

To circumvent such undesirable results, [30] considered a differential
geometry framework and the use of tensor fields to formulate $T^{\\#}$ in
Eq.(1) and (2). In this new perspective, the choice of tensor fields can then
be explicitly dictated by the dynamical equations of the underline physical
fields. $T^{\\#}$ will follow dynamical principles, and certain physical
quantities are naturally conserved during the morphing process. A similar
argument has been given for small displacement cases [28]. This further leads
to a new alignment strategy in which the displacement flow calculated from the
observed state variables can also be applied to partially correct the
displacement of the physical fields of unobserved state variables. Such a new
alignment strategy can naturally be incorporated with the classical ensemble
Kalman filter to reduce inherent difficulties arising from linear algorithms.

The constrained formulation of $T^{\\#}$, based on the concept of tensor
fields is given in section 2. This leads to a new version of the optimisation
problem in the form of (2). Existence and uniqueness of the resulting solution
are provided. A physical interpretation of the new optimisation problem is
discussed, with a brief review of some classical OF algorithms [1, 29, 16].
These OF algorithms are then compared with the proposed algorithm. The large
deformation diffeomorphism strategy [3, 4, 26, 23] and the metamorphoses
strategy [24, 25] are also discussed in section 2. The new nudging strategy
and its associated data assimilation strategy are presented in section 3.
Differences between the proposed data assimilation strategy and some existing
methods [20, 2] are also discussed. Using the thermal shallow water dynamical
framework, numerical results are serving to support the proposed developments.
Conclusion is given in Section 4. The complete code to reproduce the numerical
results in this paper is available at 10.5281/zenodo.10252176.

##  2 A differential geometry formulation of the optimisation problem

Let $(\Omega,g)$ be a compact oriented Riemannian manifold of dimension $n$
with or without boundary, in which $g$ is the Riemannian metric. For any
smooth vector field $u$, $\theta_{u}$ is the differential 1-form such that
$\theta_{u}(v)=g(u,v)=\langle u,v\rangle$ for any smooth vector field $v$. Let
$T\Omega$ be the tangent bundle of $\Omega$ and $T^{*}\Omega$ the cotangent
bundle. For any smooth map $\phi:\Omega\longrightarrow\Omega$, $\phi_{*}$
refers to the push-forward map of $T\Omega$ induced by $\phi$, and $\phi^{*}$
to the pull-back map of $T^{*}\Omega$ induced by $\phi$. We further assume
that $\phi$ is invertible, and

$\displaystyle\phi^{\\#}:T\Omega\cup T^{*}\Omega$
$\displaystyle\longrightarrow T\Omega\cup T^{*}\Omega$ (4) $\displaystyle u$
$\displaystyle\longrightarrow\phi_{*}u$ (5) $\displaystyle\omega$
$\displaystyle\longrightarrow(\phi^{-1})^{*}\omega,$ (6)

for any $u\in T\Omega$ and $\omega\in T^{*}\Omega$. It further induces an
isomorphism of tensor fields of any specific type:

$\displaystyle\phi^{\\#}:V_{1}\otimes V_{2}\otimes\cdots\otimes V_{l}$
$\displaystyle\longrightarrow V_{1}\otimes V_{2}\otimes\cdots\otimes V_{l}$
(7) $\displaystyle\alpha_{1}\otimes\alpha_{2}\otimes\cdots\otimes\alpha_{l}$
$\displaystyle\longrightarrow(\phi^{\\#}\alpha_{1})\otimes\cdots\otimes(\phi^{\\#}\alpha_{l})$
(8)

in which $V_{i}$ can either be a copy of $T\Omega$ or a copy of $T^{*}\Omega$.

For the smooth vector field $u$ such that $u\big{|}_{\partial\Omega}\in
i_{*}(T\partial\Omega)$, where $i:\partial\Omega\longrightarrow\Omega$ is the
natural embedding map of the boundary,
$\Phi_{u}:[0,\epsilon]\times\Omega\longrightarrow\Omega$ denotes the flow that
satisfies

$\displaystyle\frac{\partial\Phi_{u}(s,x)}{\partial s}=u(\Phi_{u}(s,x)).$ (9)

The Lie derivative of any tensor field $\theta$ with respect to $u$ is:

$\displaystyle\mathcal{L}_{u}\theta=\lim_{s\to
0}\frac{\Phi_{u}(s)^{\\#}\theta-\theta}{s},$ (10)

which is also a tensor field, of the same type as $\theta$. To first order, a
Taylor expansion gives:

$\displaystyle\Phi_{u}(s)^{\\#}\theta=\theta+s\mathcal{L}_{u}\theta+o(s).$
(11)

The Riemannian metric $g$ provides the inner product for the tangent space
$T_{x}\Omega$ for any $x\in\Omega$. This inner product could be generalized to
the space of tensor fields of any specific type. For $\theta_{1},\theta_{2}\in
V_{1}\otimes\cdots\otimes V_{k}$, the inner product induced by $g$ at
$x\in\Omega$ is written as $\langle\theta_{1},\theta_{2}\rangle_{x}$. The
corresponding norm writes
$|\theta_{2}|_{x}^{2}=\langle\theta_{2},\theta_{2}\rangle_{x}$. The Riemannian
metric further induces several operators. Let $*$ be the Hodge star operator,
$d$ the exterior derivative, and $\delta=(-1)^{n(k-1)+1}*d*$ the
codifferential. These operators act on differential forms. The Laplacian
operator is defined as $\Delta=d\delta+\delta d$. Let $dV$ be the volume form
on $\Omega$.

Recall that $\theta_{u}$ is a differential 1-form determined by vector field
$u$ and the Riemannian metric tensor $g$. With the above notations, we can
generalize the OF method for tensor fields on a Riemannian manifold:

######  Definition 2.0.1.

Let $\theta_{1}$ and $\theta_{2}$ be two tensor fields on $\Omega$ of the same
type. When $\Omega$ is a compact oriented Riemannian manifold without
boundary, solve for

$\displaystyle u=\operatorname*{arg\,min}_{u\in
T\Omega}\int_{\Omega}\bigg{[}\frac{\partial|\theta_{1}-\Phi_{u}(s)^{\\#}\theta_{2}|_{x}^{2}}{\partial
s}\Big{|}_{s=0}+a_{1}|d\theta_{u}|_{x}^{2}+a_{1}|\delta\theta_{u}|_{x}^{2}+a_{0}|\theta_{u}|_{x}^{2}\bigg{]}dV.$
(12)

When $\Omega$ is a compact oriented Riemannian manifold with boundary, solve
for

$\displaystyle\begin{cases}u=\displaystyle\operatorname*{arg\,min}_{u\in
T\Omega}\int_{\Omega}\bigg{[}\frac{\partial|\theta_{1}-\Phi_{u}(s)^{\\#}\theta_{2}|_{x}^{2}}{\partial
s}\Big{|}_{s=0}+a_{1}|d\theta_{u}|_{x}^{2}+a_{1}|\delta\theta_{u}|_{x}^{2}\bigg{]}dV\\\
u\big{|}_{\partial\Omega}=i_{*}v,\end{cases}$ (13)

where $v\in T\partial\Omega$ is the given boundary condition and $i$ is the
natural embedding of $\partial\Omega$.

Using the Taylor expansion of $\Phi_{u}^{\\#}\theta_{2}$, see Eq.(11), Eq.(12)
is equivalent to

$\displaystyle u=\operatorname*{arg\,min}_{u\in
T\Omega}\int_{\Omega}\bigg{[}-2\langle\theta_{1}-\theta_{2},\mathcal{L}_{u}\theta_{2}\rangle_{x}+a_{1}|d\theta_{u}|_{x}^{2}+a_{1}|\delta\theta_{u}|_{x}^{2}+a_{0}|\theta_{u}|_{x}^{2}\bigg{]}dV,$
(14)

and Eq.(13) equivalent to

$\displaystyle\begin{cases}u=\displaystyle\operatorname*{arg\,min}_{u\in
T\Omega}\int_{\Omega}\bigg{[}-2\langle\theta_{1}-\theta_{2},\mathcal{L}_{u}\theta_{2}\rangle_{x}+a_{1}|d\theta_{u}|_{x}^{2}+a_{1}|\delta\theta_{u}|_{x}^{2}\bigg{]}dV\\\
u\big{|}_{\partial\Omega}=i_{*}v.\end{cases}$ (15)

######  Theorem 2.1.

For $\theta_{1},\theta_{2},$ and $v$ with finite $H^{1}$ norm, i.e.
$\int_{\Omega}|\theta_{i}|_{x}^{2}+|d\theta_{i}|_{x}^{2}+|\delta\theta_{i}|_{x}^{2}<\infty$
for $i=1,2$, and similar for $v$, the optimisation problems (14) and (15) are
always uniquely solvable in the space $H^{1}(T\Omega)$. And the solution is
twice differentiable.

Theorem 2.1 is a direct consequence of some proven mathematical results in [9]
(or [22]). A complete demonstration is provided in the appendix. The
formulation of (12) and (13), together with theorem (2.1) provide an option
for the theory of defining and computing the displacement flow of two tensor
fields. Due to Poincare lemma, both the regularization terms in (12) and (13)
are equivalent to the $H^{1}$ norm. Still, the regularization term does not
have to be the $H^{1}$ norm. For instance,
$|d\delta\theta_{u}|_{x}^{2}+|\delta d\theta_{u}|_{x}^{2}$ is considered in
[5]. However, the vector field $u$ calculated from Eq.(12) or (13) should not
necessarily be the same as the true velocity field, largely depending on the
regularization terms. True physical laws may not always be well represented by
such regularization terms.

###  2.1 A detailed formulation of $T^{\\#}$ for physical fields

We can now write down a more explicit expression for (2). Given two physical
fields $S_{1}(x)$ and $S_{2}(x)$ that represent the same quantity $S$, we
first associate a tensor to this state variable $S$. Thus $S_{1}$ and $S_{2}$
corresponds to two tensor fields, $\theta_{1}$ and $\theta_{2}$, of the same
type. The displacement vector field $u$ is then calculated based on (12) or
(13) for $\theta_{1}$ and $\theta_{2}$. It is required that the correspondence
between $S$ and the tensor is 1-1 and onto, i.e. an unique $S$ can be inferred
from a given value of the tensor. Thus, for a given displacement map
$T:\Omega\longrightarrow\Omega$, $T^{\\#}S_{2}$ is defined in three steps:

  * 1), construct tensor field $\theta_{2}$ based on $S_{2}$;

  * 2), $\theta_{2}^{\text{new}}\leftarrow T^{\\#}\theta_{2}$;

  * 3), $T^{\\#}S_{2}$ is then inferred from $\theta_{2}^{\text{new}}$.

Under $T^{\\#}$, certain physical quantities will naturally be conserved as
long as $T$ is an diffeomorphism of $\Omega$ which preserves the orientation.
Iteratively solving (2), a displacement flow $\Phi(s,x)$ can then be
constructed. $\Phi(0,x)=x$, $\Phi(t_{N})$ will also be a diffeomorphism that
preserves the orientation. We demonstrate the conservative nature of $T^{\\#}$
using the following examples.

######  Example 2.1.1.

Suppose $S$ denotes the density of a flow. We can associate $S$ to a
differential $n-$form: $\theta_{S}=SdV$, where $dV$ is the volume form on
$\Omega$. Then for any 1-1 and onto map $T$ that preserves the orientation of
$\Omega$,

$\displaystyle(T^{\\#}\theta_{S})(x)=((T^{-1})^{*}\theta_{S})(x)=S(T^{-1}(x))(T^{-1})^{*}(dV)=S(T^{-1}(x))\alpha(x)dV,$
(16)

in which $\alpha$ is some positive function. Thus
$(T^{\\#}S)(x)=S(T^{-1}(x))\alpha(x)$. We also have:

$\displaystyle\int_{\Omega}SdV=\int_{\Omega}(T^{-1})^{*}(SdV)=\int_{\Omega}T^{\\#}SdV.$
(17)

This means that the total mass is conserved.

######  Example 2.1.2.

Suppose that $S=u=(u^{1},u^{2})$ is the velocity field on a 2 dimensional
domain in $\mathbb{R}^{2}$. We associate to $S$ a differential $1-$form
$\theta_{S}=u^{1}dx^{1}+u^{2}dx^{2}$. Then for any 1-1 and onto map $T$ which
preserves the orientation of $\Omega$, we have

$\displaystyle(T^{\\#}\theta_{S})(x)=((T^{-1})^{*}\theta_{S})(x)$
$\displaystyle=$
$\displaystyle\Big{[}u^{1}(T^{-1}(x))\frac{\partial(T^{-1})^{1}}{\partial
x^{1}}+u^{2}(T^{-1}(x))\frac{\partial(T^{-1})^{2}}{\partial
x^{1}}\Big{]}dx^{1}+\Big{[}u^{1}(T^{-1}(x))\frac{\partial(T^{-1})^{1}}{\partial
x^{2}}+u^{2}(T^{-1}(x))\frac{\partial(T^{-1})^{2}}{\partial
x^{2}}\Big{]}dx^{2}$ (18)

This shows that

$\displaystyle
T^{\\#}(u^{1},u^{2})=\bigg{(}u^{1}(T^{-1}(x))\frac{\partial(T^{-1})^{1}}{\partial
x^{1}}+u^{2}(T^{-1}(x))\frac{\partial(T^{-1})^{2}}{\partial
x^{1}},u^{1}(T^{-1}(x))\frac{\partial(T^{-1})^{1}}{\partial
x^{2}}+u^{2}(T^{-1}(x))\frac{\partial(T^{-1})^{2}}{\partial x^{2}}\bigg{)}$
(19)

Since $d(T^{-1})^{*}=(T^{-1})^{*}d$, we have

$\displaystyle\int_{\Omega}d\theta_{S}=\int_{\Omega}(T^{-1})^{*}(d\theta_{S})=\int_{\Omega}d((T^{-1})^{*}\theta_{S})=\int_{\Omega}d(T^{\\#}\theta_{S})$
(20)

Note that $d\theta_{S}=(\frac{\partial u^{2}}{\partial x^{1}}-\frac{\partial
u^{1}}{\partial x^{2}})dx^{1}\wedge dx^{2}$. Therefore

$\displaystyle\int_{\Omega}(\frac{\partial u^{2}}{\partial
x^{1}}-\frac{\partial u^{1}}{\partial x^{2}})dx^{1}\wedge
dx^{2}=\int_{\Omega}(\frac{\partial[T^{\\#}(u^{1},u^{2})]^{2}}{\partial
x^{1}}-\frac{\partial[T^{\\#}(u^{1},u^{2})]^{1}}{\partial x^{2}})dx^{1}\wedge
dx^{2},$ (21)

in which $[T^{\\#}(u^{1},u^{2})]^{i}$ denotes the $i-$th component of
$[T^{\\#}(u^{1},u^{2})]$. Hence the total vorticity is conserved.

###  2.2 A physical interpretation of the tensor field definitions

Described above for 2.1.1 and 2.1.2, the conservative nature of $T^{\\#}$
largely depends on the choice of tensor fields. The question now translates on
how to choose the adequate tensor fields. Since Eq.(11) describes how the
tensor field changes along the vector field $u$, the Lie derivative in (11)
shall closely relate to the transport equation of physical fields.

######  Example 2.2.1.

Consider $\Omega\subset\mathbb{R}^{2}$ is a two dimensional domain, and the
original dynamical equation for $S$

$\displaystyle\frac{dS}{dt}:=S_{t}+v\cdot\nabla S=0,$ (22)

in which $v(t,x)$ is the true velocity field in the system. Define the flow
$\Phi(t,x)$:

$\displaystyle\begin{cases}\frac{\partial\Phi}{\partial
t}(t,x)=v(t,\Phi(t,x))\\\ \Phi(0,x)=x\end{cases}$ (23)

Then Eq.(22) means that $S(t,\Phi(t,x))=S(0,x)$ for any $t>0$, i.e. $S$ is
conserved along the trajectory of each point. In this case,
$\Phi(t):\Omega\longrightarrow\Omega$ plays the role of a “true" displacement
map. Recall that the goal is to find a displacement map $T$ for given
snapshots $S_{1}$ and $S_{2}$. Choose the tensor $\theta_{S}=S$, which is
simply a function (or differential 0-form as a tensor field), the relation
$S(t,\Phi(t,x))=S(0,x)$ writes $\theta_{S}(t)=\Phi(t)^{\\#}\theta_{S}(0)$.
Correspondingly, we have $\mathcal{L}_{u}\theta_{2}=u\cdot\nabla S_{2}$.

If the original dynamical equation for the physical field $S$ is modified

$\displaystyle\frac{dS}{dt}:=S_{t}+v\cdot\nabla S=\text{(forcing terms)}$ (24)

$S(0,x)\neq S(t,\Phi(t,x))$. However, we can still choose $\theta_{S}=S$. The
minimization problem (2) shall search for the vector field $u$ transporting
$S_{2}$ towards $S_{1}$ the fastest along the virtual flow $\Phi_{u}$. It is
also equivalent to look for a virtual vector field $u$ so that if $S_{2}$ is
transported along $u$ by the virtual system

$\displaystyle\begin{cases}\frac{\partial\tilde{S}}{\partial
s}+u\cdot\nabla\tilde{S}=0\\\ \tilde{S}(0,x)=S_{2}\end{cases}$ (25)

in a virtual time interval $[0,\epsilon]$, then $S_{2}$ transforms the most
rapidly towards $S_{1}$.

######  Example 2.2.2.

Now, consider $\Omega\subset\mathbb{R}^{n}$ an $n-$dimensional domain with the
original dynamical equation

$\displaystyle S_{t}+\nabla\cdot(Sv)=0,$ (26)

with $v(t,x)$ the true velocity field. Similar to the previous example, define
the flow $\Phi(t,x)$:

$\displaystyle\begin{cases}\frac{\partial\Phi}{\partial
t}(t,x)=v(t,\Phi(t,x))\\\ \Phi(0,x)=x\end{cases}$ (27)

and the differential $n-$form $\theta(t,x)=S(t,x)dx^{1}\wedge
dx^{2}\wedge\cdots\wedge dx^{n}$. Distinct from example 2.2.1, Eq.(26) implies
that $\Phi(t)^{*}\theta(t,x)=\theta(0,x)$, or equivalently,
$\theta(t,x)=(\Phi(t)^{-1})^{*}\theta(0,x)=\Phi(t)^{\\#}\theta(0,x)$. Indeed,
direct calculation yields:

$\displaystyle\Phi(t+\Delta t)^{*}\theta(t+\Delta t,x)-\Phi(t)^{*}\theta(t,x)$
(28) $\displaystyle=$ $\displaystyle\Phi(t)^{*}\Phi_{t\to t+\Delta
t}^{*}\theta(t+\Delta t)-\Phi(t)^{*}\theta(t)$ (29) $\displaystyle=$
$\displaystyle\Phi(t)^{*}\Big{[}(\Phi_{t\to t+\Delta t}^{*}-Id)\theta(t+\Delta
t)+\theta(t+\Delta t)-\theta(t)\Big{]}$ (30) $\displaystyle=$
$\displaystyle\Phi(t)^{*}\Big{[}(\Phi_{t\to t+\Delta
t}^{*}-Id)\theta(t)+\theta(t+\Delta t)-\theta(t)\Big{]}+o(\Delta t)$ (31)
$\displaystyle=$ $\displaystyle\Delta
t\Phi(t)^{*}\Big{[}\mathcal{L}_{v(t)}\theta(t)+\theta_{t}(t)\Big{]}+o(\Delta
t)$ (32) $\displaystyle=$ $\displaystyle\Delta
t\Phi(t)^{*}\Big{[}\Big{(}S_{t}(t,x)+\nabla\cdot(S(t,x)v(t))\Big{)}dx^{1}\wedge\cdots\wedge
dx^{n}\Big{]}+o(\Delta t)$ (33) $\displaystyle=$ $\displaystyle o(\Delta t)$
(34)

It shows that $\Phi(t)^{*}\theta(t,x)$ does not change with time.

In both examples, the choice of tensor $\theta$ is determined by the original
transport equation (22) or (26). This choice then implies $T^{\\#}$, and the
minimization problem (2) will search for a vector field $u$ that most rapidly
transports $S_{2}$ towards $S_{1}$ by a virtual flow $\Phi_{u}$. For
geophysical fields, $T^{\\#}$ is thus essential for the application of this
method, and further discussed in section 3.

###  2.3 Comparison with OF algorithms

Given a time series of snapshots (i.e. brightness fields in
$\Omega\subset\mathbb{R}^{2}$) $\\{...,S(t-\Delta t),S(t),S(t+\Delta
t),...\\}$, the [1] OF algorithm aims at recovering the true velocity field
$u(t,x)$ by minimizing the following functional (with proper boundary
conditions):

$\displaystyle u(t)=\arg\min_{u}\|S_{t}+\langle\nabla
S,u\rangle\|^{2}+a_{1}\|\nabla
u\|^{2}=\operatorname*{arg\,min}_{u}\int_{\Omega}[S_{t}+\langle\nabla
S,u\rangle]^{2}+a_{1}|\nabla u|^{2}d^{n}x.$ (35)

We remark that the first term implies that the hidden $T^{\\#}$ in (35) is
simply $T^{\\#}S=S\circ T^{-1}$ for the brightness field $S$. Or equivalently,
$T^{\\#}S$ is inferred from $T^{\\#}\theta_{S}$, where $\theta_{S}=S$ is a
differential $0-$form. Indeed, for any vector field $u$, let
$\Phi_{u}(s,x):[0,\Delta t]\times\Omega\longrightarrow\Omega$ be the following
flow of the points in the domain:

$\displaystyle\begin{cases}\frac{\partial\Phi_{u}(s,x)}{\partial
s}=u(\Phi_{u}(s,x))\\\ \Phi_{u}(0,x)=x.\end{cases}$ (36)

The first term in (35) represents the material derivative of
$S(t+s,\Phi_{u}(s,x))$ with respect to $s$ at $s=0$. The optimisation problem
(35) searches for vector field ${\bf u}$ so that $S(t,x)\approx S(t+\Delta
t,\Phi_{u}(\Delta t,x))$. Let $S_{2}(x)=S(t,x)$, $S_{1}(x)=S(t+\Delta t,x)$,
and $T(x)=\Phi_{u}(\Delta t,x)$. Hence the first term in (35) searches for the
vector field $u$ so that $S_{1}(x)\approx S_{2}(T^{-1}(x))$. Clearly, this
implies that $T^{\\#}S_{2}=S_{2}\circ T^{-1}$. [16] generalizes this original
OF method to the Riemannian manifold context. Still, the formulation of the
problem in [16] implies that $T^{\\#}S_{2}=S_{2}\circ T^{-1}$.

In [29], an other variant is proposed, to minimize

$\displaystyle
u(t)=\operatorname*{arg\,min}_{u}\|S_{t}+\nabla\cdot(Su)\|^{2}+a_{1}\|\nabla
u\|^{2}=\operatorname*{arg\,min}_{u}\int_{\Omega}[S_{t}+\nabla\cdot(Su)]^{2}+a_{1}|\nabla
u|^{2}d^{n}x.$ (37)

Let $S_{1}(x)=S(t+\Delta t,x)$ and $S_{2}(x)=S(t,x)$. Following the analysis
in example 2.2.2, Eq.(37) aims at finding the vector field $u$ so that
$\Phi_{u}^{\\#}S_{2}\approx S_{1}$, in which $\Phi_{u}^{\\#}S$ is inferred
from $\Phi_{u}^{\\#}\theta_{S}$ for $\theta_{S}=Sdx^{1}\wedge\cdots\wedge
dx^{n}$. Thus, it can be concluded that the method of [29] implies that
$T^{\\#}$ should be defined as if the physical fields are associated to
differential $n-$forms. This is explicitly stated in [28] and corresponds to
physically-constrained choices of tensor fields in the present framework. The
proper choices of tensor fields shall indeed enable to also transport
unobserved state variables, consistent with the underlying dynamics. OF
algorithms can thus be reformulated for snapshots corresponding to snapshots
of tensor fields. For a time series $\\{\cdots,\theta(t-\Delta
t),\theta(t),\theta(t+\Delta t),\cdots\\}$ of snapshots of tensor fields, the
generalized OF method is to find the vector field $u$:

$\displaystyle
u(t)=\operatorname*{arg\,min}_{u}\int_{\Omega}|\theta_{t}+\mathcal{L}_{u}\theta|_{x}^{2}+a_{1}|\nabla
u|_{x}^{2}dV,$ (38)

in which $u$ satisfies proper boundary conditions. Here we provide a theorem
for the existence of the OF solution for Dirichlet boundary condition. Note
that in [1], $u$ is only required to have zero normal component at the
boundary.

######  Theorem 2.2.

For fixed $t$, assume that both tensor fields $\theta$ and $\theta_{t}$ have
finite $H^{2}$ norm and the given vector field $v$ on $\partial\Omega$ also
has finite $H^{1}$ norm. Further assume that
$|d\theta|_{x},|\delta\theta|_{x}$ are bounded in $\Omega$. Then for compact
oriented Riemannian manifold $\Omega$ without boundary, the following
minimisation problem has a unique solution which has a finite $H^{1}$ norm:

$\displaystyle u(t)=\operatorname*{arg\,min}_{u\in
H^{1}(T\Omega)}\int_{\Omega}|\theta_{t}+\mathcal{L}_{u}\theta|_{x}^{2}+a_{0}|\theta_{u}|_{x}^{2}+a_{1}|d\theta_{u}|_{x}^{2}+a_{1}|\delta\theta_{u}|_{x}^{2}dV.$
(39)

For a compact oriented Riemannian manifold with boundary, the following
minimisation problem has a unique solution which has a finite $H^{1}$ norm:

$\displaystyle\begin{cases}u(t)=\displaystyle\operatorname*{arg\,min}_{u\in
H^{1}(T\Omega)}\int_{\Omega}|\theta_{t}+\mathcal{L}_{u}\theta|_{x}^{2}+a_{1}|d\theta_{u}|_{x}^{2}+a_{1}|\delta\theta_{u}|_{x}^{2}dV\\\
u(t)\big{|}_{\partial\Omega}=v.\end{cases}$ (40)

The existence, uniqueness, and smoothness of order 1 of the solution is a
direct consequence of Riesz representation (or Lax-Milgram) theorem. A
mathematical proof is provided in the appendix.

While the proposed algorithm can be applied whenever two snapshots are
provided, the OF algorithms require that the time between two consecutive
images to be small, to best approximate $\theta_{t}$. To overcome such a
constrain, the standard Horn and Schunck algorithm can be incorporated within
a multi-scale strategy [18], and wavelets [8]. Note, an infinitesimal
formulation of the Horn and Schunck algorithm in the Euclidean space is
proposed in [8].

###  2.4 Comparison with other methods

[25] proposes a variational method to find the optimal way to transport a
point $S_{2}\in\mathcal{M}$ to another point $S_{1}\in\mathcal{M}$ at a
minimal cost. Here $\mathcal{M}$ is a Riemannian manifold. For geophysical
applications, this manifold $\mathcal{M}$ can be interpreted as the
configuration space. Hence each point $S_{i}\in\mathcal{M}$ represents a state
vector. It can be interpreted in a way that, to transport $S_{2}$ to $S_{1}$,
[25] does not only use vector fields, but also an external forcing to
influence the state at each time step. The solution of the variational problem
then consists of a time sequence of vector fields and a time sequence of
"external forcing". Considering an "external forcing", the eventual state
$T^{\\#}S_{2}$ can exactly match the target state $S_{1}$. This method follows
from rigorous mathematical developments [24], but turns out to be very
computationally demanding.

Some efforts about large deformation diffeomorphic metric matching (LDDMM) [3,
4] are examples of directly solving (1) to obtain a diffeomorphism $T$
belonging to the same connected component as the identity map in the group of
diffeomorphisms $\mathcal{D}(\Omega)$. Since an external forcing is not
considered, these methods can be stated to be simplified versions of that
proposed by [25]. The regularization term in (1) is chosen to be

$\displaystyle\|T\|^{2}=\int_{0}^{1}\|v(t)\|_{V}^{2}dt$ (41)

for some prescribed norm $V$. And
$T(x)\approx(Id+1/Nv_{\frac{N-1}{N}})\circ\cdots\circ(Id+1/Nv_{0})$ as
described in the introduction. Note that the optimisation problem in [3, 4] is
a specific case of Eq.(6) in [26] and Eq.(1) in [23]. Despite the fact that
$T^{\\#}$ in [3, 4, 26, 23, 25] is not particularly designed for geophysical
fields and sometimes the boundary condition is ignored, the framework proposed
in this manuscript can be stated to be a simplified version of LDDMM.

###  2.5 About the boundary condition of (13)

Boundary condition is necessary in the optimisation problem (13). In fact, the
boundary condition can be obtained by solving for (12) or (13) based on the
data on $\partial\Omega$. The process can be illustrated using the following
idealistic example.

Suppose that $\Omega$ is a three dimensional ocean, the boundary of which
consists of two parts $\Omega_{b}$ and $\Omega_{s}$, in which $\Omega_{b}$
refers to the ocean basin (i.e. the land boundaries), and $\Omega_{s}$ the sea
surface. To determine the boundary condition of $u$ on $\Omega_{b}$ and
$\Omega_{s}$, it is first natural to set $u\big{|}_{\Omega_{b}}=0$. On
$\Omega_{s}$, we must first solve for (13) where $\theta_{1}$ and $\theta_{2}$
are tensor fields on $\Omega_{s}$. Boundaries of $\Omega_{s}$ coincide with
the coast lines, and natural to set boundary conditions to be $0$. The domain
$\Omega_{s}$ is a sub-manifold of the sphere. Hence the “Riemannian" context
is necessary for this example. The complete process at each iterative step is
illustrated in Algorithm 1.

Data: tensor fields $\theta_{1}$ and $\theta_{2}$ on $\Omega$; tensor fields
$\alpha_{1},\alpha_{2}$ on $\Omega_{s}$; small positive number $\epsilon>0$;
the total number of iterations $N$; $dV_{s}$ the volume form on $\Omega_{s}$;
$dV$ the volume form on $\Omega$.

Result: vector fields $u_{1},...,u_{N}$

for  _$i=1,2,...,N$_ do

$u_{s,i}\leftarrow\displaystyle\operatorname*{arg\,min}_{u\in
T\Omega_{s}}\int_{\Omega_{s}}\bigg{[}-2\langle\alpha_{1}-\alpha_{2},\mathcal{L}_{u}\alpha_{2}\rangle_{x}+|d\theta_{u}|_{x}^{2}+|\delta\theta_{u}|_{x}^{2}\bigg{]}dV_{s}$
with boundary condition $u_{s,i}\big{|}_{\partial\Omega_{s}}=0$;

$u_{i}\leftarrow\displaystyle\operatorname*{arg\,min}_{u\in
T\Omega}\int_{\Omega}\bigg{[}-2\langle\theta_{1}-\theta_{2},\mathcal{L}_{u}\theta_{2}\rangle_{x}+|d\theta_{u}|_{x}^{2}+|\delta\theta_{u}|_{x}^{2}\bigg{]}dV$
with boundary condition $u_{i}\big{|}_{\Omega_{b}}=0$,
$u_{i}\big{|}_{\Omega_{s}}=u_{s,i}$;

$\theta_{2}\leftarrow\theta_{2}+\epsilon\mathcal{L}_{u_{i}}\theta_{2}$;

$\alpha_{2}\leftarrow\alpha_{2}+\epsilon\mathcal{L}_{u_{s,i}}\alpha_{2}$

end for

Algorithm 1 How to determine the boundary condition at each iterative step (an
example)

###  2.6 The case for localized observations

In practice, it is a common situation that some state variables are solely
observed in a subdomain $\Omega_{1}\subset\Omega$, instead of in the full
domain $\Omega$ (i.e. over a satellite swath). A weight function $W$ can be
constructed and the first terms in Eq.(12) and (13) replaced by:

$\displaystyle\frac{\partial|\theta_{1}-\Phi_{u}(s)^{\\#}\theta_{2}|_{x}^{2}}{\partial
s}\big{|}_{s=0}W(x),$ (42)

in which the pre-chosen $W(x)$ takes value $1$ for a majority of points inside
$\Omega_{1}$ but decreases to $0$ smoothly as $x$ approaches the boundary of
$\Omega_{1}$. Then Eq.(14) changes to

$\displaystyle u=\operatorname*{arg\,min}_{u\in
T\Omega}\int_{\Omega}\bigg{[}-2\langle\theta_{1}-\theta_{2},\mathcal{L}_{u}\theta_{2}\rangle_{x}W(x)+a_{1}|d\theta_{u}|_{x}^{2}+a_{1}|\delta\theta_{u}|_{x}^{2}+a_{0}|\theta_{u}|_{x}^{2}\bigg{]}dV,$
(43)

and Eq.(15) changes to

$\displaystyle\begin{cases}u=\displaystyle\operatorname*{arg\,min}_{u\in
T\Omega}\int_{\Omega}\bigg{[}-2\langle\theta_{1}-\theta_{2},\mathcal{L}_{u}\theta_{2}\rangle_{x}W(x)+a_{1}|d\theta_{u}|_{x}^{2}+a_{1}|\delta\theta_{u}|_{x}^{2}\bigg{]}dV\\\
u\big{|}_{\partial\Omega}=i_{*}v.\end{cases}$ (44)

The theorem of existence and uniqueness of solution still holds for
minimisation problems (43) and (44) following similar arguments given in the
appendix.

##  3 Towards a new nudging strategy and its application in data assimilation

###  3.1 Methodology

Let $S$ be the full state variable, $Y=h(S)$ some state variable derived from
$S$. Assume that we have the model estimate $S_{\text{model}}$, and that $Y$
is fully observed on the domain. Our target is to

  * (1), derive a displacement flow $\Phi$ so that $\Phi^{\\#}Y^{\text{model}}\approx Y^{\text{obs}}$;

  * (2), apply $\Phi$ to $S^{\text{model}}$ to correct the displacement of the full state variable $S$.

From section 2, we need to separately find the tensor fields for $Y$ and $S$.
Consistent definitions should be determined by the dynamical equation of the
system, leading to the determination of $T^{\\#}$, hence the explicit
formulation of Eq.(14) (or (15)).

Already indicated in [2] and [20], this nudging strategy can be incorporated
with ensemble Kalman filter. The basic idea is to correct the displacement of
each ensemble member before applying EnKF. Methods reported in [2] and [20]
use different cost functions to estimate the displacements. Both methods
implicitly assume $T^{\\#}S=S\circ T$ or $S\circ T^{-1}$. However, see Fig.2,
$T^{\\#}$ without considering the dynamics could possibly destroy the
intrinsic feature of the unobserved physical fields.

To demonstrate the advantage to consistently constrain $T^{\\#}$, the
following simple version of morphed EnKF is compared numerically with the
plain EnKF algorithm:

  * (1), choose tensor fields for the observed state variable $Y$ and the full state variable $X$;

  * (2), find the displacement flow $\Phi_{i}$ according to the observation $Y^{\text{obs}}$ and the model estimates $Y^{\text{model}}_{i}$ from the $i-$th ensemble member;

  * (3), apply $\Phi_{i}$ to the full state vector of the $i-$th ensemble member: $X^{f}_{i}\leftarrow\Phi_{i}^{\\#}X^{f}_{i}$;

  * (4), apply plain EnKF to the updated ensemble $\\{X^{f}_{i}\>i=1,...,N_{e}\\}$.

A morphed ensemble member based on the plain $T^{\\#}S=S\circ T^{-1}$ is also
presented for comparison.

###  3.2 Numerical results

To illustrate our purpose, a data assimilation experiment is conducted for
only one time step, using the thermal shallow water equation [27]. This model
consists of three state variables: $h-$the water height, $v=(v^{1},v^{2})-$the
velocity field, and $\Theta-$the buoyancy (or density contrast):

$\displaystyle\frac{\partial h}{\partial t}+\nabla\cdot(hv)$
$\displaystyle=0,$ (45) $\displaystyle\frac{\partial\Theta}{\partial
t}+(v\cdot\nabla)\Theta$ $\displaystyle=-\kappa(h\Theta-h_{0}\Theta_{0}),$
(46) $\displaystyle\frac{\partial v}{\partial
t}+(v\cdot\nabla)v+f\hat{z}\times v$
$\displaystyle=-\nabla(h\Theta)+\frac{1}{2}h\nabla\Theta.$ (47)

Both $\Theta$ and $h$ are assumed strictly positive at each point. It is
assumed that both the absolute vorticity $\omega=\frac{\partial
v^{2}}{\partial x}-\frac{\partial v^{1}}{\partial y}$ and $h$ are fully
observed, while $\Theta$ is completely unobserved. The details of the
experimental setup can be found in the appendix. A Python code to completely
reproduce the numerical results is provided at 10.5281/zenodo.10252176.

Figure 3: Comparison of the target fields (the first row), the original fields
(the third row), and the morphed fields (the second row).  Figure 4: Another
example of Fig.3

Fig.3 and Fig.4 show two examples of the (interpolated) truth (the target
fields), the prior estimate of one member (the original fields), and the
morphed prior estimate (the morphed fields). Expected, the phase of the
original vortex can be adjusted, to some extent, but not perfectly. This is
not surprising because the cost function in the optimisation problem (12) (or
(13)) has two parts. The first part is derived from the dynamics of the
original system. But the second part is merely for mathematical reasons. Thus,
the virtual displacement flow $\Phi$ cannot transform $\Theta_{2}$ to exactly
match $\Theta_{1}$. But the correct choice of tensor fields can maintain the
dynamical balance of the three fields during the morphing process to some
extent. Looking at Fig.3 and 4 more closely, it is found that the difference
between the morphed $h_{2}$ field and the target $h_{1}$ field is much smaller
than that between the morphed $\omega_{2}$ and the target $\omega_{1}$. We
don’t have a theoretical explanation of this fact. We suspect that it can be
attributed to the positiveness of the $h$ field.

To demonstrate the superiority of $T^{\\#}$, the nudging process without
introducing the concept of tensor fields is also considered and tested. In
this case, $T^{\\#}S=S\circ T^{-1}$, independent of $S=h$, $\omega$, $\Theta$,
or $v^{1},v^{2}$. The initial value of $\omega_{2}$ $h_{2}$, etc. are taken to
be the corresponding fields of one of the ensemble priors. We first run the
morphing process for $N=10000$ time steps. It is found that, for this specific
member, the mean-squared-error (MSE) of $h$ is always decreasing, while the
MSE for $\omega$ decreases till $N\approx 6000$ and then starts to increase.
The fields at $N=6000$ are plotted in Fig.5. Apparently, the morphed
$\omega_{2}$ has lost its intrinsic feature.

Figure 5: The target field (first row), the original field (third row), and
the field morphed by $T^{\\#}S=S\circ T^{-1}$ for
$S=h,\omega,\Theta,v^{1},v^{2}$ (second row). Figure 6: The truth (first row),
the prior mean of ensemble (the second row), and the prior mean of the morphed
ensemble (the third row).

Fig. 6 presents the truth, the prior mean of the ensemble, and the prior mean
of the morphed ensemble. The direct prior mean has completely lost the small
scale structures inside the vortex, while the prior mean of the morphed
ensemble still maintain the small scale features to some extent. Stated in the
introduction, our interpretation is that $y^{o}-x^{b}$ is no longer a good
representation of the error in this case. Instead, the location error, or more
generally the displacement flow, can better represent the error of each
ensemble member. Thus the space of ensemble members is not a linear subspace
in the Euclidean space of the state vector, but could be a curved manifold. In
this case, it is not surprising to see that the arithmetic mean lies outside
of the manifold. To address this problem, the Fréchet mean instead of the
arithmetic mean should be used to define the ensemble mean.

Figure 7: The truth (first row), the posterior estimate of one member (the
second row), and the posterior estimate of one morphed member (the third row).

Fig.7 presents the posterior estimate of one ensemble member using EnKF and
morphed EnKF (denoted by mEnKF for short). Apparently, the plain EnKF results
in erroneous estimates, while the mEnKF still produces reasonable estimates
with fine scale features inside the vortex. Again, this is the disadvantage of
linear algorithms when the ensemble members do not lie on a linear subspace of
the Euclidean space of state vectors.

##  4 Conclusion

A complete framework is introduced to estimate the displacement flow between
two tensor fields on a compact oriented Riemannian manifold with or without
boundary. It is proved that the solution always exists and is unique. When a
time sequence of snapshots of geophysical fluids is available, a reformulation
of the generalized OF algorithm [28] can also be obtained from this
differential geometry perspective.

The novelty of the proposed framework lies in the definition of $T^{\\#}S$ for
a given displacement map $T$ (not necessarily the displacement generated by
the true velocity field) acting on a given state variable $S$: $T^{\\#}S$ must
be determined by the original dynamics of $S$. Technically, each state
variable $S$ is associated to a tensor field $\theta_{S}$, and $T^{\\#}S$ is
determined by $T^{\\#}\theta_{S}$. A key advantage is to consistently correct
the displacement of some state variable while best maintaining the intrinsic
structure of the underline physical fields.

This new framework can then be used as a nudging process when only part of the
state variables are observed. It can also be used to correct the displacement
of each ensemble member before applying the linear EnKF algorithm. Numerical
results with a double-vortex model show that the morphed EnKF algorithm
produces more reliable posterior estimates than the plain EnKF algorithm. The
preferable usage of $T^{\\#}$ is well demonstrated by comparing fields morphed
by algorithms derived or not from the concepts of differential geometry.

It must be pointed out that all physical fields could not always be associated
to a tensor field. The use of a tensor field is only required to define the
correct $T^{\\#}$. For an arbitrary derived physical field $S$, $T^{\\#}S$
should be determined by the forward operator $S=F(X)$, where $X$ is the full
state vector. Thus, for a specific derived variable $S$, whether a
displacement vector field can be determined by the optimisation procedure (12)
or (13) needs to be further studied.

## Acknowledgments

The authors would like to express their gratitude to Long Li at the department
of mathematics of ShanghaiTech University and Wei Pan at department of
mathematics of Imperial College London for their insightful and helpful
discussion, and and Etienne Mémin for his introduction on relevant works. YZ’s
research is supported by Hohai university (Grant No. 1074-423101 and
522-020512). VR’s research is supported by the company SCALIAN DS and by
France Relance through the MORAANE project. BC’s research is supported by ERC
EU SYNERGY Project No. 856408-STUOD, and the the support of the ANR Melody
project.

##  Appendix A Proof of theorem 2.1

Recall that $d$ is the exterior derivative operator, $*$ the Hodge star
operator on differential forms, and $\delta=(-1)^{n(k-1)+1}*d*$ is the so-
called co-differential operator. Let $i:\partial\Omega\rightarrow\Omega$ be
the natural embedding of $\partial\Omega$. For convenience but without loss of
generality, we set the regularization parameter $a_{0}=a_{1}=1$.

First, consider the case when $\Omega$ is a compact oriented Riemannian
manifold with boundary. Assume that $u$ is a solution of (15). Then for any
vector field $h\in T\Omega$ such that $h\big{|}_{\partial\Omega}=0$, we have

$\displaystyle\frac{\partial}{\partial\epsilon}\Big{[}\int_{\Omega}-2\langle\theta_{1}-\theta_{2},\mathcal{L}_{u+\epsilon
h}\theta_{2}\rangle_{x}+|d\theta_{u+\epsilon
h}|_{x}^{2}+|\delta\theta_{u+\epsilon
h}|_{x}^{2}dV\Big{]}\Big{|}_{\epsilon=0}=0$ (48)

This implies that for any $h\in T\Omega$ such that
$h\big{|}_{\partial\Omega}=0$,

$\displaystyle\int_{\Omega}-2\langle\theta_{1}-\theta_{2},\mathcal{L}_{h}\theta_{2}\rangle_{x}+\langle
d\theta_{u},d\theta_{h}\rangle_{x}+\langle\delta\theta_{u},\delta\theta_{h}\rangle_{x}dV=0$
(49)

Applying the Green’s formula for differential forms (see for instance Eq.(2.1)
in [9]), we have that

$\displaystyle\int_{\Omega}\langle
d\theta_{h},d\theta_{u}\rangle_{x}dV-\int_{\Omega}\langle\theta_{h},\delta
d\theta_{u}\rangle_{x}dV=\int_{\partial\Omega}(i^{*}\theta_{h})\wedge*(i^{*}d\theta_{u})=0$
(50)

since $\theta_{h}\big{|}_{\partial\Omega}=0$. Similarly,

$\displaystyle\int_{\Omega}\langle
d\delta\theta_{u},\theta_{h}\rangle_{x}dV-\int_{\Omega}\langle\delta\theta_{u},\delta\theta_{h}\rangle_{x}dV=\int_{\partial\Omega}(i^{*}\delta\theta_{u})\wedge*(i^{*}\theta_{h})=0.$
(51)

Thus

$\displaystyle\int_{\Omega}\langle
d\theta_{u},d\theta_{h}\rangle_{x}+\langle\delta\theta_{u},\delta\theta_{h}\rangle
dV=\int_{\Omega}\langle\theta_{h},(d\delta+\delta
d)\theta_{u}\rangle_{x}dV=\int_{\Omega}\langle\theta_{h},\Delta\theta_{u}\rangle_{x}dV,$
(52)

where $\Delta=d\delta+\delta d$ is the Hodge Laplacian operator.

Let $\alpha=\theta_{1}-\theta_{2},\beta=\theta_{2}$. By partition of unity,
$\langle\alpha,\mathcal{L}_{h}\beta\rangle$ can be decomposed into a finite
sum:

$\displaystyle\langle\alpha,\mathcal{L}_{h}\beta\rangle=\sum_{i}\langle\alpha_{i},\mathcal{L}_{h}\beta_{i}\rangle,$
(53)

in which $\alpha_{i}$’s and $\beta_{i}$’s have the same regularity as $\alpha$
and $\beta$, and are all of the form $\eta_{1}\otimes\cdots\otimes\eta_{l}$,
with $\eta_{j}\in V_{j}$ and $V_{j}=H^{1}(T\Omega)$ or $H^{1}(T^{*}\Omega)$.
Note that

$\displaystyle\mathcal{L}_{h}(\eta_{1}\otimes\cdots\otimes\eta_{l})=\mathcal{L}_{h}\eta_{1}\otimes\eta_{2}\otimes\cdots\otimes\eta_{l}+\cdots+\eta_{1}\otimes\cdots\otimes\eta_{l-1}\otimes\mathcal{L}_{h}\eta_{l},$
(54)

Thus there exists, finitely, many known numbers $m_{i}$ and known
$\xi_{i},\eta_{i}\in T\Omega$ or $\xi_{i},\eta_{i}\in T^{*}\Omega$, so that

$\displaystyle\langle\alpha,\mathcal{L}_{h}\beta\rangle=\sum_{i}m_{i}\langle\xi_{i},\mathcal{L}_{h}\eta_{i}\rangle.$
(55)

Since for any $\xi,\eta\in T\Omega$,
$\langle\xi,\mathcal{L}_{h}\eta\rangle=-\langle\theta_{\eta},\mathcal{L}_{h}\theta_{\xi}\rangle$,
in which $\theta_{\xi}$ (or $\theta_{\eta}$) refers to the differential 1-form
such that $\theta_{\xi}(v)=\langle\xi,v\rangle$ (or
$\theta_{\eta}(v)=\langle\eta,v\rangle$ resp. ) for any vector field $v\in
T\Omega$, without loss of generality we can assume that all the $\xi_{i}$’s
and $\eta_{i}$’s in Eq. (55) are differential 1-forms. With Cartan’s formula
$\mathcal{L}_{h}=d\iota_{h}+\iota_{h}d$ and Stokes’ formula, we have

$\displaystyle\int_{\Omega}-2\langle\theta_{1}-\theta_{2},\mathcal{L}_{h}\theta_{2}\rangle_{x}dV=\int_{\Omega}-2\langle\alpha,\mathcal{L}_{h}\beta\rangle_{x}dV=-2\int_{\Omega}\sum_{i}m_{i}\langle\xi_{i},\mathcal{L}_{h}\eta_{i}\rangle_{x}dV$
$\displaystyle=$
$\displaystyle-2\int_{\Omega}\sum_{i}m_{i}\mathcal{L}_{h}\eta_{i}\wedge*\xi_{i}=-2\sum_{i}\int_{\Omega}[(d\iota_{h}+\iota_{h}d)\eta_{i}]\wedge*(m_{i}\xi_{i})$
$\displaystyle=$
$\displaystyle-2\sum_{i}\int_{\Omega}d[\iota_{h}\eta_{i}\wedge*(m_{i}\xi_{i})]-\iota_{h}\eta_{i}\wedge
d*(m_{i}\xi_{i})+\iota_{h}d\eta_{i}\wedge*(m_{i}\xi_{i})$ $\displaystyle=$
$\displaystyle-2\Big{\\{}\sum_{i}\int_{\partial\Omega}\iota_{h}\eta_{i}\wedge*(m_{i}\xi_{i})+\int_{\Omega}-\iota_{h}\eta_{i}\wedge
d*(m_{i}\xi_{i})+\iota_{h}d\eta_{i}\wedge*(m_{i}\xi_{i})\Big{\\}}$
$\displaystyle=$ $\displaystyle 2\sum_{i}\int_{\Omega}\iota_{h}\eta_{i}\wedge
d*(m_{i}\xi_{i})+\iota_{h}d\eta_{i}\wedge*(m_{i}\xi_{i})$ (56)

For any differential 1-form $\theta$, denote by $X_{\theta}\in T\Omega$ the
vector field such that $\langle\eta,\theta\rangle=\eta(X_{\theta})$ for any
$\eta\in T^{*}\Omega$. Then

$\displaystyle\iota_{h}d\eta_{i}\wedge*(m_{i}\xi_{i})=\langle\iota_{h}d\eta_{i},m_{i}\xi_{i}\rangle_{x}dV$
$\displaystyle=$
$\displaystyle(\iota_{h}d\eta_{i})(X_{m_{i}\xi_{i}})dV=d\eta_{i}(h,X_{m_{i}\xi_{i}})dV=-d\eta_{i}(X_{m_{i}\xi_{i}},h)dV$
$\displaystyle=$
$\displaystyle-(\iota_{X_{m_{i}\xi_{i}}}d\eta_{i})(h)dV=-\langle\iota_{X_{m_{i}\xi_{i}}}d\eta_{i},\theta_{h}\rangle_{x}dV,$
(57)

and

$\displaystyle\iota_{h}\eta_{i}\wedge
d*(m_{i}\xi_{i})=\langle\iota_{h}\eta_{i},*d*m_{i}\xi_{i}\rangle_{x}dV=\langle\theta_{h}(X_{\eta_{i}}),*d*(m_{i}\xi_{i})\rangle_{x}dV$
$\displaystyle=$
$\displaystyle\theta_{h}(X_{\eta_{i}})*d*(m_{i}\xi_{i})dV=\langle\theta_{h},[*d*(m_{i}\xi_{i})]\eta_{i}\rangle
dV$ (58)

It thus exists a differential 1-form $\mu\in L^{2}(T^{*}\Omega)$, which is
determined by $\eta_{i},m_{i},\xi_{i}$, so that

$\displaystyle\int_{\Omega}2\langle\theta_{1}-\theta_{2},\mathcal{L}_{h}\theta_{2}\rangle_{x}dV=\int_{\Omega}\langle\theta_{h},\mu\rangle_{x}dV.$
(59)

Combining Eq.(52) and (59), Eq.(49) is equivalent to

$\displaystyle\int_{\Omega}-\langle\theta_{h},\mu\rangle_{x}+\langle\theta_{h},\Delta\theta_{u}\rangle_{x}dV=0.$
(60)

Thus the optimisation problem (13) is equivalent to solving the following
equation for $\theta_{u}$:

$\displaystyle\begin{cases}\Delta\theta_{u}=\mu\\\
\theta_{u}\big{|}_{\partial\Omega}=\theta_{\widetilde{i_{*}v}}\big{|}_{\partial\Omega},\end{cases}$
(61)

in which $\widetilde{i_{*}v}$ is a smooth extension of $v$ from
$\partial\Omega$ to $\Omega$. Then the existence, uniqueness, and smoothness
of the solution to (61) is then guaranteed by theorem (3.4.10) of [22].

For the case, $\Omega$ a compact and oriented Riemannian manifold without
boundary, the optimisation problem (12) is equivalent to solving the following
equation for differential 1-forms:

$\displaystyle(1-\Delta)\theta_{u}=\mu,$ (62)

where $\mu$ is a differential 1-form that can be derived from the given data.
The spectral theory of Laplacian operator on a Riemannian manifold (see for
instance theorem (1.30) and (1.31) in [21]) states that the space of square
integrable differential $k-$forms $W^{0,2}(\bigwedge^{k}T^{*}\Omega)$ has an
orthonormal basis
$\\{\phi_{i}:\Delta\phi_{i}=\lambda_{i}\phi_{i},\displaystyle\int_{\Omega}\langle\phi_{i},\phi_{i}\rangle_{x}dV=1,0\leq
i<\infty\\}$, and that all eigen-forms $\phi_{i}$ are smooth on $\Omega$. Note
that $\lambda_{i}\leq 0$ $\forall i$. Thus we can assume the decomposition:
$\mu=\sum_{i\geq 0}a_{i}\phi_{i}$. Then
$\theta_{u}=\displaystyle\sum_{i}\frac{a_{i}}{1-\lambda_{i}}\phi_{i}$ is a
solution to Eq.(62). The identity can be verified directly. But in addition we
need to show that $\theta_{u}\in
W^{2,2}=\\{\theta:\displaystyle\int_{\Omega}\langle
d\delta\theta,d\delta\theta\rangle_{x}+\langle\delta d\theta,\delta
d\theta\rangle_{x}dV<\infty\\}$. First we show that $\theta_{u}\in
W^{1,2}=\\{\theta:\displaystyle\int_{\Omega}\langle
d\theta,d\theta\rangle_{x}+\langle\delta\theta,\delta\theta\rangle_{x}dV<\infty\\}$.
This can be verified directly:

$\displaystyle\int_{\Omega}\langle
d\theta_{u},d\theta_{u}\rangle_{x}+\langle\delta\theta_{u},\delta\theta_{u}\rangle_{x}dV=\lim_{N\to\infty}\sum_{i=0}^{N}\frac{a_{i}^{2}}{1-\lambda_{i}}\int_{\Omega}\langle
d\phi_{i},d\phi_{i}\rangle_{x}+\langle\delta\phi_{i},\delta\phi_{i}\rangle_{x}dV$
$\displaystyle=$
$\displaystyle\lim_{N\to\infty}\sum_{i}\frac{a_{i}^{2}}{1-\lambda_{i}}\int_{\Omega}\langle\delta
d\phi_{i},\phi_{i}\rangle_{x}+\langle
d\delta\phi_{i},\phi_{i}\rangle_{x}dV=\lim_{N\to\infty}\sum_{i}\frac{\lambda_{i}a_{i}^{2}}{1-\lambda_{i}}\leq\sum_{i}|a_{i}|^{2}$
(63)

This shows that $\theta_{u}\in W^{1,2}$. Next, direct calculation yields that
$\displaystyle\int_{\Omega}\langle\Delta\theta,\Delta\theta\rangle_{x}dV=\int_{\Omega}\langle
d\delta\theta,d\delta\theta\rangle_{x}+\langle\delta d\theta,\delta
d\theta\rangle_{x}dV$ for any $\theta$. Hence

$\displaystyle\int_{\Omega}\langle
d\delta\theta_{u},d\delta\theta_{u}\rangle_{x}+\langle\delta
d\theta_{u},\delta
d\theta_{u}\rangle_{x}dV=\int_{\Omega}\langle\Delta\theta_{u},\Delta\theta_{u}\rangle_{x}dV=\sum_{i}\frac{\lambda_{i}^{2}}{(1-\lambda_{i})^{2}}a_{i}^{2}\leq\sum_{i}|a_{i}|^{2}.$
(64)

This shows that $(1-\Delta)^{-1}\mu$ is a well-defined twice-differentiable
$1-$form.

##  Appendix B Proof of theorem 2.2

First assume that the manifold has $C^{1}-$boundary. Then the boundary
condition $v$ can be extended to $u_{0}\in H^{1}(\Omega)$ such that
$u_{0}\big{|}_{\partial\Omega}=v$. Thus the optimisation problem (40) is
equivalent to

$\displaystyle\begin{cases}u(t)=\displaystyle\operatorname*{arg\,min}_{u\in
H^{1}(T\Omega)}\int_{\Omega}|\theta_{t}+\mathcal{L}_{u_{0}}\theta+\mathcal{L}_{u}\theta|_{x}^{2}+|d\theta_{u}+d\theta_{u_{0}}|_{x}^{2}+|\delta\theta_{u}+\delta\theta_{u_{0}}|_{x}^{2}\\\
u(t)\big{|}_{\partial\Omega}=0.\end{cases}$ (65)

The first term in the above functional can be rewritten as

$\displaystyle\langle\mathcal{L}_{u}\theta,\mathcal{L}_{u}\theta\rangle_{x}+2\langle\mathcal{L}_{u}\theta,\theta_{t}+\mathcal{L}_{u_{0}}\theta\rangle_{x}+|\theta_{t}+\mathcal{L}_{u_{0}}\theta|_{x}^{2}$
(66)

Since $\theta\in H^{1},u_{0}\in H^{1}$, $\mathcal{L}_{u_{0}}\theta\in L^{2}$.
Note that for the space of vector fields vanishing on $\partial\Omega$, the
$L^{2}$ norm of $u$ is bounded by $|d\theta_{u}|^{2}+|\delta\theta_{u}|^{2}$
up to a constant depending on the domain only. For $u,w$ of finite $H^{1}$
norm and vanishing on the boundary, let

$\displaystyle B(u,w)$
$\displaystyle=\int_{\Omega}\langle\mathcal{L}_{u}\theta,\mathcal{L}_{w}\theta\rangle_{x}+\langle
d\theta_{u},d\theta_{w}\rangle_{x}+\langle\delta\theta_{u},\delta\theta_{w}\rangle_{x}dV$
(67) $\displaystyle a(u)$
$\displaystyle=\int_{\Omega}\langle\mathcal{L}_{u}\theta,\theta_{t}+\mathcal{L}_{u_{0}}\theta\rangle_{x}dV$
(68)

Obviously $B$ is symmetric and coercive due to Poincare lemma. We will show
that $B$ is bounded and $a$ is continuous with respect to the $H^{1}$ norm.
Then $B$ gives an equivalent norm as the common $H^{1}$ norm, denoted by
$\|\cdot\|_{B}$. Then by Riesz representation theorem, there exists $f\in
H^{1}(T\Omega)$ which vanishes on $\partial\Omega$, such that $a(u)=\langle
u,f\rangle_{B}$. Then

$\displaystyle B(u,u)+a(u)=\langle u,u\rangle_{B}+2\langle
u,f\rangle_{B}\geq\langle u+f,u+f\rangle_{B}-\|f\|_{B}^{2},$ (69)

implying that $-f$ is the unique solution. Next we show that $B$ and $a$ are
continuous forms.

For the continuity of $a$, following a similar argument as those in appendix
A, there exists finitely many $\eta_{i}\in H^{1}(T^{*}\Omega)$ with bounded
$|d\eta_{i}|_{x}$ and $\xi_{i}\in L^{2}(T^{*}\Omega)$ depending only on the
given data, such that

$\displaystyle
a(u)=\sum_{i}\int_{\Omega}\langle\mathcal{L}_{u}\eta_{i},\xi_{i}\rangle_{x}dV$
(70)

We can further assume that $\eta_{i},\xi_{i}$ are compactly supported inside a
local coordinate $(\Omega_{i},x)$. Since $|d\eta_{i}|_{x}$ is bounded,

$\displaystyle\|\mathcal{L}_{u}\eta_{i}\|_{L^{2}}=\|d\iota_{u}\eta_{i}+\iota_{u}d\eta_{i}\|_{L^{2}}=\|d(\eta_{i}(u))+\iota_{u}d\eta_{i}\|_{L^{2}}\lesssim\|u\|_{H^{1}}.$
(71)

Thus

$\displaystyle|a(u)|\lesssim\sum_{i}\|\mathcal{L}_{u}\eta_{i}\|_{L^{2}}\|\xi_{i}\|_{L^{2}}\lesssim\|u\|_{H^{1}}.$
(72)

In order to show that $B$ is bounded, following the same argument in appendix
A, there exists finitely many differential 1-forms $\eta_{i}$ in $H^{1}$ with
bounded $|d\eta_{i}|_{x}$, and bounded functions $m_{i}$, such that

$\displaystyle\langle\mathcal{L}_{u}\theta,\mathcal{L}_{u}\theta\rangle_{x}=\sum_{i}m_{i}\langle\mathcal{L}_{u}\eta_{i},\mathcal{L}_{u}\eta_{i}\rangle_{x}.$
(73)

Therefore $\|\mathcal{L}_{u}\theta\|^{2}_{L^{2}}\lesssim\|u\|_{H^{1}}^{2}$,
meaning that $B$ is bounded. The proof for the case when $\Omega$ is a compact
oriented Riemannian manifold without boundary is similar, thus omitted.

##  Appendix C Details of the numerical experiments

###  C.1 Model and domain

The data assimilation experiment is conducted using the thermal shallow water
equation [27]. This model consists of three state variables: $h-$the water
height, $v=(v^{1},v^{2})-$the velocity field, and $\Theta-$the buoyancy (or
density contrast):

$\displaystyle\frac{\partial h}{\partial t}+\nabla\cdot(hv)$
$\displaystyle=0,$ (74) $\displaystyle\frac{\partial\Theta}{\partial
t}+(v\cdot\nabla)\Theta$ $\displaystyle=-\kappa(h\Theta-h_{0}\Theta_{0}),$
(75) $\displaystyle\frac{\partial v}{\partial
t}+(v\cdot\nabla)v+f\hat{z}\times v$
$\displaystyle=-\nabla(h\Theta)+\frac{1}{2}h\nabla\Theta.$ (76)

Both $\Theta$ and $h$ are assumed strictly positive at each point.

The quantity $h$ follows a dynamics similar to the state variable $S$ in
example 2.2.2, and transportation terms for $\Theta$ follow those in example
2.2.1. A natural choice is then $\theta_{h}=hdx^{1}\wedge dx^{2}$, a
differential 2-form, and $\theta_{\Theta}=\Theta$ a differential 0-form. The
test case is a double-vortex case. Hence $\theta_{v}=v^{1}dx^{1}+v^{2}dx^{2}$
is chosen to ensure vorticity conservation during the morphing process. These
choices of differential forms differ from those presented in [30]. Less
constrained by the underlying dynamics, [30] discussed the derivation of a
perturbation scheme able to conserve particular quantities. Here, the choice
for the tensor fields should obey to the prescribed dynamics of the system to
maximally maintain the dynamical balance during the morphing process.

The data assimilation process is solely conducted for one time step. It is
assumed that both the vorticity field $\omega=\frac{\partial v^{2}}{\partial
x^{1}}-\frac{\partial v^{1}}{\partial x^{2}}$ and the $h$ field are fully
observed. Since $\theta_{v}=v^{1}dx^{1}+v^{2}dx^{2}$, naturally $\omega$ is
associated to a differential 2-form $\theta_{\omega}=\omega dx^{1}\wedge
dx^{2}=d\theta_{v}$.

The domain is 2-dimensional doubly periodic:
$\Omega=[0,5000\text{km}]_{\text{per}}\times[0,5000\text{km}]_{\text{per}}$,
which is a compact Riemannian manifold without boundary. In this case, a
numerical solution of (14), or equivalently the numerical solution of Eq.(62),
with $a=1$, can be derived in the Fourier space. In fact, the vector field $u$
is separately calculated for $h$ observations and for $\omega$ observations.
From the two observables, the explicit expressions of $u$ are:

$\displaystyle u_{\omega}$
$\displaystyle=(u_{\omega}^{1},u_{\omega}^{2})=2(I-\Delta)^{-1}\big{[}\omega_{2}\nabla(\omega_{1}-\omega_{2})\big{]},$
(77) $\displaystyle u_{h}$
$\displaystyle=(u_{h}^{1},u_{h}^{2})=2(I-\Delta)^{-1}\big{[}h_{2}\nabla(h_{1}-h_{2})\big{]}.$
(78)

The final $u$ for each iterative step is then chosen to be
$\frac{1}{2}(\frac{u_{\omega}}{\|u_{\omega}\|_{1}}+\frac{u_{h}}{\|u_{h}\|_{1}})$,
in which

$\displaystyle\|u\|_{1}^{2}=\int_{\Omega}\langle\theta_{u},\theta_{u}\rangle_{x}+\langle
d\theta_{u},d\theta_{u}\rangle_{x}+\langle\delta\theta_{u},\delta\theta_{u}\rangle_{x}dV=\displaystyle\int_{\Omega}|u|^{2}+|\nabla
u^{1}|^{2}+|\nabla u^{2}|^{2}dx^{1}dx^{2}.$ (79)

###  C.2 Numerical methods and experimental parameters

The units in the code are set to be km and 100s. The initial condition and
model parameters are taken from the numerical experiment in subsection (5.3)
of [10]. The 3-step Adams-Bashforth method (see for instance chapter 3.1 of
[11] ) is used in model integration. Additionally, the one dimensional Hou-Li
spectral filter [12] $\exp\\{-36[(k_{x}/k_{max})^{a}+(k_{y}/k_{max})^{a}]\\}$,
with $a=12$, is applied to the Fourier modes of the state vector at the end of
each model integration step. The truth is generated by running the model
forward for 2750 time units. To generate the ensemble members, the center of
the initial vortex, $(ox,oy)$, is perturbed:

$\displaystyle ox\sim\mathcal{N}(0.1,0.01),\hskip
5.69054ptoy\sim\mathcal{N}(0.1,0.01).$ (80)

The ensemble members are then generated by running the model forward starting
from perturbed initial condition for 2000 time units. For the explicit meaning
of $ox$ and $oy$, please refer to section 5.3 of [10]. The ensemble size
$N_{e}=20$. Both the ensemble members and the truth are generated using a
$256\times 256$ grid. But before starting the morphing process or data
assimilation, the ensemble members and the observations are both projected to
a coarse-grid (64 $\times$ 64) and then interpolated back to the original grid
($256\times 256$).

The morphing process is implemented with the 5-step Adams-Bashforth method and
the spectral Hou-Li filter with $a=36$. For clarity, the pseudo-code of the
complete morphing process is shown in Algorithm 2. We choose
$\epsilon=0.000033$ and $N=10000$.

Data: $\omega^{\text{true}}$, $h^{\text{true}}$, $h^{\text{model}}$,
$v^{\text{model}}$, and $\Theta^{\text{model}}$ on a $256\times 256$ grid of
$\Omega$; small positive number $\epsilon>0$; the total number of iterations
$N$.

Result: morphed model estimates $h^{\text{morphed}},\Theta^{\text{morphed}}$,
and $v^{\text{morphed}}$.

$\omega_{1}\leftarrow\omega^{\text{true}}$; $h_{1}\leftarrow h^{\text{true}}$;
$h_{2}\leftarrow h^{\text{model}}$;
$\Theta_{2}\leftarrow\Theta^{\text{model}}$; $v_{2}\leftarrow
v^{\text{model}}$; $\omega_{2}\leftarrow\displaystyle\frac{\partial
v^{2}_{2}}{\partial x^{1}}-\frac{\partial v^{1}_{2}}{\partial x^{2}}$ ;

Set
$\mathtt{d}h,\mathtt{d}h_{-1},\mathtt{d}h_{-2},\mathtt{d}h_{-3},\mathtt{d}h_{-4},\mathtt{d}\Theta,\mathtt{d}\Theta_{-1},\mathtt{d}\Theta_{-2},\mathtt{d}\Theta_{-3},\mathtt{d}\Theta_{-4},\mathtt{d}v,\mathtt{d}v_{-1},\mathtt{d}v_{-2},\mathtt{d}v_{-3},\mathtt{d}v_{-4}$
all as 0 arrays of the proper sizes ;

for  _$i=1,2,...,N$_ do

Calculate $u_{\omega}$ using Eq.(77);

Calculate $u_{h}$ using Eq.(78);

$u_{i}\leftarrow\frac{1}{2}(\frac{u_{\omega}}{\|u_{\omega}\|_{1}}+\frac{u_{h}}{\|u_{h}\|_{1}})$;

$\mathtt{d}h_{-4}\leftarrow\mathtt{d}h_{-3}$,
$\mathtt{d}h_{-3}\leftarrow\mathtt{d}h_{-2}$,
$\mathtt{d}h_{-2}\leftarrow\mathtt{d}h_{-1}$,
$\mathtt{d}h_{-1}\leftarrow\mathtt{d}h$, do the same for $\mathtt{d}\Theta$
and $\mathtt{d}v$;

$\mathtt{d}h\leftarrow\nabla\cdot(h_{2}u_{i})$, $\mathtt{d}\Theta\leftarrow
u_{i}\cdot\nabla\Theta_{2}$,
$\mathtt{d}v=\displaystyle(v_{2}\cdot\frac{\partial u_{i}}{\partial
x^{2}}+u_{i}\cdot\nabla v_{2}^{1},v_{2}\cdot\frac{\partial u_{i}}{\partial
x^{2}}+u_{i}\cdot\nabla v_{2}^{2})$;

$h_{2}\leftarrow
h_{2}+\epsilon(\frac{1901}{720}\mathtt{d}h-\frac{2774}{720}\mathtt{d}h_{-1}+\frac{2616}{720}\mathtt{d}h_{-2}-\frac{1274}{720}\mathtt{d}h_{-3}+\frac{251}{720}\mathtt{d}h_{-4})$,
do the same for $\Theta_{2}$ and $v_{2}$;

Apply the spectral Hou-Li filter to $h_{2},\Theta_{2}$, and $v_{2}$;

$\omega_{2}\leftarrow\displaystyle\frac{\partial v^{2}_{2}}{\partial
x^{1}}-\frac{\partial v^{1}_{2}}{\partial x^{2}}$;

end for

$h^{\text{morphed}}\leftarrow h_{2}$,
$\Theta^{\text{morphed}}\leftarrow\Theta_{2}$, $v^{\text{morphed}}\leftarrow
v_{2}$

Algorithm 2 The pseudo-code for the nudging process with Adams-Bashforth
method for thermal shallow water model on a doubly periodic domain.

The observation is exactly the same as the interpolated value of the truth.
There is no error in the observation. However, the matrix $R$ in ensemble
Kalman filter is set to be a diagonal matrix with diagonal elements equal to

$\displaystyle
R_{\omega}=\frac{0.01}{64^{2}}\sum_{i=1}^{64^{2}}(\omega^{\text{obs}}_{i})^{2},\text{\hskip
2.84526pt or \hskip
2.84526pt}R_{h}=\frac{0.01}{64^{2}}\sum_{i=1}^{64^{2}}(h^{\text{obs}}_{i})^{2},$
(81)

where $\omega_{i}^{\text{obs}}$ and $h_{i}^{\text{obs}}$ refer to the $\omega$
value and $h$ value at the $i-$th grid-point in the $64\times 64$ grid. Data
assimilation is conducted on the 64$\times$64 grid.

## References

  * [1] Determining optical flow.  Artificial Interlligence, 17:185–203.
  * [2] J. D. Beezley and J. Mandel.  Morphing ensemble kalman filters.  Tellus A: Dyn. Meteorol. Oceanogr., 60:131 – 140, 2007.
  * [3] M. F. Beg, M. I. Miller, A. Trouvé, and L. Younes.  Computing large deformation metric mappings via geodesic flows of diffeomorphisms.  Int. J. Comput. Vision, 61:139–157, 2005.
  * [4] M. Bruveris, F. Gay-Balmaz, D. D. Holm, and T. S. Ratiu.  The momentum map representation of images.  J. Nonlinear Sci., 21:115–150, 2009.
  * [5] T. Corpetti, É. Mémin, and P. Pérez.  Dense estimation of fluid flows.  IEEE Trans. Pattern Anal. Mach. Intell., 24:365–380, 2002.
  * [6] P. Courtier, J.-N. Thépaut, and A. Hollingsworth.  A strategy for operational implementation of 4d-var, using an incremental approach.  Q. J. R. Meteorolog. Soc., 120(519):1367–1387, 1994.
  * [7] P. Del Moral.  Nonlinear filtering: Interacting particle resolution.  Comptes Rendus de l’Académie des Sciences - Series I - Mathematics, 325(6):653–658, 1997.
  * [8] P. Dérian, P. Héas, C. Herzet, and É. Mémin.  Wavelets and optical flow motion estimation.  Numer. Math. theory Methods Appl., 6:116–137, 2013.
  * [9] G. F. D. Duff and D. C. Spencer.  Harmonic tensors on riemannian manifolds with boundary.  Ann. Math., 56:128, 1952.
  * [10] C. Eldred, T. Dubos, and E. Kritsikis.  A quasi-hamiltonian discretization of the thermal shallow water equations.  J. Comput. Phys., 379:1–31, 2019.
  * [11] E. Hairer, G. Wanner, and S. P. Nørsett.  Solving Ordinary Differential Equations I.  Springer series in computational mathematics. Springer Berlin, Heidelberg, 1993.
  * [12] T. Y. Hou and R. Li.  Computing nearly singular solutions using pseudo-spectral methods.  J. Comput. Phys., 226:379–397, 2007.
  * [13] P. Houtekamer and H. L. Mitchell.  A sequential ensemble kalman filter for atmospheric data assimilation.  Mon. Wather Rev., 129:123–137, 2001.
  * [14] P. L. Houtekamer and H. L. Mitchell.  Data assimilation using an ensemble Kalman filter technique.  Mon. Wather Rev., 126(3):796–811, 1998.
  * [15] R. E. Kalman.  A New Approach to Linear Filtering and Prediction Problems.  Journal of Basic Engineering, 82(1):35–45, 03 1960.
  * [16] J. Lefévre and S. Baillet.  Optical flow and advection on 2-riemannian manifolds: A common framework.  IEEE Trans. Pattern Anal. Mach. Intell., 30:1081–1092, 2008.
  * [17] A. C. Lorenc.  Analysis methods for numerical weather prediction.  Quarterly Journal of the Royal Meteorological Society, 112:1177–1194, 1986.
  * [18] E. Meinhardt, J. S. Pérez, and D. Kondermann.  Horn-schunck optical flow with a multi-scale strategy.  Image Process. Line, 3:151–172, 2013.
  * [19] E. Ott and B. R. Hunt.  A local ensemble Kalman filter for atmospheric data assimilation.  Tellus, 56A:415–428, 2004.
  * [20] S. C. Ravela, K. A. Emanuel, and D. B. McLaughlin.  Data assimilation by field alignment.  Physica D, 230:127–145, 2007.
  * [21] S. Rosenberg.  The Laplacian on a Riemannian Manifold: An Introduction to Analysis on Manifolds.  London Mathematical Society Student Texts. Cambridge University Press, 1997.
  * [22] G. Schwarz.  Hodge Decomposition - A Method for Solving Boundary Value Problems.  Springer Berlin, Heidelberg, 1995.
  * [23] A. Trouvé.  Diffeomorphisms groups and pattern matching in image analysis.  International Journal of Computer Vision, 28:213–221, 1998.
  * [24] A. Trouvé and L. Younes.  Local geometry of deformable templates.  SIAM J. Math. Anal., 37:17–59, 2005.
  * [25] A. Trouvé and L. Younes.  Metamorphoses through lie group action.  Found. Comput. Math., 5:173–198, 2005.
  * [26] A. Trouvé.  An infinite dimensional group approach for physics based models in pattern recognition.  International Journal of Computer Vision - IJCV, 01 1995.
  * [27] E. S. Warneford and P. J. Dellar.  The quasi-geostrophic theory of the thermal shallow water equations.  J. Fluid Mech., 723:374 – 403, 2013.
  * [28] R. P. Wildes, M. J. Amabile, A.-M. Lanzillotto, and T.-S. Leu.  Recovering estimates of fluid flow from image sequence data.  Comput. Vis. Image Underst., 80:246–266, 2000.
  * [29] L. Younes, I. Herlin, and D. Béréziat.  A generalized optical flow constraint and its physical interpretation.  In Proceedings IEEE Conference on Computer Vision and Pattern Recognition. CVPR 2000, volume 3, page 2487, Los Alamitos, CA, USA, jun 2000\. IEEE Computer Society.
  * [30] Y. Zhen, V. Resseguier, and B. Chapron.  Physically constrained covariance inflation from location uncertainty.  Nonlinear Process. Geophys., 2023.