Mesoscale Ocean Dynamics Modelling
All ocean circulation models, from highly idealized to very comprehensive, operate with viscosity and diffusivity coefficients,
which are by many orders of magnitude larger than those for the sea water. It is hypothesized that artificial increase of these 
parameters accounts for effects of oceanic turbulence phenomena that are not properly captured by the circulation models. 
Some turbulence phenomena, such as breaking internal waves, are so poorly known that we really don't have a better option
for dealing with them. There is a hope, supported by the last 30 years of vigorous research, that we can do much better
for modelling mesoscale oceanic eddies, which can be viewed as “internal weather” of the ocean. By reduction of
the viscosity and diffusivity coefficients, and by simultaneous refinement of the spatial resolution, numerical models
become capable of resolving the mesoscale eddies in a dynamically consistent way. Because of the spatial resolution
requirement, the eddy-resolving ocean circulation models are computationally very expensive. This precludes not only
 proper analysis of these models but also using them for predicting climate variability.
Satellite image of the Gulf Stream ocean current (in red) flowing from west to east across the Atlantic ocean Satellite image of the Gulf Stream current (from Google)
There are 3 ways to move forward with eddy-resolving models: (1) by enhancing computational resources;                         
(2) by improving simple models of the eddy effects; and (3) by improving numerical algorithms. 
The third way, which is in focus of this proposal, requires significant improvement of the existing advection schemes,
which are arguably the main stumbling block. Over the last 10 years I with collaborators have been designing
a new generation of efficient numerical schemes for convection-dominated flows, i.e. for computations of  non-linear flow transport
and propagation at low-viscosity flow regimes (e.g., Karabasov and Goloviznin 2007).The new generation of computational schemes
is needed to replace classical finite-difference schemes, which have large phase and/or amplitude errors,and their high-order extensions, 
which have high computational cost and implementation problems. The new method, which is called “CABARET”, is non-dissipative,
with very low phase error, and based on a very compact computational stencil. The CABARET method for for mesoscale ocean dynamics
is based on the extension of the original second-order Upwind Leapfrog three-time-level advection scheme (e.g., Iserles, 1986; Thomas and Roe, 1993)
to the family of two-time-level non-oscillatory schemes for quasilinear hyperbolic conservation laws. Introduction of the dual variables
staggered in space and time allows CABARET to emphasize several properties of the governing equations, e.g. conservation, low-dispersion and
low-dissipation without numerical oscillations, at once, with a modest computational cost.
 
At present the CABARET has been implemented into an eddy-resolving quasigeostrophic (QG)  model, which simulates the large scale dynamics
of west-boundary eastward jet extension. This is a classical model of wind-driven double-gyre circulation in a midlatitude closed basin of 3840 x 3840 km.
The stratification is represented by stacked layers of constant salinity each (isopycnal), which are
dynamically coupled at the interfaces.The governing equations, which constitute a system of quasi-linear conservation laws with source terms
on the right-hand side due to the effect of viscous forces, the Coriolis force and the wind, are cast in the vorticity-stream function variables. 
At the boundaries a 'partial-slip' boundary condition is applied based on a velocity wall-function with a prescribed boundary layer thickness. 
No explicit subgrid-scale model is used. More information about this type of models, which are the working horse of the modern studies of the mesoscale ocean 
dynamics, can be found at the homepage of my collaborator  Dr.Pavel Berloff in Woods Hole Oceanographic Institution.

Examples of flow animation at the same eddy viscosity regime for the top isopycnal layer of the model (the same 30 contours are plotted in each case) 
The movies show the effect of hydrodynamic instability cascade, from bottom (small) to top (large) scales of the boundary layer, which occurs at
the left computational domain boundary (western boundary), leading to its separation and development into a free jet in the eastward direction.
The physical size of the computational box approximately corresponds to the satellite image of the Gulf Stream ocean current shown above. The lateral eddy viscosity parameter is 12 m^2/s.
Simulation over 100 days with 1 day/frame (3-layer QG model)
CABARET scheme: Potential Vorticity (grid 257x257) 
CABARET scheme: Stream function (grid 257x257) 
Simulation over 1000 days with 10 days/frame (3-layer QG model)
Potential Vorticity anomaly (grid 257x257): left - conventional Arakawa method, right - CABARET 
Stream function (grid 257x257): left - conventional Arakawa method, right - CABARET 
Conventional Arakawa method: Potential Vorticity anomaly (grid 1025x1025)
Conventional Arakawa method: Stream function (grid 1025x1025)
Examples of flow animation in the upper ocean with the conventional 2-nd order scheme at larger eddy viscosity parameters on 257x257 grid, 1000 days with 10 days/frame (3-layer QG model) 
Without refining the grid, the parameterization leads to only a partial improvement of the Eastern Jet solution. 
Conventional Arakawa method: Potential Vorticity anomaly,  viscosity parameter is 30 m^2/s
Conventional Arakawa method: Potential Vorticity anomaly,  viscosity parameter is 50 m^2/s
Conventional Arakawa method: Potential Vorticity anomaly,  viscosity parameter is 100 m^2/s
Conventional Arakawa method: Potential Vorticity anomaly,  viscosity parameter is 200 m^2/s
Conventional Arakawa method: Potential Vorticity anomaly,  viscosity parameter is 800 m^2/s

10-layer QG model: Examples of flow animation in the upper ocean at the eddy viscosity parameter of 100 m^2/s over 1000 days with 10 days/frame  
CABARET scheme: Potential Vorticity anomaly,  257x257 grid
CABARET scheme: Potential Vorticity anomaly,  513x513 grid
Conventional Arakawa method: Potential Vorticity anomaly,  257x257 grid
Without refining the grid, even a heavy parameterization in the case of QG model with 10 vertical layers doesn't lead to any improvement of the EJ solution of the conventional scheme. 
This suggests that the high-resolution of numerical advection scheme is a must if the vertical stratification needs to be resolved in the model. 
Below is  the conventional scheme solution at the latteral viscosity parameter of 800 m^2/s, simulation over 1000 days with 10 days/frame:
Conventional Arakawa method: Potential Vorticity anomaly,  257x257 grid