Enriched evolution of global sea surface height via generalized Schrodinger bridge and Fokker-Planck solver

The Schrödinger bridge problem (SB) was first introduced by Schrödinger in 1932 [17] and now becomes a fundamental framework in the physics, mathematics, engineering and information geometry to modeling the ensemble statistical transition path between fixed initial and final distributions.

Given initial density $\rho_{0}$ and final density $\rho_{1}$, denote (SB) as the Schrödinger bridge problem

(2.1)		$\displaystyle\min_{b,\rho}\frac{1}{2}\int\int|b|^{2}\rho\,\mathrm{d}x\,\mathrm{d}t$
	$\displaystyle\operatorname{s.t.~{}}\partial_{t}\rho+\nabla\cdot(\rho b)=\gamma\Delta\rho,\quad\rho_{0}=\mu,\,\,\rho_{1}=\nu.$

When $\gamma\to 0$, the optimal $b$ is the optimal velocity field in Optimal transport.

It worth to notice the above Schrödinger Bridge problem can be generalized to include an energy functional $E$ beyond the entropy. An internal energy showing an underlying driving potential can be used to guide the evolution of SSH patterns, which may be affected by multiple factors as a result of global warming. The combined driven force of those factors is usually vague and difficult to distinguish from each other. Thus the art of choosing proper underlying driving potential is case by case.

For simplicity, we consider distributions in the probability space $\mathcal{P}$ on $\mathbb{T}^{d}$ or $\mathbb{R}^{d}$. Denote the final distribution of SSH as $\pi$. Based on the Boltzmann analysis, we can regard the final distribution as a probability that is proportional to the exponential of a potential difference

(2.2)

$$\pi\propto e^{-\frac{1}{\gamma}U}.$$

Then we use the relative entropy as a typical internal energy

(2.3)

$$E=\gamma\int\rho\ln\frac{\rho}{\pi}\,\mathrm{d}x,\quad\pi\propto e^{-\frac{1}{\gamma}U}.$$

It is easy to calculate the first variation of $E$ is

(2.4)

$$\frac{\delta E}{\delta\rho}=\gamma\ln\frac{\rho}{\pi}+\gamma=\gamma\ln\rho+U+\gamma.$$

One have the general Schrödinger bridge problem (SBg)

(2.5)			$\displaystyle\min_{b,\rho}\frac{1}{2}\int_{0}^{1}\int|b|^{2}\rho\,\mathrm{d}x\,\mathrm{d}t$
		$\displaystyle\operatorname{s.t.~{}}\partial_{t}\rho+\nabla\cdot(\rho b)=\nabla\cdot\left(\rho\nabla\frac{\delta E}{\delta\rho}\right),\quad\rho_{0}=\mu,\,\,\rho_{1}=\nu.$

It is equivalent to the so called general Yasue problem [19]

(2.6)		$\displaystyle\min\int_{0}^{1}\int\frac{1}{2}|v|^{2}\rho\,\mathrm{d}x\,\mathrm{d}t+\frac{1}{2}|\nabla\frac{\delta E}{\delta\rho}|^{2}\rho\,\mathrm{d}x\,\mathrm{d}t$
	$\displaystyle\operatorname{s.t.~{}}\partial_{t}\rho+\nabla\cdot(\rho v)=0,\quad\rho_{0}=\mu,\,\,\rho_{1}=\nu.$

Indeed, this can be seen by the changing variable $b=v+\nabla\frac{\delta E}{\delta\rho}$ in (2.5). Then one can directly compute the running cost with Lagrangian $L$,

	$\displaystyle L:=$	$\displaystyle\frac{1}{2}\int|b|^{2}\rho\,\mathrm{d}x$
	$\displaystyle=$	$\displaystyle\int\left(\frac{1}{2}|v|^{2}\rho+\frac{1}{2}|\nabla\frac{\delta E}{\delta\rho}|^{2}\rho+\rho v\cdot\nabla\frac{\delta E}{\delta\rho}\right)\,\mathrm{d}x$
	$\displaystyle=$	$\displaystyle\int\left(\frac{1}{2}|v|^{2}\rho+\frac{1}{2}|\nabla\frac{\delta E}{\delta\rho}|^{2}\rho+\partial_{t}\rho\frac{\delta E}{\delta\rho}\right)\,\mathrm{d}x,$

where we used the constraint $\partial_{t}\rho+\nabla\cdot(\rho v)=0$. Notice the last term above is $\frac{\,\mathrm{d}E}{\,\mathrm{d}t}$ which after the time integration is a constant depending only on the fixed initial distribution and the ending distribution.

We remark that there is another way to recast the corresponding Lagrangian function $L$ as

(2.7)		$\displaystyle L:=$	$\displaystyle\frac{1}{2}\int|b|^{2}\rho\,\mathrm{d}x=\frac{1}{4}\|\nabla\cdot(\rho b)\|_{H^{-1}(\rho)}^{2}$
(2.8)		$\displaystyle=$	$\displaystyle\frac{1}{4}\|\partial_{t}\rho-\nabla\cdot\left(\rho\nabla\frac{\delta E}{\delta\rho}\right)\|^{2}_{H^{-1}(\rho)}.$

This $L$ indeed measures the residual of gradient flow of $E$ in the $H^{-1}(\rho)$ norm. This $L$ is also the $L$-function in the Large derivation principle of the average process $X_{\varepsilon}=\frac{1}{N}X_{i}=\gamma X_{i}$, with $X_{i}$ being the i.i.d Markov process with generator $Qf=\Delta f-\nabla U\cdot\nabla f$ [5, 9, 12].

In the next section, we will use the generalized Schrödinger Bridge problem with various energy landscape to design algorithms in image morphing and transition state theory.

In this section, we convert the optimization problem (2.5) as a coupled PDE system with two boundary conditions.

First we derive the Euler-Langrange-Bellman equation for the generalized Schrödinger bridge problem (2.5).

One can derive the Euler-Lagrange equations using Lagrangian multiplier $\phi_{t}(x)$

(2.9)

$$\sup_{b_{t},\rho_{t}}\inf_{\Phi_{s}}\left(\int_{0}^{1}\left(\int\left(\frac{1}{2}|b|^{2}\rho_{t}(x)+\phi_{t}(x)[\partial_{t}\rho_{t}+\nabla\cdot(\rho_{t}b_{t})-\nabla\cdot\left(\rho\nabla\frac{\delta E}{\delta\rho}\right)]\right)\,\mathrm{d}x\right)\,\mathrm{d}s\right).$$

Notice $\nabla\frac{\delta E}{\delta\rho}=\gamma\frac{\nabla\rho}{\rho}+\nabla U$. Then the constraint becomes the Fokker-Planck equation for the drift-diffusion process

(2.10)

$$\partial_{t}\rho_{t}+\nabla\cdot(\rho_{t}b_{t})=\nabla\cdot\left(\gamma\nabla\rho+\rho\nabla U\right).$$

Thus the Euler-Lagrange equations are

(2.11)		$\displaystyle\partial_{t}\rho_{t}+\nabla\cdot(\rho_{t}b_{t})=\nabla\cdot\left(\gamma\nabla\rho+\rho\nabla U\right),$
		$\displaystyle\frac{1}{2}|b|^{2}-\partial_{t}\phi_{t}-b\cdot\nabla\phi+\nabla\phi\cdot\nabla U-\gamma\Delta\phi=0,$
		$\displaystyle b(x)-\nabla\phi_{t}=0.$

Therefore, the optimal control velocity is given by

$$b_{t}(x)=\nabla\phi_{t}(x)$$

and the Euler-Lagrange-Bellman equation becomes

(2.12)		$\displaystyle\partial_{t}\phi+\frac{1}{2}|\nabla\phi|^{2}-\nabla\phi\cdot\nabla U+\gamma\Delta\phi=0,$
		$\displaystyle\partial_{t}\rho_{t}+\nabla\cdot\left(\rho_{t}(\nabla\phi_{t}-\nabla U)\right)=\gamma\Delta\rho,$
		$\displaystyle\rho_{t=0}=\mu,\quad\rho_{t=1}=\nu.$

Now we introduce a Hamiltonian functional $\mathcal{H}:\mathcal{P}(\mathbb{T}^{d})\times C^{1}(\mathbb{T}^{d})\to\mathbb{R}$

(2.13)

$$\mathcal{H}(\rho(\cdot),\phi(\cdot)):=\int_{\mathbb{T}^{d}}\left(\frac{1}{2}|\nabla\phi|^{2}-\langle\nabla\phi,\nabla\frac{\delta E}{\delta\rho}\rangle\right)\rho(x)\,\mathrm{d}x.$$

Then by elementary computations, we have

	$\displaystyle\frac{\delta\mathcal{H}}{\delta\phi}(\rho,\phi)=-\nabla\cdot\left(\rho\nabla\phi\right)+\nabla\cdot\left(\rho\nabla\frac{\delta E}{\delta\rho}\right),$
	$\displaystyle\frac{\delta\mathcal{H}}{\delta\rho}(\rho,\phi)=\frac{1}{2}|\nabla\phi|^{2}+\gamma\Delta\phi-\nabla\phi\cdot\nabla U.$

Thus the Euler-Lagrange-Bellman equation (2.12) becomes a Hamiltonian system in infinite dimensional space

(2.14)		$\displaystyle\partial_{t}\rho=\frac{\delta\mathcal{H}}{\delta\phi}(\rho,\phi),$
		$\displaystyle\partial_{t}\phi=-\frac{\delta\mathcal{H}}{\delta\rho}(\rho,\phi).$

Then by the Hopf-Cole transformation (c.f. [3, 6, 9]): $(\eta^{\text{bw}},\eta^{\text{fw}})\to(\rho,\phi)$ defined as

(2.15)		$\displaystyle\delta E(\rho)=\delta E(\eta^{\text{bw}})+\delta E(\eta^{\text{fw}}),$
(2.16)		$\displaystyle\phi=2\delta E(\eta^{{\text{bw}}}),$

where $\delta E(\rho)=\gamma\log\rho+\gamma+U$ is the first variation of $E$.

Then we have

(2.17)

$\displaystyle\eta^{\text{bw}}=e^{\frac{\phi}{2\gamma}}e^{-\frac{U}{\gamma}-1},\quad\eta^{\text{fw}}=\rho e^{-\frac{\phi}{2\gamma}}.$

Or equivalently

(2.18)		$\displaystyle\phi=2\gamma\ln\eta^{\text{bw}}+2U+2\gamma,$
(2.19)		$\displaystyle\rho=e^{\frac{\phi}{2\gamma}}\eta^{\text{fw}}=\eta^{\text{bw}}\eta^{\text{fw}}e^{\frac{U}{\gamma}+1}.$

With this transformation, one can check

	$\displaystyle\frac{\nabla\phi}{2}=\gamma\frac{\nabla\eta^{\text{bw}}}{\eta^{\text{bw}}}+\nabla U,$
	$\displaystyle\gamma\frac{\nabla\rho}{\rho}-\frac{\nabla\phi}{2}=\gamma\frac{\nabla\eta^{\text{fw}}}{\eta^{\text{fw}}}.$

Then in terms of the $(\eta^{\text{fw}},\eta^{\text{bw}})$ the ELB equations (2.12) are

(2.20)		$\displaystyle\partial_{t}\eta^{\text{bw}}+\gamma\Delta\eta^{\text{bw}}+\nabla\eta^{\text{bw}}\cdot\nabla U+\eta^{\text{bw}}\Delta U=0,$
		$\displaystyle-\partial_{t}\eta^{\text{fw}}+\gamma\Delta\eta^{\text{fw}}+\nabla\eta^{\text{fw}}\cdot\nabla U+\eta^{\text{fw}}\Delta U=0.$

or equivalently

(2.21)		$\displaystyle\partial_{t}\eta^{\text{bw}}+\gamma\Delta\eta^{\text{bw}}+\nabla\cdot(\eta^{\text{bw}}\nabla U)=0,$
		$\displaystyle-\partial_{t}\eta^{\text{fw}}+\gamma\Delta\eta^{\text{fw}}+\nabla\cdot(\eta^{\text{fw}}\nabla U)=0.$

Define the Fokker-Planck operator

(2.22)

$$FP^{U}(\eta):=\gamma\Delta\eta+\nabla\cdot(\eta\nabla U)=\nabla\cdot\left(\eta(\gamma\nabla\ln\eta+\nabla U)\right).$$

With a given equilibrium $\pi=e^{-U}$, denote $m:=\frac{1}{\gamma}$ and $\pi^{m}=e^{-\frac{U}{\gamma}}$. With $\pi^{m}=e^{-\frac{U}{\gamma}}$, one can rewrite the FP operator as

(2.23)		$\displaystyle FP^{\pi^{m}}(\eta):=$	$\displaystyle\gamma\Delta\eta+\nabla\cdot(\eta\nabla U)=\gamma\nabla\cdot\left(\eta\left(\nabla\ln\eta+\frac{1}{\gamma}\nabla U\right)\right)$
		$\displaystyle=$	$\displaystyle\gamma\nabla\cdot\left(\eta\nabla\ln\frac{\eta}{\pi^{m}}\right)=\gamma\nabla\cdot\left(\pi^{m}\nabla\frac{\rho}{\pi^{m}}\right).$

Therefore the $(\eta^{\text{fw}},\eta^{\text{bw}})$ equations become

(2.24)		$\displaystyle-\partial_{t}\eta^{\text{bw}}=FP^{U}(\eta^{\text{bw}}),\quad\eta^{\text{bw}}|_{t=1}=\eta^{\text{bw}}_{1},$
		$\displaystyle\partial_{t}\eta^{\text{fw}}=FP^{U}(\eta^{\text{fw}}),\quad\eta^{\text{fw}}|_{t=0}=\eta^{\text{fw}}_{0}.$

The gSB problem is to find the data $\eta^{\text{fw}}_{0}$ and $\eta^{\text{bw}}_{1}$ such that

(2.25)		$\displaystyle\rho_{0}=\eta^{\text{fw}}_{0}\eta^{\text{bw}}(0)e^{\frac{U}{\gamma}+1},$
(2.26)		$\displaystyle\rho_{1}=\eta^{\text{fw}}(1)\eta^{\text{bw}}_{1}e^{\frac{U}{\gamma}+1}.$

As a result, the obtained $\eta^{\text{fw}}(t)$ and $\eta^{\text{bw}}(t)$ can be used to construct the gSB solution

(2.27)

$$\rho(t)=\eta^{\text{fw}}(t)\eta^{\text{bw}}(t)e^{\frac{U}{\gamma}+1}.$$

The Fokker-Planck equation can be solved via the finite volume method introduced in [16, 4, 14, 8], which can also be generalized to the solver for Fokker-Planck equations on a manifold. The first difficulty here is the initial data and the ending data for $\eta^{{\text{fw}}}$ and for $\eta^{bw}$ are unknown. Thus some iterative shooting method for finding $\eta^{{\text{fw}}}$ and $\eta^{bw}$ based on given $\rho_{0}$ and $\rho_{1}$ are necessary. Some existing ones are [18, 2]. On the other hand, how to design the potential $U$ in the forward and the backward equations is an art based on the application in the current context. We refer to [10, 11, 12] for various design of either a reversible or an irreversible drift to adjust the optimally controlled trajectories in the Fokker-Planck equation.

In this section, we introduce the double control method for computing the optimal transformation path. The idea is illustrated in figure 1.

Given the initial/ending distributions $p_{0}=\mu,\,\,p_{1}=\nu$, one first solves the forward control problem with the driving potential $U_{1}=-\gamma\log\nu.$ Then with the same initial/ending distributions $q_{0}=\mu,\,\,q_{1}=\nu$, one solves the backward control problem with the driving potential $U_{0}=-\gamma\log\mu.$

Precisely, define the forward control problem

	$\displaystyle\min_{u,p}$	$\displaystyle\frac{1}{2}\int\int|u|^{2}p\,\mathrm{d}x\,\mathrm{d}t$
(2.28)		$\displaystyle\operatorname{s.t.~{}}$	$\displaystyle\partial_{t}p+\nabla\cdot(pu)=\nabla\cdot\left(p\nabla\frac{\delta E}{\delta p}\right)=\gamma\Delta p+\nabla\cdot\left(p\nabla U_{1}\right)$
		$\displaystyle p_{0}=\mu,\,\,p_{1}=\nu.$

Define the backward control problem

	$\displaystyle\min_{v,q}$	$\displaystyle\frac{1}{2}\int\int|v|^{2}q\,\mathrm{d}x\,\mathrm{d}t$
(2.29)		$\displaystyle\operatorname{s.t.~{}}$	$\displaystyle-\partial_{t}q+\nabla\cdot(qv)=\nabla\cdot\left(q\nabla\frac{\delta E}{\delta q}\right)=\gamma\Delta q+\nabla\cdot\left(q\nabla U_{0}\right)$
		$\displaystyle q_{0}=\mu,\,\,q_{1}=\nu.$

Then the Euler-Lagrange-Bellman equation to the forward control problem (2.4) is

(2.30)		$\displaystyle\partial_{t}\phi+\frac{1}{2}|\nabla\phi|^{2}+\gamma\Delta\phi-\nabla\phi\cdot\nabla U_{1}=0,$
		$\displaystyle\partial_{t}p+\nabla\cdot\left(p(\nabla\phi-\nabla U_{1})\right)=\gamma\Delta p$

with $u=\nabla\phi$. And the Euler-Lagrange-Bellman equation to the backward control problem (2.4) is

(2.31)		$\displaystyle-\partial_{t}\varphi+\frac{1}{2}|\nabla\varphi|^{2}+\gamma\Delta\varphi-\nabla\varphi\cdot\nabla U_{0}=0,$
		$\displaystyle-\partial_{t}q+\nabla\cdot\left(q(\nabla\varphi-\nabla U_{0})\right)=\gamma\Delta q$

with $v=\nabla\varphi$.

Under the same Hopf-Cole transformation (2.15) for $(p,\phi)$,

(2.32)		$\displaystyle\phi=2\gamma\ln\eta^{\text{bw}}+2U_{1}+2\gamma,$
(2.33)		$\displaystyle p=e^{\frac{\phi}{2\gamma}}\eta^{\text{fw}}=\eta^{\text{bw}}\eta^{\text{fw}}e^{\frac{U_{1}}{\gamma}+1},$

we have the $(\eta^{\text{fw}},\eta^{\text{bw}})$ equations

(2.34)		$\displaystyle-\partial_{t}\eta^{\text{bw}}=FP^{U_{1}}(\eta^{\text{bw}}),\quad\eta^{\text{bw}}(1)=\eta^{\text{bw}}_{1},$
		$\displaystyle\partial_{t}\eta^{\text{fw}}=FP^{U_{1}}(\eta^{\text{fw}}),\quad\eta^{\text{fw}}(0)=\eta^{\text{fw}}_{0}.$

Then forward control problem (2.4) is to find the data $\eta^{\text{fw}}_{0}$ and $\eta^{\text{bw}}_{1}$ such that

(2.35)		$\displaystyle p_{0}=\eta^{\text{bw}}(0)\eta^{\text{fw}}_{0}e^{\frac{U_{1}}{\gamma}+1},$
(2.36)		$\displaystyle p_{1}=\eta^{\text{bw}}_{1}\eta^{\text{fw}}(1)e^{\frac{U_{1}}{\gamma}+1}.$

Particularly, when $p_{1}=e^{-\frac{U_{1}}{\gamma}-1}$, we have $1=\eta_{1}^{\text{bw}}\eta^{\text{fw}}(1).$ Thus the data for $\eta$ equation is

(2.37)		$\displaystyle p_{0}p_{1}=\eta^{\text{bw}}(0)\eta^{\text{fw}}_{0},$
(2.38)		$\displaystyle p_{1}^{2}=\eta^{\text{bw}}_{1}\eta^{\text{fw}}(1).$

As a result, the obtained $\eta^{\text{fw}}(t)$ and $\eta^{\text{bw}}(t)$ can be used to construct the solution to the forward control problem (2.4)

(2.39)

$$p(t)p_{1}=\eta^{\text{bw}}(t)\eta^{\text{fw}}(t).$$

Similarly, under the Hopf-Cole transformation (2.15) for $(q,\varphi)$,

(2.40)		$\displaystyle\varphi=2\gamma\ln\xi^{\text{fw}}+2U_{0}+2\gamma,$
(2.41)		$\displaystyle q=e^{\frac{\varphi}{2\gamma}}\xi^{\text{bw}}=\xi^{\text{bw}}\xi^{\text{fw}}e^{\frac{U_{0}}{\gamma}+1},$

we have the equations for $(\xi^{\text{fw}},\xi^{\text{bw}})$

(2.42)		$\displaystyle\partial_{t}\xi^{\text{fw}}=FP^{U_{0}}(\xi^{\text{fw}}),\quad\xi^{\text{fw}}(0)=\xi^{\text{fw}}_{0},$
		$\displaystyle-\partial_{t}\xi^{\text{bw}}=FP^{U_{0}}(\xi^{\text{bw}}),\quad\xi^{\text{bw}}(1)=\xi^{\text{bw}}_{1}.$

Then backward control problem (2.4) is to find the data $\xi^{\text{fw}}_{0}$ and $\xi^{\text{bw}}_{1}$ such that

(2.43)		$\displaystyle q_{0}=\xi^{\text{bw}}(0)\xi^{\text{fw}}_{0}e^{\frac{U_{0}}{\gamma}+1},$
(2.44)		$\displaystyle q_{1}=\xi^{\text{bw}}_{1}\xi^{\text{fw}}(1)e^{\frac{U_{0}}{\gamma}+1}.$

Particularly, when $q_{0}=e^{-\frac{U_{0}}{\gamma}-1}$, we have $q_{0}^{2}=\xi_{0}^{\text{fw}}\xi^{\text{bw}}(0).$ Thus the data for $\xi$ equation is

(2.45)		$\displaystyle q_{0}^{2}=\xi^{\text{bw}}(0)\xi^{\text{fw}}_{0},$
(2.46)		$\displaystyle q_{1}q_{0}=\xi^{\text{bw}}_{1}\xi^{\text{fw}}(1).$

As a result, the obtained $\xi^{\text{fw}}(t)$ and $\xi^{\text{bw}}(t)$ can be used to construct the solution to the forward control problem (2.4)

(2.47)

$$q(t)q_{0}=\xi^{\text{bw}}(t)\xi^{\text{fw}}(t).$$

We now design an algorithm based on the generalized Sinkhorn algorithms [18]. The convergence of this algorithm is an interesting question which need some analysis. In the example for the evolution of SSH, we numerically show the convergence of this iteration method.

Denote the numerical solution at $k$-th iteration as $\eta^{{\text{fw}},k},\eta^{{\text{bw}},k}$. Below we describe the algorithm.   
Input: Begin image $p_{0}$, End image $p_{1}$, Initial guess $\eta^{{\text{fw}},0}_{0}$.
Step 1. Given data $\eta^{{\text{fw}},k}_{0}$, solve FP (3.1) to time $t=1$. Denote the obtained solution as $\eta^{{\text{fw}},k}(1)$.
Step 2. Based on (3.4), solve $\xi^{{\text{bw}},k}_{1}$ by

(3.6)

$$p_{1}^{2}=\xi^{{\text{bw}},k}_{1}\eta^{{\text{fw}},k}(1).$$

Step 3. Using data $\xi^{{\text{bw}},k}_{1}$, solve FP (3.2) to time $t=0$. Denote the obtained solution as $\xi^{{\text{bw}},k}(0).$
Step 4. Based on (3.3), solve $\eta^{{\text{fw}},k+1}_{0}$ by

(3.7)

$$p_{0}^{2}=\xi^{{\text{bw}},k}(0)\eta^{{\text{fw}},k+1}_{0}.$$

Output: The converged $\eta_{0}^{{\text{fw}}}:=\lim_{k\to\infty}\eta_{0}^{{\text{fw}},k}$, $\xi_{1}^{{\text{bw}}}:=\lim_{k\to\infty}\xi_{1}^{{\text{bw}},k}$. Then with these $\eta_{0}^{\text{fw}},\xi^{{\text{bw}}}_{1}$ solve $\eta^{\text{fw}}(t)$,$\xi^{\text{bw}}(t)$ and construct

(3.8)

$$\rho(t)^{2}=\eta^{\text{fw}}(t)\xi^{\text{bw}}(t),\quad 0\leq t\leq 1.$$

In this paper, the global SSH data are derived from the HYCOM (HYbrid Coordinate Ocean Model) re-analysis dataset. The HYCOM consortium is a multi-institutional effort sponsored by the National Ocean Partnership Program (NOPP), as part of the U. S. Global Ocean Data Assimilation Experiment (GODAE), to develop and evaluate a data-assimilative hybrid isopycnal-sigma-pressure (generalized) coordinate ocean model. The GODAE objectives of three-dimensional depiction of the ocean state at fine resolution in real time, provision of boundary conditions for coastal and regional models, and provision of oceanic boundary conditions for a global coupled ocean-atmosphere prediction model, are being addressed by a partnership of institutions that represent a broad spectrum of the oceanographic community.

HYCOM is a product from the HYCOM Consortium that forecasts global ocean conditions. The HYCOM $+$ NCODA Global $1/12$ degree Analysis (GLBy0.08/expt_93.0) provides five variables including Sea Surface Height as well as eastward velocity, northward velocity, in-situ temperature, and salinity at 40 vertical depth levels ranging from 0 m at the ocean surface to 5,000 m. It combines the model results and satellite observation data with the horizontal grid to be 0.08 degree along longitude and 0.08 degree along latitude.

In the experiments of this paper, the starting point is set to be the monthly averaged SSH in August, 1994 and the end point is set to be the monthly averaged SSH in August, 2014. The goal for the experiment is to simulate the optimal evolution path from the starting point to the end. In order to monitor the simulation results in detail, results in the northern Pacific region will be viewed and analyzed. This region is specially selected because the existence of Kuroshio current. Similar to the Gulf Stream in the North Atlantic, the Kuroshio Extension is a dynamic but relatively unstable system, with variability in the associated bifurcation latitude occurring on interannual time scales. The cause of these variations and their effects on the surface flow and total transport of waters has been studied extensively, SSH can be a direct indicator for Kuroshio so with recent advances in sea surface height satellite altimetry methods allowing for observational studies on larger timescales. As can be seen in figure 2 and figure 3, both SSH patterns clearly indicate the Kuroshio current and its extension.