Negative Salinity for Seawater Intrusion Model and Instability Issues

I am trying to simulate the seawater interface in a highly conductive aquifer, yet get massive numerical instability and negative mass concentrations!

I have tried to include the upwinding, but need a sharper mixing zone (i.e. lower dispersivity?) to even closely match the collected data.
If I understand it correctly, upwinding will increase the dispersivity to reduce the gradient over the cells and should resolve some of the instability issues?
The other option I have seen recommended is to refine the mesh more, however my model already has 1.5 million nodes with element diameter of < 0.4 metres, and Peclet number is less than 0.5 on average, and I can't wait all year for even finer models to run...

so there are two questions;
[list type=decimal]
[li]How can I stop getting negative solute concentrations (i.e. is there a hard limit available?)[/li]
[li]How do people [b]realistically [/b]model saline water intrusion in reasonable time-frames?[/li]
[/list]

Any real-world examples would be greatly appreciated - I'm struggling to find anything in the literature...

There is no hard limit available to stop negative concentrations, because the amount and magnitude of negative concentrations depend on model specific settings.

A short discussion in our forum is provided here: http://forum.mikepoweredbydhi.com/index.php?topic=2263.msg5187#msg5187

A more general, but interesting and critical discussion is provided by the White Paper II, Chapter 1. You may find this and other White Papers on the DHI website: https://www.mikepoweredbydhi.com/download/product-documentation

Great - Thanks Bjorn.

I am fairly confident I have covered all of the points mentioned in that thread - I will cover them here again for posterity..

[list][li][b]Is the mesh fine enough?[/b][/li][/list]
I think so! Assuming the Peclet number is a valid measure of mesh approriate-ness (I think we discussed in another thread).
[list][li][b]Does the mesh contain bad shaped elements?[/b][/li][/list]
0.0% Delaunay-violating triangles
6.5% > 90 degrees
0.0% > 120 degrees
[list][li][b]What temporal discretization are you using?[/b][/li][/list]
Automatic, Maximum of 1 day.
[list][li][b]Fully implicit schemes?[/b][/li][/list]
FE/BE
[list][li][b]What error norm are you using?[/b][/li][/list]
First.

The only way I've been able to stabilise it is using the upwinding, and I've used the PGLS (least-squares) one for no other reason that it seems to work best (I didn't have much luck with Shock upwinding).

There are just so-many nodes now, it's incredibly slow to run - and it's not even 3D yet...!

[quote]There is no hard limit available to stop negative concentrations, because the amount and magnitude of negative concentrations depend on model specific settings. [/quote]

Could you elaborate on the model-specific settings...? Is this mesh-related or another parameter?

Model-specific settings involve all settings you adopt in the Problem Settings, the mesh, type of boundary conditions, contrasts in material properties etc. In fact, the settings are related to your individual model. There is no way to provide general guidelines. However, if you like you could contact the FEFLOW Support. The support can have a have a look on your model and provide suggestions on how to improve it. https://www.mikepoweredbydhi.com/support 

Support's not willing to help due to academic licensing restrictions...

I guess I get to keep begging here in the forums!