Numerical dispersion in a simple rotating interface model (GWT+BUY)

[ifariasg]

Together with @dzamrsky we have been investigating how MF6+BUY results compare to SEAWAT due to a general idea that MF6 (GWT specifically) tends to show more numerical dispersion. Since UTVD was added to MF6 v6.7.0 we thought this was acknowledged and improved so we started migrating some existing models to MF6 to compare.

For more advective problems results seem to be quite good, but when looking at more density-driven problems we noticed discrepancies. For example, we tested a simplified version of the rotating interface problem.

We found that when dispersion parameters are set to 0, there area noticeable differences in concentration results, with MF6 showing consistently more numerical (?) dispersion. No matter what setup I used I couldn't improve the results. Some of the things I tried:

• Change solver parameters
• Reduce time step size
• Use ATS with very small time step size
• Change initial head formulation
• Change (or remove) the CHD added for stability
• Use the RHS formulation for BUY
• Change the advection scheme (UTVD,TVD,central)

What's also worrying, the MF6 version was ~10 times slower to run, no matter the settings.

At this stage I can only think of two possibilities, there is an option I didn't implement that makes both models essentially different, or there is something strange going on internally. To spark some discussion, I'm sharing the script used to see if anybody has an idea on what could be causing this. Also it would be interesting to know if anybody has been experiencing similar problems.
benchmark_seawat_vs_mf6.py

This is the comparison with no dispersion

...and strangely enough, with dispersion the results look a lot more alike...

[mjr-deltares]

Hi @ifariasg , thanks for starting up this discussion, and for posting the script along with the findings. It's very helpful for us to further improve the code. We (the development team, @Manangka) will definitely have a closer look at this. But since the Groundwater Modeling and More conference in Princeton is coming up, it might take a few weeks until you'll hear from us :-).

Just two things for now. You mention that for more advective problems, the results are quite good. Can you share what type of problems you have looked at? And about the performance results, A quick look in the listing file tells me that there are convergence issues with the transport model. I didn't have time to go into it any further, but I think it will distort the run time performance numbers significantly. I guess we will need to straighten that out first.

[ifariasg]

Hi @ifariasg , thanks for starting up this discussion, and for posting the script along with the findings. It's very helpful for us to further improve the code. We (the development team, @Manangka) will definitely have a closer look at this. But since the Groundwater Modeling and More conference in Princeton is coming up, it might take a few weeks until you'll hear from us :-).

Just two things for now. You mention that for more advective problems, the results are quite good. Can you share what type of problems you have looked at? And about the performance results, A quick look in the listing file tells me that there are convergence issues with the transport model. I didn't have time to go into it any further, but I think it will distort the run time performance numbers significantly. I guess we will need to straighten that out first.

Hi @mjr-deltares,

As always, thanks for your answer.
The advective problem I'm referring to is a small clipped area from the national fresh-salt model of the Netherlands where I test an extraction/reinjection pair of wells. Now that I think of it, this model also has (mechanical) dispersion, so the fact that its advective might not be the reason both models look comparable.

This is a comparison of that "more advective problem" I tested for reference. The difference scale is much narrower, I would consider them reasonably equal. The top row has the SEAWAT version, MF6 + TVD and MF6 +uTVD. Second row are the differences assuming SEAWAT is the "base"

Curious to hear if anybody else at the conference talks about similar issues!

[Manangka]

Hi @ifariasg.

Thank you for sharing your model and great to hear that for some of your problems the new utvd scheme shows improved results.
I've recreated your model in a notebook. I removed the dispersion stuff to simplify things and really focus on the numerical dispersion.
There were some changes i needed to make. For instance, SEAWAT by default takes more timesteps for the transport model than for the flow model. To be able to compare it to MODFLOW i fixed the transport timestep to be equal to the flow timestep.
rotating interface.ipynb

The first thing i noticed is that with the notebook the SEAWAT and MODFLOW takes about the same amount of time (around 2s on my machine). I haven't dived into the reason why there is a difference in runtime between your script and my notebook

The second thing is that I indeed see that there is a difference between the SEAWAT concentration and that of MODFLOW. The SEAWAT model has less numerical dispersion. A good way to see this is by plotting the contour lines of the concentration. The width between them says something about the spread.

The reason why MODFLOW displays more numerical dispersion is due to the internal limiter we use under the hood. By default this is fixed to use the van Leer flux limiter. This limiter is known for its stability but is also more diffusive than some of the other limiter. When I switch to a more compressive limiter, like the Superbee limiter, the results of SEAWAT and MODFLOW are comparable. Unfortunately this setting is not something you can set through an option though but it is hardcoded.

Van Leer

Superbee

And the frames at the last time step

Van Leer

Superbee

[ifariasg]

Hi @Manangka,

Thanks for the thorough explanation. I can't but feel a bit worried about your findings. Is this something that the user of MF6 should be worried about? Or did I just catch an edge case? Is there any pointers about when this could be a problem in real world scenarios? My mind immediately thinks in models with very sharp (density driven) concentration interfaces but I'm not sure if it would manifest like in this example.

[ifariasg]

As a side note,
Regarding runtimes, I managed to find why my Mf6 model was much slower than SEAWAT:

1. Changing a bit solver parameters (based on @Manangka's notebook) brought the runtime down to 3-4x times slower. I was missing inner_dvclose.
2. In my script I activated the dispersion package always, even when all its parameters are set to zero. In MF6 by default, this activates the XT3D formulation for transport, which in this case slows the models by another factor of 3. This can be turned off with the flag xt3d_off=True in the dispersion package.

With those two changes both MF6 and SEAWAT have comparable runtimes.