Wrong results with models that use time series

[pya]

Approach

I ran all models from MODFLOW 6 examples: built 02/06/2025 19:17.

Code

This code was used to run all models:

from modflowapi import Callbacks, run_simulation

def callback(sim, callback_step):
    pass

run_simulation(
    dll='path/to/libmf6.dylib',
    sim_path='path/tomodel',
    callback=callback,
)

Versions

• modflowapi: 0.3.0.dev0
• xmipy: 1.3.1
• MODFLOW: 6.6.1

Results

Most model runs produce the same output compared with the standalone MODFLOW executable. But these models:

• ex-gwf-csub-p03a
• ex-gwf-fhb
• ex-gwf-advtidal
• ex-gwf-csub-p03b

produce significantly different results in the balance values as listed in each <mode-name>.lst file. All of these models use of time series.

These models also use time series:

• ex-gwf-csub-p01
• ex-gwf-sfr-pindersauera
• ex-gwf-sfr-pindersauerb

but show no differences in results compared to the standalone runs.

Potential problem

The problems seems to be in the solution loop in modflowapi.extensions.runner.run_simulation().

Replacing this loop:

    for sol_id, slnobj in sorted(sim.solutions.items()):
            models = {}
            maxiter = slnobj.mxiter
            solution = {sol_id: slnobj}
            for model in sim.models:
                if sol_id == model.solution_id:
                    models[model.name.lower()] = model

            sim_grp = ApiSimulation(
                mf6, models, solution, sim._exchanges, sim.tdis, sim.ats
            )
            mf6.prepare_solve(sol_id)
            if sim.kper != kperold[sol_id - 1]:
                callback(sim_grp, Callbacks.stress_period_start)
                kperold[sol_id - 1] += 1
            elif current_time == 0:
                callback(sim_grp, Callbacks.stress_period_start)

            kiter = 0
            callback(sim_grp, Callbacks.timestep_start)

            if sim_grp.ats_period[0]:
                mindt = sim_grp.ats_period[-1]
                while sim_grp.delt > mindt:
                    sim_grp.iteration = kiter
                    callback(sim_grp, Callbacks.iteration_start)
                    has_converged = mf6.solve(sol_id)
                    callback(sim_grp, Callbacks.iteration_end)
                    kiter += 1
                    if has_converged and sim_grp.allow_convergence:
                        break

            else:
                while kiter < maxiter:
                    sim_grp.iteration = kiter
                    callback(sim_grp, Callbacks.iteration_start)
                    has_converged = mf6.solve(sol_id)
                    callback(sim_grp, Callbacks.iteration_end)
                    kiter += 1
                    if has_converged and sim_grp.allow_convergence:
                        break

            callback(sim_grp, Callbacks.timestep_end)
            mf6.finalize_solve(sol_id)

with:

mf6.do_time_step()

fixes the problem and all models produce the same results as the standalone runs.

[jlarsen-usgs]

@pya

This is something that I'll need to explore more to see what's going on. Unfortunately, we can't simply remove/replace the solve loop with do_time_step() in the extensions source code because it will remove the callback entry points from user access. Instead, it looks like we'll have to find out 1) which modflow function isn't being called and 2) create either an api call that executes the function or write some python code to mirror the function in the api extensions.

It'll go on the list of needed improvements. Thanks for identifying this issue!

[pya]

@jlarsen-usgs

The "do_time_step() shortcut hack" is more a debugging tool than a solution. I am not sure if this could be more a problem in xmipy than in modflowapi.

[jlarsen-usgs]

Tested the listed timeseries models with modflow version 6.7.0.dev1 and the corresponding libmf6.dll on windows 11. Was not able to reproduce the issue. List budget outputs for models run with mf6.exe and modflowapi.run_simulation() had identical budget volumes.

[pya]

Sounds good. Maybe the problem disappeared. I will test with the latest libmf6 version on macarm. What version of xmipy did you use? I use 1.3.1 from conda-forge.

[jlarsen-usgs]

I'm using xmipy version 1.3.1.

Here's the code I check the issue with. I visually inspected the cumulative and incremental budget results, but that could also be done with np.allclose().

def callback(*args, **kwargs):
    pass


if __name__ == "__main__":
    mname = "ex-gwf-csub-p03b"
    sim_ws = Path(f"./{mname}")
    new_ws = Path(f"./debug_ts")
    # for addr construction
    mname = mname.upper()

    fsim = flopy.mf6.MFSimulation.load(sim_ws=sim_ws)
    fsim.set_sim_path(new_ws)
    fsim.write_simulation()
    fsim.run_simulation()
    fgwf = fsim.get_model()

    heads = fgwf.output.head()
    list = fgwf.output.list()
    vinc, vcum = list.get_dataframes()


    run_simulation(dll="libmf6", sim_path=str(sim_ws), callback=callback, verbose=True)

    asim = flopy.mf6.MFSimulation.load(sim_ws=sim_ws)
    agwf = asim.get_model()
    alist = agwf.output.list()
    ainc, acum = list.get_dataframes()

[pya]

I still get differences for these models:

• ex-gwf-csub-p03a
• ex-gwf-fhb
• ex-gwf-advtidal
• ex-gwf-csub-p03b

I run with these versions:

xmipy version:	1.3.1
modflowapi version:	0.3.0.dev0 (develop May 20, 2025)
MODFLOW version:	6.7.0.dev2 (nightly May 20, 2025)

Let's take the model ex-gwf-fhb as an example.
I get these differences between the ex-gwf-fhb.lst of both model runs:

	mf6	modflowapi
lineno
ex-gwf-fhb.lst_351	STO-SS = 0.0000 STO-SS = 0.0000 STORAGE	STO-SS = 3.6826E-12 STO-SS = 0.0000 STORAGE
ex-gwf-fhb.lst_352	WEL = 1850572.8592 WEL = 5807.8512 WEL	WEL = 800000.0000 WEL = 2000.0000 WEL
ex-gwf-fhb.lst_355	TOTAL IN = 1850572.8592 TOTAL IN = 5807.8512	TOTAL IN = 800000.0000 TOTAL IN = 2000.0000
ex-gwf-fhb.lst_359	STO-SS = 100.0215 STO-SS = 0.2045 STORAGE	STO-SS = 18.0000 STO-SS = 3.2399E-12 STORAGE
ex-gwf-fhb.lst_361	CHD = 1850472.8377 CHD = 5807.6467 CHD	CHD = 799981.9994 CHD = 2000.0000 CHD
ex-gwf-fhb.lst_363	TOTAL OUT = 1850572.8592 TOTAL OUT = 5807.8512	TOTAL OUT = 799999.9994 TOTAL OUT = 2000.0000
ex-gwf-fhb.lst_365	IN - OUT = -2.6643E-06 IN - OUT = -2.3920E-10	IN - OUT = 5.9481E-04 IN - OUT = 1.5929E-06
ex-gwf-fhb.lst_367	PERCENT DISCREPANCY = -0.00 PERCENT DISCREPANCY = -0.00	PERCENT DISCREPANCY = 0.00 PERCENT DISCREPANCY = 0.00
ex-gwf-fhb.lst_447	STO-SS = 0.0000 STO-SS = 0.0000 STORAGE	STO-SS = 3.6826E-12 STO-SS = 0.0000 STORAGE
ex-gwf-fhb.lst_448	WEL = 2960490.2146 WEL = 5394.6281 WEL	WEL = 1200000.0000 WEL = 2000.0000 WEL
ex-gwf-fhb.lst_451	TOTAL IN = 2960490.2146 TOTAL IN = 5394.6281	TOTAL IN = 1200000.0000 TOTAL IN = 2000.0000
ex-gwf-fhb.lst_455	STO-SS = 140.9305 STO-SS = 0.2045 STORAGE	STO-SS = 18.0000 STO-SS = 1.8977E-12 STORAGE
ex-gwf-fhb.lst_457	CHD = 2960349.2840 CHD = 5394.4236 CHD	CHD = 1199981.9991 CHD = 2000.0000 CHD
ex-gwf-fhb.lst_459	TOTAL OUT = 2960490.2146 TOTAL OUT = 5394.6281	TOTAL OUT = 1199999.9991 TOTAL OUT = 2000.0000
ex-gwf-fhb.lst_461	IN - OUT = -2.7022E-06 IN - OUT = -1.9372E-10	IN - OUT = 8.9541E-04 IN - OUT = 1.4563E-06
ex-gwf-fhb.lst_463	PERCENT DISCREPANCY = -0.00 PERCENT DISCREPANCY = -0.00	PERCENT DISCREPANCY = 0.00 PERCENT DISCREPANCY = 0.00
ex-gwf-fhb.lst_565	STO-SS = 43.0601 STO-SS = 0.2153 STORAGE	STO-SS = 3.6826E-12 STO-SS = 0.0000 STORAGE
ex-gwf-fhb.lst_566	WEL = 5541858.9256 WEL = 9000.0000 WEL	WEL = 2000000.0000 WEL = 2000.0000 WEL
ex-gwf-fhb.lst_569	TOTAL IN = 5541901.9857 TOTAL IN = 9000.2153	TOTAL IN = 2000000.0000 TOTAL IN = 2000.0000
ex-gwf-fhb.lst_573	STO-SS = 178.0613 STO-SS = 0.0000 STORAGE	STO-SS = 18.0000 STO-SS = 1.0503E-12 STORAGE
ex-gwf-fhb.lst_575	CHD = 5541723.9244 CHD = 9000.2153 CHD	CHD = 1999981.9986 CHD = 2000.0000 CHD
ex-gwf-fhb.lst_577	TOTAL OUT = 5541901.9857 TOTAL OUT = 9000.2153	TOTAL OUT = 1999999.9986 TOTAL OUT = 2000.0000
ex-gwf-fhb.lst_579	IN - OUT = -3.6778E-06 IN - OUT = -4.0709E-09	IN - OUT = 1.4449E-03 IN - OUT = 1.3228E-06
ex-gwf-fhb.lst_581	PERCENT DISCREPANCY = -0.00 PERCENT DISCREPANCY = -0.00	PERCENT DISCREPANCY = 0.00 PERCENT DISCREPANCY = 0.00

Some values are significantly different.

[pya]

These are the excerpts from ex-gwf-csub-p03b.lst:

MF6

  VOLUME BUDGET FOR ENTIRE MODEL AT END OF TIME STEP   22, STRESS PERIOD 353
  ---------------------------------------------------------------------------------------------------

     CUMULATIVE VOLUME      L**3       RATES FOR THIS TIME STEP      L**3/T          PACKAGE NAME    
     ------------------                 ------------------------                     -------------------                       

           TOTAL OUT =           3.7933             TOTAL OUT =       5.7583E-05

            IN - OUT =      -1.4365E-05              IN - OUT =      -4.7444E-10

 PERCENT DISCREPANCY =          -0.00     PERCENT DISCREPANCY =          -0.00

modflowapi

  VOLUME BUDGET FOR ENTIRE MODEL AT END OF TIME STEP   22, STRESS PERIOD 353
  ---------------------------------------------------------------------------------------------------

     CUMULATIVE VOLUME      L**3       RATES FOR THIS TIME STEP      L**3/T          PACKAGE NAME    
     ------------------                 ------------------------                     ----------------

     ...                    

           TOTAL OUT =           2.4318             TOTAL OUT =       6.7702E-05

            IN - OUT =      -6.1185E-07              IN - OUT =      -5.7306E-12

 PERCENT DISCREPANCY =          -0.00     PERCENT DISCREPANCY =          -0.00

The numbers should be identical. For all other model runs, except the four listed above, they are identical.

[jesshe]

I've run into this issue recently and basically get the same results as @pya above.

@jlarsen-usgs looks like that check you ran with the code snippet above is comparing results from the same simulation... ainc and acum should be reading from your "alist" variable not "list"?

[jlarsen-usgs]

@jesshe
good catch, I'll look back into this issue and see what I can track down

[jesshe]

@jlarsen-usgs
Hey wondering if there's any news on this? I think the whole issue can definitely be bypassed if the api is used to handle time variation instead of time series files but I'm wondering if its worth re-configuring my models this way (it its going to get patched relatively soon).

[jlarsen-usgs]

@jesshe
Sorry that work on this has been put to the side for so long. Investigating and trying to patch this is the next priority on my todo list, so I expect to be starting to work on this issue sometime in the next week or so. I'm not exactly sure where the issue is occurring at this point; on the MODFLOW side, python side, or communication between MODFLOW and python via the API.