Dynamic Update of Hydraulic Conductivity During Simulation Has No Effect (MODFLOW API)

[zongxunchai]

Description

I am working on reactive transport modeling where chemical reactions alter the hydraulic conductivity (K) of the porous medium during simulation.

To explore this, I coupled MODFLOW 6 with PhreeqcRM, verified several PHT3D benchmark cases, and attempted to dynamically update k11 using the modflowapi interface. However, modifying k11 during runtime does not appear to affect the simulation results.

Minimal Example

A minimal reproducible 1D example:

import time
import flopy
import numpy as np
import modflowapi
import matplotlib.pyplot as plt

initial_k = np.ones((1, 1, 50)) * 10.0

sim = flopy.mf6.MFSimulation(sim_name='model', sim_ws='./mymodel', exe_name='mf6')
tdis = flopy.mf6.ModflowTdis(sim, nper=1, perioddata=[(365, 100, 1.0)])
ims = flopy.mf6.ModflowIms(sim)
gwf = flopy.mf6.ModflowGwf(sim, modelname='gwf_model', save_flows=True)
dis = flopy.mf6.ModflowGwfdis(gwf, nrow=1, ncol=50)
ic = flopy.mf6.ModflowGwfic(gwf)
npf = flopy.mf6.ModflowGwfnpf(gwf, k=initial_k, save_specific_discharge=True)
chd = flopy.mf6.ModflowGwfchd(gwf, stress_period_data=[[(0, 0, 0), 1.], [(0, 0, 49), 0.]])
oc = flopy.mf6.ModflowGwfoc(gwf,
    budget_filerecord='gwf_model.bud',
    head_filerecord='head.hds',
    saverecord=[('HEAD', 'ALL'), ('BUDGET', 'ALL')])

sim.write_simulation()
sim.run_simulation()

Then, attempting to update k11 dynamically:

modflow_api = modflowapi.ModflowApi("C:\\ProgramFiles\\MODFLOW\\libmf6.dll", working_directory="./mymodel")
modflow_api.initialize()
sim = modflowapi.extensions.ApiSimulation.load(modflow_api)
gwf_model = sim.get_model('gwf_model')
npf = gwf_model.get_package("npf")

hk = npf.get_array("k11")
current_time = modflow_api.get_current_time()
end_time = modflow_api.get_end_time()

while current_time < end_time:
    dt = modflow_api.get_time_step()
    modflow_api.prepare_time_step(dt)
    modflow_api.update()

    hk[0, 0, 0:10] += 1  # increase K in first 10 cells
    npf.set_array("k11", hk)

    current_time = modflow_api.get_current_time()

modflow_api.finalize()

Expected Behavior

Changes to k11 during runtime should alter flow results (e.g., head distribution).

Actual Behavior

The simulation results remain identical regardless of modifications to k11.

Additional Context

The following statement from Hughes et al., 2022 (Environmental Modelling & Software) may be relevant:

“… modifying k11 after initialization would have no effect on the simulation results. Other parameters, such as those related to the dimensionality of the coefficient matrices or to the time discretization, should probably not be modified to avoid undefined behavior and possible program failure.”

This suggests that k11 may be fixed after initialization, preventing runtime updates.

Questions

Is it currently possible to dynamically update hydraulic conductivity during a simulation using modflowapi (e.g., within each time step)?

Has any new mechanism been introduced—perhaps via the Time-Varying Hydraulic Conductivity (TVK) capability or another feature—that would allow K to evolve dynamically according to an external process (e.g., chemical reactions computed by PhreeqcRM)?

If not, is there any recommended workaround for coupling reactive transport processes that alter flow properties mid-simulation?

Any guidance or documentation references would be greatly appreciated.

[emorway-usgs]

Hello @zongxunchai, you will want to try updating K11 through the TVK package. Through the use of the TVK package, when an updated value of K11 is detected, it will trigger a re-calculation of inter-cell conductance values. Accessing K11 directly through the API won't trigger the recalculation of inter-cell conductance values.

[zongxunchai]

Dear @emorway-usgs , Thank you for your suggestion. I explored the TVK package in detail, but unfortunately TVK does not appear to support the type of dynamic property update I need.

In my application, hydraulic conductivity (K) is not known in advance and cannot be prescribed as a time-series. Instead, K evolves during the simulation based on geochemical reactions computed by PhreeqcRM at each time step. This requires updating K11 in response to external processes at runtime.

However, TVK does not currently enable dynamic, externally-driven K updates such as those arising from reactive transport coupling.

[christianlangevin]

Another thought, which might be a workaround for limited use cases, would be to simulate the problem as being unconfined (icelltype=1 and h<top) and NOT use the Newton formulation. I think in this particular case, conductance would be calculated directly from K11 in the NPF package, which you could change using the API. It might be worth a shot, though I recognize this is not ideal.

[emorway-usgs]

@zongxunchai, hopefully I'm not sending you down a dead-end path. I understand that you don't know the reaction-dependent K11 values you want to use in advance, but you still need to enter initial values that MODFLOW will read even if you manually override those values through the API later in the simulation.

If, in attempting to test what I was suggesting, your test model(s) didn't show any difference between a base run and a "scenario" run with updated K11 values, it might be because in addition to altering K11 through the API, you may also need to update a variable called nodekchange (see line 430 in screen shot below). By setting nodekchange equal to 1 (presumably a flag that alerts the code to a user-updated change in K11) for the affected cells the code will call the calc_condsat(..) to update the appropriate intercell conductance values stored in memory that otherwise remain static throughout the simulation.

My apologies if I'm still not understanding your specific use case. I think the recent shutdown contributed to mental atrophy. I'd sure like to see your MODFLOW 6/PhreeqcRM coupling succeed.

[zongxunchai]

@emorway-usgs @christianlangevin Thank you both for your valuable insights and suggestions.

I am happy to report that your advice regarding the internal flags was the key to solving the problem. Following your guidance, I used the modflowapi to set nodekchange to 1.

Additionally, I found that I also needed to update kchangestp and kchangeper to the current time step and period within the simulation loop. After implementing these changes, the hydraulic conductivity is now updating dynamically, and the conductance values are successfully being recalculated during runtime.

The coupling between MODFLOW 6 and PhreeqcRM is now working as expected. Thank you again!

[emorway-usgs]

Happy to hear your PhreeqcRM-MODFLOW 6 coupling is working out. As you get further along, if you have a small working example that you'd be willing to share with the MF6 community, I invite you to zip it up and attach it here. I've little doubt others will be interested to see how you did the coupling.

I'm glad your persistence has returned some early dividends!

[zongxunchai]

Hello @emorway-usgs , I am more than happy to share the working example with the community.

Actually, I have already uploaded the coupled code to a GitHub repository, which can be found here:

https://github.com/wangzitao21/mf6pqc

I have named this MODFLOW 6 and PhreeqcRM coupling script MF6PQC. I have tested it against several PHT3D benchmark cases, as well as the examples provided in this paper: Xie et al., 2015.

So far, these cases seem to run successfully with correct results. However, I have currently only implemented a loose coupling approach using the update() method; I haven't yet managed to figure out a way to implement tight coupling.

Please feel free to take a look. Since the code is still under active development and limited by my current coding proficiency, there are likely many areas that could be improved. I would be incredibly grateful if you have any advice or guidance to offer.

Just a quick note on the username mismatch: @wangzitao21 is my main account (Zitao Wang is my real name), whereas @zongxunchai is just a random alias I use.