[Question] Comparison stationary heads for IMOD Python run vs IMOD Batch function

[CasvHWSBrabantseDelta]

Hi Developers,

I am currently trying to import a regional gwf model from imod5 into imod python. I started with writing and running a stationary model and managed to get output. I did a check where I compared the heads from the stationary model in imod Python with heads from the stationary model runned with the RUNFILE batch function. I encountered some differences up to 30 cm on some places in the model, especially close to rivers. I was wondering if I entered exactly the same input settings for both model runs. The imported .prj file is exactly the same, so it has nothing to do with the input files. I also checked the solver settings and could not find what was different. Below I provide you the settings for both model runs. I also uploaded the projectfile and the IMS files from both runs.

Question
Can you help me with figuring out if all input settings for both model runs are exactly the same, or that I missed something? Could it otherwise just be a model uncertainty? I doubt that.

Kind regards,

Cas van Haperen, WS Brabantse Delta

RUNFILE Batch Function

REM Running PRJ file, creating MODFLOW NAM file and writing results to output folder

ECHO FUNCTION=RUNFILE >model3.ini                                                                                                                                                    
ECHO PRJFILE_IN="E:\brabantmodel\AaofWeerijs\OPLEVERING\model\Initial_model_STAT\Turfvaart_ZoW_model_2025_STAT1.prj" >>model3.ini
ECHO SIM_TYPE=3 >>model3.ini
ECHO ISOLVE=1 >>model3.ini
ECHO ISTEADY=0 >>model3.ini
ECHO OUTPUT_FOLDER= "E:\brabantmodel\AaofWeerijs\OPLEVERING\OUTPUT\Turfvaart_ZoW_model_2025_STAT_new">>model3.ini
ECHO NAMFILE_OUT="E:\brabantmodel\AaofWeerijs\OPLEVERING\OUTPUT\Turfvaart_ZoW_model_2025_STAT\Turfvaart_ZoW_model_2025_STAT_new.NAM" >>model3.ini

REM Settings for simulation

ECHO SDATE=20130101000000 >>model3.ini
ECHO EDATE=20240101000000 >>model3.ini

ECHO ISS=0 >>model3.ini                                                                                                                                                                    
ECHO ITT=3 >>model3.ini                                                                                                                                                                        
ECHO IDT=1 >>model3.ini
ECHO NLAY = 43 >>model3.ini
ECHO NSTEP=1 >>model3.ini
ECHO NMULT=1 >>model3.ini                                                                                                                                                                    
ECHO SSYSTEM=1 >>model3.ini
ECHO ISAVEENDDATE=0 >>model3.ini
ECHO ICONCHK=0 >>model3.ini
ECHO IDOUBLE=0 >>model3.ini 
REM ECHO IDEFLAYER=2 >>model3.ini
ECHO SAVESHD=-1 >>model3.ini
ECHO WINDOW=92000.000,368000.000,113000.000,400000.000 >>model3.ini
ECHO CELLSIZE=100 >>model3.ini
ECHO BUFFER=0.0  >>model3.ini
ECHO BUFFERCS = 100 >>model3.ini

REM RUN WITH MODFLOW 6

ECHO MODFLOW6=E:\brabantmodel\AaofWeerijs\IMOD_V5_6\bin\MODFLOW6_v6.4.2.exe >>model3.ini
                                                                                                                                         
"E:\brabantmodel\AaofWeerijs\IMOD_V5_6\iMOD_V5_6.EXE" model3.ini
Pause

IMOD python settings

directory = Path("E:/brabantmodel/AaofWeerijs/OPLEVERING")
model_dir = directory / "model/Initial_model_STAT"
Dbase_dir = directory / "Dbase"
prj_path = model_dir / "Turfvaart_ZoW_model_2025_STAT1.prj"
log_file = model_dir / "build_MODFLOW6_model.log"
logger = get_logger(name="mijn script", log_level=logging.INFO, log_file=log_file) ### verander "mijn script" in de gewenste Log naam. 

mf6_sim = from_imod5_to_mf6_model(prj_path=prj_path,times = pd.DatetimeIndex([pd.Timestamp("2013-01-01")]))
#%%
logger.info("clip model data")
mf6_dir = model_dir / "mf6_STAT"
gwf_model = mf6_sim["imported_model"]

xmin,ymin,xmax,ymax = 92000.000,368000.000,113000.000,400000.000
gwf_model = gwf_model.clip_box(x_min=xmin, y_min=ymin, x_max=xmax, y_max=ymax)
#%%
mf6_sim["imported_model"] = gwf_model

dis = gwf_model["dis"]
idomain = dis["idomain"]
mask = (idomain != 0)
mf6_sim.mask_all_models(idomain)

gwf_model["riv-2riv"].cleanup(dis=gwf_model["dis"])

gwf_model["oc"] = imod.mf6.OutputControl(save_head="last", save_budget="last",)

times_str =["2013-01-01", "2013-01-02","2013-01-03", "2013-01-04"]
mf6_sim.create_time_discretization(additional_times=times_str)

outer_dvclose = 0.001
outer_maximum = 500
inner_maximum = 25
inner_dvclose = 0.001
inner_rclose = 0.1
linear_acceleration = "cg"
mf6_sim["ims"] = imod.mf6.Solution(modelnames=["imported_model"],outer_dvclose=outer_dvclose, outer_maximum=outer_maximum,
                                    inner_maximum=inner_maximum, inner_dvclose=inner_dvclose, inner_rclose=inner_rclose,
                                    linear_acceleration=linear_acceleration)

if mf6_dir.exists():
    shutil.rmtree(mf6_dir)

mf6_sim.write(directory = mf6_dir, validate=True)

mf6_exe = r"E:\brabantmodel\AaofWeerijs\IMOD_V5_6\bin\MODFLOW6_v6.4.2.exe"
mf6_sim.run(mf6path=mf6_exe)

Projectfile_STAT_model.txt

MFSIM_IMS6_batchfunction.txt

MFSIM_IMS_IMOD_PYTHON.txt

[JoerivanEngelen]

Hi Cas,

Thanks for posting this issue. From what I see, your settings are fine here. There are two moments where differences can manifest between iMOD5 and iMOD Python input:

• Upon importing the model river systems
• When the cleanup of river sys 2 is called (probably because you ran into a ValidationError)

I think I need to debug through your model to see what is the cause of the differences. Could you share a small copy of the model with us around an area where the differences manifest a lot? You can create this with the iMOD5 MODELCOPY. Please keep clip such that the grids are small, e.g. a clipout of 30 by 30 cells. The resulting folder might be too large to share here as a zip. If so, you can email the model to imod.support@deltares.nl with a link to this Github Issue.

Good luck!
Joeri

[CasvHWSBrabantseDelta]

SUBMODEL_Aa_of_Weerijs_STAT.zip

Hi Joeri, I uploaded the submodel as you asked, I only copied the given .prj file into a 30x30 grid. Do you also need a .RUN file in this folder, or can you make this yourself?

Thanks for trying to help me!

Cas

[JoerivanEngelen]

Hi Cas, thanks for sharing the model. Final thing that I missed to ask you last time: in the script snippet you share, you call a function that is not part of iMOD Python: from_imod5_to_mf6_model(prj_path=prj_path,times = pd.DatetimeIndex([pd.Timestamp("2013-01-01")]))

I need to see what that function does to better pinpoint differences, as I need to know what it does internally. Can you share that function with me as well? If it refers to other custom functions please also share these.

[CasvHWSBrabantseDelta]

Hi Joeri, Here is the function I wrote, it just uses your functions, nothing special extra.

from pathlib import Path
import pandas as pd
import logging
from imod.mf6 import Modflow6Simulation
from imod import prj

logger = logging.getLogger(name)

def from_imod5_to_mf6_model(prj_path: Path, times: list[pd.Timestamp])->Modflow6Simulation:
"""Reads a .prj file and creates MODFLOW6 input.

Args:
    prj_path (Path): The input .prj file containing the model input file path locations to the imod5 DBASE.
    times (list[pd.Timestamp]): list of timesteps for the model run.

Returns:
    Modflow6Simulation: MODFLOW6 simulation.
"""
logger.info(f"Reading imod5 data from {prj_path}")
imod5_data, period_data = prj.open_projectfile_data(prj_path)
logger.info("Create Modflow6Simulation")
return Modflow6Simulation.from_imod5_data(imod5_data, period_data, times)

If you have other questions, feel free to ask!

Thanks for helping!