Sigma coordinates ?[👆]

[jmbeckers]

Following the discussion in https://github.com/gher-uliege/DIVAnd.jl/issues/123#issuecomment-1410619102 where the coarse vertical resolution with large steps in the topography was mentioned, I wondered if we should implement sigma-coordinates?

Neglecting geometric corrections in the operators (because the grid would not be exactly orthogonal), one could do the following: Take the topography at the resolution of the planned analysis. Generate an interpolating function bfun(x,y) using this gridded topography at the resolution of the subsequent analysis. Define the sigma function ztosigma(x,y,z,bfun) and its inverse function sigmatoz(x,y,s,bfun) Create the data input by applying ztosigma(x,y,z,bfun) to its z coordinate. Create the analysis grid by applying the surface topology mask to a regular grid in sigma space. Analyse in sigma space. Create zi=sigmatoz.(xi,yi,si,bfun) for plotting purposes (e.g. vertical sections pcolor(xi,zi,fi)).

That would be interesting for processes taking place near the bottom, but as for modelling could lead to problems where you have strong stratifications (which might be taken out by defining a reference field in z space first though).

Maybe just a notebook with an example to start with and make a feasibility test ? It would need no modifications in the DIVAnd source.

[ctroupin]

I like the idea. I am not sure I understand well what bfun(x,y) is or does.

[Alexander-Barth]

slightly related, there is currently also the option to make an analysis at fixed distance from the ocean floor:

https://github.com/gher-uliege/DIVAnd.jl/blob/master/src/diva.jl#L52

But this is not a true sigma coordinate, and it not meant to cover the whole water column. The risk I see is that we would "mix" water at very different depth levels (but at the same sigma level), but maybe we can reduce the correlation length near steep topography to mitigate this.

I think within the DIVAndrun function, we can already handle sigma coordinates (by handling it as almost orthogonal) except for the function localize_separable_grid:

https://github.com/gher-uliege/DIVAnd.jl/blob/master/src/localize_separable_grid.jl#L57

But a call to this function can be by-passed by providing the fractional indices (option fracindex of DIVAndrun) corresponding to the location of the observation (directly related to your ztosigma function).

[Alexander-Barth]

I guess bfun(x,y) is the bottom depth at the location x, y....

[jmbeckers]

I like the idea. I am not sure I understand well what bfun(x,y) is or does.

It simply calculates bottom depth at any location (x,y) since one needs that to calculate the sigma position of all data points (whose x,y coordinates will not fall exactly on a grid point). Using Interpolations one could make a simple LinearInterpolation

[jmbeckers]

slightly related, there is currently also the option to make an analysis at fixed distance from the ocean floor:

https://github.com/gher-uliege/DIVAnd.jl/blob/master/src/diva.jl#L52

But this is not a true sigma coordinate, and it not meant to cover the whole water column. The risk I see is that we would "mix" water at very different depth levels (but at the same sigma level), but maybe we can reduce the correlation length near steep topography to mitigate this.

I think within the DIVAndrun function, we can already handle sigma coordinates (by handling it as almost orthogonal) except for the function localize_separable_grid:

https://github.com/gher-uliege/DIVAnd.jl/blob/master/src/localize_separable_grid.jl#L57

But a call to this function can be by-passed by providing the fractional indices (option fracindex of DIVAndrun) corresponding to the location of the observation (directly related to your ztosigma function).

Yes, one could do it directly within DIVAnd but maybe it is easier to transform the coordinates than to prepare the fractional indexes? What would probably change is the interpretation of the vertical length scale. If done in sigma space before calling DIVAnd, one would spread over a fraction of the water column (unless a horizontally variable Lv is used), whereas for the curvilinear grid approach within DIVAnd, the vertical correlation length remains in the real space (more natural).

If one would keep the "within DIVAnd" approach, one would need to create the 3D grid (actually easy by using the sigmatoz function) and calculate the fractional indexes (maybe by using https://github.com/gher-uliege/DIVAnd.jl/blob/master/src/localize_separable_grid.jl#L57 but in sigma space!)

[jmbeckers]

Whatever version, one most certainly would need an Lv which changes spatially to take into account that for the same sigma level, one can be in a very different vertical position.

[jmbeckers]

Also note that once the analysis done (whatever version), if you want to plot really horizontal plots, once could use a fake set of data (x,y) on the regular grid and constant z, calculate the fractional_indices (by transforming into separable sigma space) and then use sparse_interp(mask, I, iscyclic) to get the operator H which allow to interpolate from the 3D grid in sigma space into the desired z-level (the pseudo data). Hope it makes sense.

[jmbeckers]

Proof of concept:

https://github.com/gher-uliege/Diva-Workshops/blob/master/notebooks/3-Analysis/nn-Sigmalayers.ipynb

Both approaches demand some work outside of DIVAndrun and probably it is better to work in the real domain and variable pm. The price to pay (fractional index calculation) is not too high and working in the real domain might be better for other parts where integrals and so on are calculated. One then does not need to bother about the coordinate change. Otherwise both methods are normally identical (in the sense that L*pm remains constant when switching from one representation to the other and the localisation in sigma space is exactly the one used to calculate the fractional indexes for the real space version).

[jmbeckers]

The risk I see is that we would "mix" water at very different depth levels (but at the same sigma level), but maybe we can reduce the correlation length near steep topography to mitigate this.

Or we can use a modified sigma approach:

https://getm.eu/files/GETM/doc/html/node14.html

and keep some layers roughly horizontal. It would not be much more complicated than the classical sigma approach.

[jmbeckers]

One could even use a fake topography (deeper and smoother) to create the vertical grid and once created apply the mask based on real topography, to have a mixture between sigma and z behaviour

[jmbeckers]

Very general vertical coordinate change proof of concept: https://github.com/gher-uliege/Diva-Workshops/blob/master/notebooks/3-Analysis/nn-GeneralizedVerticalCoordinate.ipynb