[🐞] When running DIVAnd, artifacts appear which seem to be linked to the bathymetry

[qcazenave]

No problem with running DIVAnd so I am not sure whether the "bug" type is relevant.

My problem is that when looking at the maps generated by diva3d, we noticed some gradients that seem to exactly correspond to change in the bottom depth as if information from below were transmitted to the layer above but not spread over nearby regions where there is no information available from below since it is already the bottom depth... Is there a way to suppress this effect ?

Below is the way we run diva3d as well as an example of the ouput maps, first without and then, with the isobaths :

`mask = falses(size(b,1),size(b,2),length(depthr)); #b is bathymetry obtained with the "extract_bath" function for k = 1:length(depthr) mask[:,:,k] .= (b .>= depthr[k]) end;

dbinfo = diva3d((lonr,latr,depthr,TS), (obslon,obslat,obsdepth,obstime), obsval, (lenx,leny,lenz2), epsilon2, filename,varname, bathname = bathname, bathisglobal = bathisglobal, mask=mask,

plotres = plotres_timeindex,

ncvarattrib = ncvarattrib,
ncglobalattrib = ncglobalattrib,
timeorigin = timeorigin,
transform = Anam.loglin(slog),
memtofit = 100
)`

[jmbeckers]

Can you try with ;alphabc=0 as additional keyword in diva3d? Normally it is accessible directly from diva3d.

I think the problem is due to vertical windowing (do you see the output proceeding by windows moving upward ?), associated with the alphabc=1 applied as boundary condition now on the "top" and "bottom" layer of the window.

[qcazenave]

I tried with alphabc = 0 but it does not change the result. I don't really understand what this parameters does, nor what youn mean by "do you see the output proceeding by windows moving upward ?" Sorry :/

[jmbeckers]

alphabc changes the way boundary conditions are applied at the edges of the domain. What I mean with the windows moving upward is the output of diva3d which tells you something like

Info: number of windows: 

and if you have more than 1 window it means a domain decomposition is at work.

[jmbeckers]

I tried with alphabc = 0 but it does not change the result.

Not even on the near bottom layers (I mean not necessarily correcting the artefacts but slightly changing values) ?

[qcazenave]

Ok... then it says : Info: number of windows: 1 with alphabc=0 (but it said the same with the default parameterization...)

[qcazenave]

OK my mistake, I went too fast, it is indeed different.

[jmbeckers]

OK my mistake, I went too fast, it is indeed different.

But not eliminating the artefacts ?

[qcazenave]

unfortunately not...

[jmbeckers]

Thanks for the plots (I assume the dots on the vertical indicate the grid position on the vertical). So it is not even a layer closest to the "bottom". Very strange. Would it be possible to share the setup (or better a simplified one, for example same grid but toy data where the problem still exists) ?

[jmbeckers]

Can you also have a look at a vertical section at -3 longitude ?

[jmbeckers]

Are the plots done using plotres in diva3d or are they done afterwards via a netcdf reading (matlab?)

[qcazenave]

Plots are done afterwards with a netcdf reading (pyplot): it is a map (lon,lat) fo the field at 75m deep (depth levels are from 0 to 125m and the last 3 layers are 75, 100 and 125m) but this artefact is visible at other depths as well, above and below (at 40m, you can even see the isobaths 50m and 75m in the field features). The features are also visible on panoply but it is not that obvious (due to the different colormaps, probably...). I will have a look at the cross section.

[jmbeckers]

Plots are done afterwards with a netcdf reading (pyplot): it is a map (lon,lat) fo the field at 75m deep (depth levels are from 0 to 125m and the last 3 layers are 75, 100 and 125m) but this artefact is visible at other depths as well, above and below (at 40m, you can even see the isobaths 50m and 75m in the field features). The features are also visible on panoply but it is not that obvious (due to the different colormaps, probably...). I will have a look at the cross section.

Can you do the plots within diva3d with plotres to exclude the possibility of a problem when creating the netcdf file ? That would make sure the artefact is due to the analysis itself.

[Alexander-Barth]

Can you also make a plot of the background profile dbinfo[:background_profile] (2d array: depth x time) and what is the vertical correlation length?

[qcazenave]

Result of plotres :

[jmbeckers]

Thanks for the effort. So definitively due to the analysis itself.

[qcazenave]

background profiles (corrected from log transform...) :

[jmbeckers]

Log transformation via anamorphosis ? I think the background field is stored without the inverse transform applied (probably we should correct that storing anyway, but I do not think it is the problem here, since the calculations continue with anomalies in the log space.)

Ok, I see you edited your post by applying the inverse transform.

[jmbeckers]

I try to recreate something similar with a toy example. The only situation where I could recreate something into this direction is a situation with a strong vertical gradient, low vertical resolution and large vertical correlation length.

What is lenz2 in your case ?

[qcazenave]

Yeas, sorry I forgot to give you the vertical correlation length -> so I re-used the correlation lengths that had been calculated with fithorzlen et fitvertlen in the 4th phase. As you can see, the vertical correlation length is quite large compared to the vertical profile... so this might explain that... I will try a new run with much smaller values for lenz to see.

Lenz (one column for each season)

Vertical resolution

[jmbeckers]

Well then I think the vertical length scale is definitively too high ! You basically ask to mix the whole water column which then naturally leads to steps above places where the topography changes.

[qcazenave]

So I used a new vertical correlation length (see image below) the maps are different but the effect is still visible.

[qcazenave]

Vertical cross-section:

[jmbeckers]

Vertical sections actually seem quite ok but they also show the lack of vertical resolution of the grid compared to the horizontal resolution. There are large regions of flat topography simply due to lack of vertical resolution. So my approach would use higher vertical resolution (and then possibly even smaller vertical correlation length; with 50m you still mix a lot near the bottom)

[qcazenave]

OK, I thought we had to use correlation length = min(2*dz) which is why for the deepest layers the correlation length is 50m (resolution = 25m from 50m to 125m)

[jmbeckers]

OK, I thought we had to use correlation length = min(2*dz) which is why for the deepest layers the correlation length is 50m (resolution = 25m from 50m to 125m)

Yep, that is also why I suggest to increase resolution :-)

[jmbeckers]

Unless proved differently, I think we can consider that the artefacts are due to coarse vertical resolution, (too) large vertical resolution in a stratified situation. I do not think that we can provide corrections to that other than increasing resolution and I would close the topic for the time being. Comments can still be added of course.

For DIVAnd I suggest to a) adapt the background field information by applying the inverse transformation b) add some sanity checks in vertical correlation length by show a warning if it's value is close to the domain size

[qcazenave]

Hi, coming back to you. So I need to increase the vertical resolution, OK. But then if I want the correlation length (horizontal or vertical) to be determined with fitlen, I get an error because there is not enough data in the region compared to the vertical resolution:

apparently, the fithorzlen function cannot deal with an empty "fitinfos" output from ther fitlen function...

Is there an option to prevent this error from happening ?

[jmbeckers]

Can you run with

ENV["JULIA_DEBUG"] = "DIVAnd".

Do you see an output of fitting which works well near the surface and then fails in deeper parts ?

Also, before fitting L, does the higher vertical resolution with imposed L improve the situation compared to the coarser resolution ?

[qcazenave]

Yes, it only fails for deeper parts (below 75m) Ok I will run with the 2 different resolutions but an imposed L to confirm it works better.

[jmbeckers]

Normally when there is a problem in one layer, there is a vertical filtering of the estimations later anyway. It seems that it however request the :rqual to be set, which is not the case when there is no data at all.

So

we would need to add

dbinfo = Dict{Symbol,Any}(
        :covar => covar,
        :fitcovar => fitcovar,
        :distx => distx,
        :rqual => rqual,
        :range => range,
        :covarweight => covarweight,
        :distx => distx,
        :sn => SN,
        :meandist => meandist,
    )

instead of the empty Dict{Symbol,Any}() when catching

https://github.com/gher-uliege/DIVAnd.jl/blob/df4b25bbf92988f61333e943ebf322a313456c2d/src/fit.jl#L428

And using the dict with appropriate default values when there are not data. Then the vertical filtering should be able to proceed.

[jmbeckers]

Now until the fix is implemented and tested, for your case, I would suggest to only fit down to 70m and once happy with it, copy the values of the correlation length scale found at 70m into the deeper layers (since you have anyway very few data, that seems a reasonable assumption)

[qcazenave]

I tried with a finer resolution (32 depth levels : depth = [0:2:20 ..., 25:5:125 ...] ) and I get: Unhandled Task ERROR: SuiteSparse.CHOLMOD.CHOLMODException("out of memory") T-T With 22 depth levels, it worked. I would say that the results are similar in the way that the transition is still visible but it is less obvious since the gap between 2 successive layers is smaller so the transition is closer to the layer boundary.

I don't know if it is a big problem since it is not that visible (but I could also see the same effect at the bottom layers of the Atlantique domain) and no one said anything about it during the previous EMODNET phases (we checked and the artifact was already there on the Loire DIVAnd maps of phase 4 but had remained unnoticed...)

In case you are interested, I'll join to the current message

• the netcdf file with the observations used
• the bathymetry file
• the notebook test_Loire.zip