Inconsistent river flow values from SIMVALS("RIV") vs .lst/.cbc outputs

[manupy-hub]

Hi, I’m learning to use the MODFLOW‑6 API and trying to integrate it with other platforms, but I’ve run into some inconsistencies—probably because I’m doing something wrong. I’m extracting river flows at each time step with:

Q_out = gwf.get_value(
    gwf.get_var_address("SIMVALS", gwf_name, "RIV")
)

Although this appears to work, the values I get from SIMVALS("RIV") don’t match those printed in the .lst and .cbc files at the end of the simulation (which agree with each other).

Why does SIMVALS("RIV") return different values? Am I calling it at the wrong point in the time step, or is there another transformation applied to the .lst/.cbc outputs that I’m missing? Any pointers would be greatly appreciated.

Thank you for your help!

`while ctime < etime:
gwf.update()
river_nodes =gwf.get_value(gwf.get_var_address("nodelist", gwf_name, "RIV"))
river_nodes

# number of iterations allowed for the linear solver
max_iter = gwf.get_value(gwf.get_var_address("MXITER", "SLN_1"))
# dt       = gwf.get_time_step()

# # prepare this time step
# gwf.prepare_time_step(dt)
gwf.prepare_solve(1)

# fetch one‐based stress period and time step indices
stress_period = int(gwf.get_value(gwf.get_var_address("KPER", "TDIS"))[0])
time_step     = int(gwf.get_value(gwf.get_var_address("KSTP", "TDIS"))[0])
# if we're in the user‐specified SP, zero out Riv
if stress_period == sp_to_update:
    # get the correct CHD array

    river_nodes_curr =gwf.get_value(gwf.get_var_address("nodelist", gwf_name, "RIV"))
    addr_bound = gwf.get_var_address("BOUND", gwf_name, "RIV")
    bound=gwf.get_value(addr_bound)

    
    ###############assign Conduct 0 to bound in SP 3
    bound[:,1]=0
    
    bound=np.ascontiguousarray(bound, dtype=np.float64)
    
    
    gwf.set_value(addr_bound, bound)

    

    kiter = 0
    # print(hcof)
    while kiter < max_iter:
        # print(kiter)
        if gwf.solve(1):
            print(f"SP {stress_period}, TS {time_step} converged in {kiter} iterationssssssss")
            
            break
        kiter += 1
    gwf.finalize_solve(1)
    print('Q_out: ',np.sum(gwf.get_value(gwf.get_var_address("SIMVALS", gwf_name, "RIV"))))




else:

    kiter = 0
    while kiter < max_iter:
        if gwf.solve(1):
            print(f"SP {stress_period}, TS {time_step} converged in {kiter} iterations")
            
            break
        kiter += 1
    gwf.finalize_solve(1)
    print('Q_out: ',np.sum(gwf.get_value(gwf.get_var_address("SIMVALS", gwf_name, "RIV"))))

ctime = gwf.get_current_time()

gwf.finalize()`

Outputs from the coode (simvals)
SP 1, TS 3 converged in 0 iterations
Q_out: -248.10936681052772

SP 2, TS 3 converged in 8 iterations
Q_out: -1050.4023165637536

SP 3, TS 3 converged in 2 iterations
Q_out: -67.97356541245426

SP 4, TS 3 converged in 6 iterations
Q_out: -1040.2878996694453

SP 5, TS 3 converged in 17 iterations
Q_out: -986.5099712653791

outputs in cbc
Q_out (SP 1) : -248.08753756890292
Q_out (SP 2) : -1001.3297650143686
Q_out (SP 3) : 0.0
Q_out (SP 4) : -992.6270730347387
Q_out (SP 5) : -923.2945007972521

outputs from lst

Sp1

SP3

[mjr-deltares]

Hi @manupy-hub , would you mind sharing the full example so we can try to reproduce? Maybe by adding it as a file to this issue.