Non-linear sorption in IST

[p-ortega]

Hey dev team,

I'm wondering if there are any short-term plans to add non-linear sorption to the dual-domain IST package.

Following on that, can sorption be 'on' in the MST but 'off' in the IST?
Currently, I get the following error with soprtion = "LINEAR" in MST and sorption = None in IST

Sorption is active for the IST Package but it is not compatible with the sorption option selected for the MST Package. If sorption is active for the IST Package, then SORPTION LINEAR must be specified in the options block of the MST Package.

Thanks,
P

[p-ortega]

I reckon this kind of answers my second question

https://github.com/MODFLOW-USGS/modflow6/blob/develop/src/Model/GroundWaterTransport/gwt-ist.f90?plain=1#L190

[langevin-usgs]

I think you can shut off sorption in the immobile domain (but leave it on in the mobile domain) by specifying zero for the distribution coefficient in the IST Package. Have you tested to see if that works?

Regarding nonlinear sorption for IST, I think we would consider adding it. It might be a few weeks before we get to it. Is this important for your problem?

[p-ortega]

Hi @langevin-usgs ,

Thanks for the quick response.

Yeah, that's an option that works with linear sorption in MST and IST! Thanks for that.

My issue now is that I have (unwillingly 😅 ) to do the same, but using non-linear sorption in MST but still have IST on (with or without sorption).
A workaround could be to read the concentrations with the API and run the non-linear sorption outside of mf6 ..

As for the non-linear sorption in IST: it's not really key for my model 😄 , but I thought I could test it if it were close to deployment, so no worries there.

Thanks again!

[p-ortega]

Hi @langevin-usgs,

It's me again, haha. I kind of failed on this, so...

Could I kindly ask if the team could look into implementing non-linear (Freundlich) sorption in IST, or possibly having the option to run MST with nonlinear sorption and IST without?

I tried handling this via the API by calculating sorption in the outer loop, but without much success when comparing the results with mf6 (got kind of close, but mf6 struggles to converge).

At the moment, I have a dodgy solution (at the risk of opening a can of worms) where I cloned the mf6 repo, commented out the part where IST checks if the sorption definition matches MST, and then recompiled it. It runs okay-ish, but I suspect some mass balance/conservation issues could arise from this ad-hoc solution?

Thank you,
P

[langevin-usgs]

Hi @p-ortega, I will take a look and see what it would take to implement the non-linear isotherms in the IST package. If it seems straightforward, I'll see if I can get it implemented this week.

[langevin-usgs]

There is a pull request up at #2027 that has the nonlinear Freundlich and Langmuir isotherms implemented for the immobile domain (using the IST Package). There is more work to do on it, but if you are capable of pulling down the pull request branch, and building MODFLOW 6, it would be helpful if you could give it a try while I work on some of the documentation. No problem if you don't get to it; I will work to finish it up.

After checking out the pull request branch on my personal fork, you will need to run the update_flopy.py script in the modflow6/autotest folder. This will make flopy aware of the new options. You will need to run the update_flopy.py script using python from whatever environment you typically work in. The modflow6/autotest/test_gwt_ist03.py has the syntax for using the new capability. Note that there is also a new option to save the immobile domain sorbate concentrations.

[p-ortega]

Awesome @langevin-usgs ! Thank you so much .. I'll give it a try and come back to you.

[p-ortega]

Hi @langevin-usgs

I can confirm that this works in my model! Thank you!

I am also having a harder time achieving convergence (or not converging) for one of my components, as mentioned in the PR 2027. My current IMS setup is:

    nouter, ninner = 500, 400 
    rclose, relax = 1e-6, 1.0
    inner_dvclose = 1e-6
    outer_dvclose = 1e-5 
    under_relaxation = "dbd"
    under_relaxation_gamma = 0.2 # thumb suck
    under_relaxation_theta = 0.7

    ims = flopy.mf6.ModflowIms(
        sim,
        complexity = 'COMPLEX',
        print_option="SUMMARY",
        outer_maximum=nouter,
        under_relaxation=under_relaxation,
        # under_relaxation_gamma = under_relaxation_gamma,
        inner_maximum=ninner,
        inner_dvclose=inner_dvclose,
        outer_dvclose=outer_dvclose,
        rcloserecord=rclose,
        under_relaxation_theta=under_relaxation_theta,
        linear_acceleration="BICGSTAB",
        scaling_method="NONE",
        reordering_method="NONE",
        relaxation_factor=relax,
        filename=f"{model_name}.ims",
    )

[langevin-usgs]

Glad to hear it's working. Will close this. About the convergence -- well that's a downside. As a long shot you could try decreasing relax down from 1.0 to .9 and all the way down to 0. I suspect your convergence problems are in the outer iterations so probably futzing with the DBD settings are your best bet. Could also decrease the number of inner iterations. 400 seems like a lot for what I'm guessing is a rather nonlinear problem.

We are finalizing the documentation and other fine details on #2027 but it should get merged into develop soon, and will be released with the next version. Keep me posted on performance and lessons learned.