I would probably need to sit down and thing about this some more. Also, I have never used delft-3d, so I can only guess the shallow water equations it solving. However, I will try to give you a few points to consider.

1) It sounds like you are modeling fully turbulent, 3D, flow using k-e in Fluent and then trying to compare it to an analytical relationship and the shallow water result. Given this, do you have a good reason to believe Fluent should be converging on the same answer as the analytical and shallow water solution in the first place? Did the analytical study include turbulent fluxes? What about the shallow water equations?

2) Is the actual flow you are trying to model fully turbulent? As far as I know, Fluent doesn't have a transition model for it's VOF algorithm. Therefore, if the actual flow isn't fully turbulent, the k-e model is going to make it fully turbulent regardless (it may be beneficial to switch to a one equation turbulence model like Spalart if not fully turbulent). This could easily give false mass fluxes in all directions.

3) Is the actual flow you are trying to model steady?

4) Does Fluent offer LES in conjunction with VOF? Might be worth a shot for unsteadiness and/or better flux calculations.

5) Have you done a grid resolution study at all? Are the Fluent fluxes still under-predicted if the grid resolution gets finer?

6) And finally, as a CFD pro, it is fairly well known throughout the community that Fluent is shite

. Especially when it comes to resolving steep gradients (like what you see between an air/water interface). Fluent is going to smooth, diffuse, wash away, flux correct, and numerical viscosity its way through any steep gradient it feels will go unstable and crash the code. For them, the wrong answer is much better than NaN.