Mimicking MFUSG solver and properties approaches for Unconfined Aquifers with MF6

[javgs-bd]

I'm trying mimic as much as possible the settings we have been using for MFUSG models with MF6 using Flopy, for a model where the water table crosses multiple layers, and will possibly change significantly between natural and stressed scenarios. It also includes multiple combinations of highly non-linear BCs (such as depth dependent evap).

I've been trying to map all the parallels that need to be checked between the LPF/NPF/STO packages and the SMS and IMS solvers. I'm clear that the laytype option of the LPF needs to be replaced through icelltype in NPF and iconvert in STO. I can also correlate the majority of the IMS options with those of the SMS. There are some variables that are bit puzzling though, because I'm not clear on how redundant they are and what considerations should be made. These are the following:

• newtonoptions in the ModflowGwf object and under_relaxation in mfims. I would guess both of them need to be activated at the same time if one is under relaxing in the IMS solver. Is there a reason to have them in two places? Or what happens if only one of them is activated?.
• SOLVEACTIVE and DAMPBOT seem to have gone without analogies in IMS. No problem with that.
• dev_modflowusg_upstream_weighted_saturation in NPF. I would understand that by activating this the NPF horizontal conductance weighting will follow the principles of MFUSG LPF. Would this overwrite any setting of the alternative_cell_averaging in case it would have been setup??
• cvoptions and perched in the NPF. I would understand that these would relate to CONSTANTCV and NOVFC of the LPF (I have the MF6 manual checking pending).

Beside wondering if I'm right on those specific question, I guess I could narrow down the main question to:
if a LAYTYP of 4 is used for layers in a MFUSG model plus CONSTANTCV and NOVFC are activated in the LPF package. Should I rely that dev_modflowusg_upstream_weighted_saturation, cvoptions and perched are the variables I need to setup to reproduce the conductance calculation scheme?

I would greatly appreciate any shared ideas on this...
Regards
Javier

[javgs-bd]

Ok, when running MF6 it complains if cvoptions='DEWATERED' or perched='True' are used with the Newton Raphson scheme activated. I guess dev_modflowusg_upstream_weighted_saturation and the IMS settings are the ones doing the trick. I better go and check the IMS docs..

[javgs-bd]

Just in case somebody is looking into the same, also dev_modflowusg_upstream_weighted_saturation in the NPF package is an invalid option with MF6 6.4.2

[hjia1005]

@javgs-bd I am very interested in this because I recently did the same - to translate a mf6-disu model to a mfusg model to reduce the runtime. But I didn't go as deep in detail as what you did.

• I used equivalent packages and data for both models, such as IC - BAS6, NPF/STO - LPF, TVK/TVS - TVM, etc.
• I used the same TDIS settings. Since this is a recovery model I don't have a convergency issue. To simplify, I removed all ATS options and included only one timestep in each stress period.
• I used icelltype = 1 for all cells in mf6 model and LAYTYP = 4 for all layers in mfusg model.
• I used Newton-Raphson Under_relaxation in both models.

The mfusg model does speed up significantly. However, to my surprise, there are large discrepancies in both heads and water budgets between the two models, especially in the areas with packages that cause big hydraulic gradients. Would like to know if you have experienced the same.

[langevin-usgs]

I'd be curious to know if @javgs-bd or @hjia1005 find the cause of the differences between mf6 and mfusg. We've done quite a bit of testing on our end and find the performance to be similar. The solvers and solver settings have some differences and some convergence checks may be different as well, but the underlying equations for flow should be very close, especially when using the DEV_MODFLOWUSG_UPSTREAM_WEIGHTED_SATURATION option in mf6. Note that this option is available in mf6 but you have to use the "develop" version. This dev option is not available in the official release. You can get the most recent development version of mf6 here.

[hjia1005]

@langevin-usgs just a quick report back here. I added the DEV_MODFLOWUSG_UPSTREAM_WEIGHTED_SATURATION option to the NPFpackage of mf6 model to mimick the upstream weighting I used in mfusg model (LAYTYP = 4 for all layers). But now the model is struggling to converge, still stuck at the first stress period after 1.5 hours. Looking at the listing file I don't see any possibility that it will converge. The model converges well without the DEV_MODFLOWUSG_UPSTREAM_WEIGHTED_SATURATION .

Below are the options I used in NPF and LPF respectively for you to check. If needed, I can share the model files with you privately via my work email, my boss has approved this.

# File generated by Flopy version 3.5.0.dev0 on 12/22/2023 at 09:26:32.
BEGIN options
  SAVE_SPECIFIC_DISCHARGE
  TVK6  FILEIN  model.npf.tvk
  DEV_MODFLOWUSG_UPSTREAM_WEIGHTED_SATURATION
END options

# LPF package for MODFLOW-USG generated by Flopy 3.5.0.dev0
        50    -1E+30         0         0  
         4         4         4         4         4         4         4         4         4         4         4         4         4         4         4
         0         0         0         0         0         0         0         0         0         0         0         0         0         0         0
         1         1         1         1         1         1         1         1         1         1         1         1         1         1         1
         0         0         0         0         0         0         0         0         0         0         0         0         0         0         0
         0         0         0         0         0         0         0         0         0         0         0         0         0         0         0

[javgs-bd]

Hello, and thanks @hjia1005 and @langevin-usgs for following up on this thread.

To be rigorous, I have not really begun the side-by-side comparison journey of mf6 and mfusg. I've changed to mf6 in our most recent model because the ability to activate BCs (actually it was the HFB pkg) in different stress periods at will was a necessity. Also the flopy full support has become a bliss.

So the only contribution I can make to this discussion are impressions on the ability of the IMS solver to keep the model stable. In this regard, I'm actually impressed on how stable is behaving, but this has only happened after I dropped the attempt to mimick mfusg. The key step for me was to drop the DEV_MODFLOWUSG_UPSTREAM_WEIGHTED_SATURATION option.

I followed Sorab's solvers course recommendations on the solver tweaking. The following IMS setup is working for me in a really challenging steady state setup. I know these setting are case specific but maybe they can be tested and see if convergence is improved. They have also been defined following the flopy and MF6 manual documentation.

ims = flopy.mf6.ModflowIms(
                            sim,
                            print_option                    = 'SUMMARY', 
                            csv_outer_output_filerecord     = 'ims_outer.csv', 
                            csv_inner_output_filerecord     = 'ims_inner.csv',                             
                            outer_dvclose                   = 1e-4, 
                            outer_maximum                   = 600,
                            under_relaxation                ='DBD', 
                            under_relaxation_gamma          = 0.001,
                            under_relaxation_theta          = 0.9,
                            under_relaxation_kappa          = 0.001,
                            under_relaxation_momentum       = 0.0,
                            backtracking_number             = 200,
                            backtracking_tolerance          = 1.5, 
                            backtracking_reduction_factor   = 0.9,
                            backtracking_residual_limit     = 1,
                            inner_maximum                   = 600,
                            inner_dvclose                   = 1e-5, 
                            rcloserecord                    = [1e-5, 'RELATIVE_RCLOSE'],
                            linear_acceleration             = 'BICGSTAB', 
                            relaxation_factor               = 0.97,
                            preconditioner_levels           = 7,
                            preconditioner_drop_tolerance   = 1e-7,
                            number_orthogonalizations       = 10,
                            scaling_method                  = 'NONE',
                            reordering_method               = 'NONE', 
                            pname                           = root_name
                          )

If I'm not wrong the backtracking residual limit at 1 might be unnecessarily low but I have run out of time for testing. Somehow I got to extremely low values of gamma and kappa.
The model has gone from more than 10 minutes to only 2. Levels and north settings have been key for linear solver convergence and rcloserecord, have helped. I've also seen improved robustness by lowering the drop tolerance from 1e-5 to 1e-7, I guess that matrix needed those terms. I've also tried different reordering methods.

I have to say that Sorab's course have been extremely valuable in knowing what to try and why in the solver settings. Now the IMS documentation is another major help.

I hope this is useful.

Regards!

[langevin-usgs]

Thanks for the updates @javgs-bd and @hjia1005!

If you are using Newton with a steady state stress period, the PTC capability in IMS in MODFLOW 6 in on by default. We've generally found PTC (psuedo-transient continuation) to be effective at improving convergence for steady state stress periods. However, there are cases where it can slow convergence. MODFLOW-USG does not have a PTC capability, so this can be another important difference when comparing MODFLOW-USG and MODFLOW 6 solution behavior. You can turn off PTC for MODFLOW 6 using the NO_PTC option in IMS. This may be worth trying.

[hjia1005]

@langevin-usgs I recently had another comparison between mf6 and mfusg for your information.

It's an 8-layer model with full mesh for each layer. Used Gridgen to create the mesh and applied refinement up to level 5 on some important features. Extensive pilot points were used in each layer and the structural overlay from PLPROC was applied to introduce some faults into the model. Pilot points were also applied on the vertices of the fault to adjust both the hydraulic properties and the positions of the vertices. I specially used very high hydraulic property values to the fault to see the difference.

The model was initially developed as mf6 with DISV grid, it converged well with the initial uniform hydraulic properties per layer but was a bit slow. Then tried mfusg. It converged too and the runtime was around 1/3 of the mf6 model. Therefore I decided to use mfusg until the pest run (used pestpp-ies) was kicked off. NONE of the initial realizations could converge with many different SMS settings collected from the experts.

I had to change it back to mf6 and surprisingly convergence was not an issue at all no matter how crazy the parameter field is. Because it's the initial pest run, I left the parameter ranges quite wide.

I heard that mfusg does not work well with quadtree grids, can confirm it from my own experience. By the way, the two converged models with the same initial parameter setting also had quite different results. That could be partly due to the different solver settings.

[langevin-usgs]

Hello @hjia1005, thanks for the information. I'm surprised that the MODFLOW 6 and MODFLOW-USG models behave so differently. While they do have some underlying differences in their formulations and solver capabilities, the numerics in the programs are more similar than they are different. Let us know if you ever discover why the behavior is so different between the two. Thanks!