Loading an existing model from MODFLOW 6 input files

[phaehnel]

Dear developers,

since I have several MODFLOW 6 models which were not created with imod-python, I was wondering if there is a way to load an existing MODFLOW 6 model with imod-python from the MODFLOW 6 input files directly. I couldn't find anything in the examples (https://deltares.github.io/imod-python/examples/index.html) or the API documentation (https://deltares.github.io/imod-python/api/index.html) regarding this.

So far I only found the option to load an existing model with the sim.from_file() method that loads the existing model via a TOML file and netCDF files that represent the different MODFLOW packages' input files. This works fine if I have a model create with imod-python and used sim.dump() beforehand to get the data required by sim.from_file().

There are two scenarios where I think this could be useful:

• The abovementioned model that was not created with imod-python. I might want to switch to imod for performance reasons (e.g., model splitting is more efficient with imod-python than with flopy; making a large-scale steady-state model transient with a very large amount of stress periods) but not recreate the entire model setup from scratch.
• I created the model with imod-python but it is quite large and input files already require a lot of disk space. I don't want to use more disk space to save a dump version to be able to load the model later on again without re-running the script that created the model. I see that the netCDF files use way less space than ASCII-based MODFLOW 6 input files but I can imagine highly-dynamic transient boundary conditions where this approach could become very disk-space demanding.

Maybe I also have not quite fully understood how to work with the dump and TOML option or I have overseen some other functionality. However, I would be interested if you see any way in approaching this, e.g., the two scenarios mentioned above?

Cheers,
Patrick

[Huite]

Hi @phaehnel,

Models not written by imod-python

Unsurprisingly, we run into similar issues when it comes to models made with alternative tools! I have posted a basic "Flopy-converter" here which we needed for a specific project:
#1678

This uses Flopy to read the files (so at least we can deal with Python objects). An alternative approach can be found on this (very, very stale) PR:
#688
(Incidentally, I needed this because the model was written with Fortran and could not be read with Flopy due to the use of repetition factors in the grid input among things -- so something like 365*0 in the ibound expands to an array with 365 zeros; this is only supported by the Fortran read intrinsic procedure, so MODFLOW 6 will accept it but not Flopy ...)

You might wonder why we don't "just" support reading MODFLOW 6 input files. The reason is that, unfortunately, MODFLOW 6 has a bewildering variety in its inputs. To name a few aspects:

• Inline data vs data in external files
• In external files: ascii and binary input
• fully irregular time series support (through time series files)
• constant / layered directives
• grid input versus list-based input

We based imod-python on xarray and grid representations, because these are more transparent and easier to work with. Unfortunately, this means that there's no room for time series with alternatively spaced time indexes within a single package, or list-based inputs that have multiple boundaries for the same cell (since the grid representation can only hold a single value for a cell); as a consequence, list-inputs have to be expanded into separate package instances when imported into imod-python.

For this reason, we've been hesitant to merge efforts into a release like a Flopy converter or the Modflow 6 input parsing. We could a read a certain subset of MODFLOW 6 inputs relatively easily, but e.g. list-based input with multiple boundaries in a single cell is extremely common and already requires a somewhat ugly workaround. And before you know, it's starts getting pretty complex with a lot of potential for bugs.(Worth noting: MODFLOW 6 is adding support for netCDF input with fully gridded data which matches our data model exactly.)

The links above probably have useful code, but we haven't found the time to test them properly. So my recommendation would be: feel free to use it, but be aware that you need to check its correctness carefully and that you likely need to adapt it for your specific model (again, because MODFLOW 6 input has so many flavors).

Disk space

Casting everything a grid data indeed increases storage sizes, especially if the input is relatively sparse (e.g. a river package with boundaries on only 1% of the cells). Fortunately, this has a relatively simple workaround: enable compression on the netCDF writing. This is currently pretty tedious though, I guess: it requires setting the encoding on every package dataset (via ds.encoding["zlib"] = True if memory serves me). It might make sense to add a compression argument to something like that dump method. I'll open an issue.

Writing to Zarr instead of netCDF will automatically enable (very performant) compression. Incidentally, I've added an engine argument in this PR #1686.
(If you're looking at model splitting, that PR is relelvant too...)
Zarr has the advantage of being generally quite a bit faster than netCDF (up to something like a factor of four), but the zarr spec is still a work in progress and and so it's still unsuitable for long term storage.

[phaehnel]

Dear @Huite, thank you for the detailed reply and the valubale insights! That helps me a lot!

Models not written in imod-python:
I've seen the netCDF capabilities in the MODFLOW 6 Extended version, so I should test this, I reckon. Sounds like the netCDF approach will provide the missing link in the future, so that the entire model can be defined and run based on netCDF input?

Features in general:
And one follow-up question: If your data model does not match the variety of options possible with MODFLOW 6, does that mean that there are certain things that I can absolutely not represent with imod-python that could however generally be done with MODFLOW 6? You mentioned some workarounds are possible, like defining different packages to stay consistent with the data model. The overarching question I guess is: If I decide to use imod-python, do I limit myself to a subset of the vast variety of functionalities that MODFLOW 6 provides? I know that imod-python is still under development and there are probably a lot of features still to come.

I ask this because I am part of a team that develops a larg-scale groundwater model and we need to decide on the tools we will utilize. Imod-python is very interesting due to the very good performance with large inputs and creating large models in general. However, we are currently a bit lost regarding the availability of features we might need in the future (e.g., streamflow routing (I've seen its under development), water mover package, or particle tracking are potentially of interest for us in the future) and if we'd rather compromise speed in creating the model for feature availability.

Disk space:
The zlib option seems like a good approach, I've used that in the past to save extremely large datasets with netCDF and the compression was quite impressive. And yes, if I remember correctly, a compression argument had to be defined for every created dataset as one might not want to use the compression on all datasets.

[Huite]

Sounds like the netCDF approach will provide the missing link in the future, so that the entire model can be defined and run based on netCDF input?

Yes. Also worth noting: the Flopy developers are looking into rewriting the MODFLOW 6 functionality, basing it around xarray as this makes netCDF input very straightforward, and provides the opportunity to benefit of some of the lessons from imod-python (support larger than memory models, fix performance). At the same time, unlike imod-python, Flopy has to keep feature parity with MODFLOW 6, because it plays an essential role in the MODFLOW 6 testing framework; such a rewrite is obviously a substantial task.

And one follow-up question: If your data model does not match the variety of options possible with MODFLOW 6, does that mean that there are certain things that I can absolutely not represent with imod-python that could however generally be done with MODFLOW 6? You mentioned some workarounds are possible, like defining different packages to stay consistent with the data model. The overarching question I guess is: If I decide to use imod-python, do I limit myself to a subset of the vast variety of functionalities that MODFLOW 6 provides? I know that imod-python is still under development and there are probably a lot of features still to come.

This is an excellent question. I think the answer is no: everything MODFLOW 6 offers should be possible -- although perhaps in a slightly opinionated way. I will give some concrete examples which hopefully illustrates our philosopy.

• Time series: In my view, MODFLOW 6 time series are mostly a workaround about bad tooling. E.g. if you're editing models by hand, it's very convenient to be able to just adjust a text table (possibly even with Excel). I would also argue that time series are error prone: they linearly interpolate, so it's less obvious what is being fed into the model. (It degrades performance on the MODFLOW 6 Fortran side too, because the input parser never knows what to expect and cannot specialize.) What's the workaround? Xarray makes it trivial to broadcast a time axis, with strong date time handling provided by e.g. a pandas TimeIndex.
• Observations: we have not prioritized these at all, because post-processing tooling easily provides the same results. Just load the heads grab the values at the required locations.
• List-based input: the primary innovation of the original (Fortran) iMOD is that it's grid focused. This increases transparency and ease of data inspection (just call .plot()), and makes the manipulating the data much easier because the location of a boundary element is directly encoded by the position in the array/grid. Workaround: in our Netherlands models, this means that we separate our river packages according to primary, secondary, tertiary water bodies (and have done so for decades).

However, we are currently a bit lost regarding the availability of features we might need in the future (e.g., streamflow routing (I've seen its under development), water mover package, or particle tracking are potentially of interest for us in the future).

That's understandable. Flopy supports all of MODFLOW 6, but expects you to provide the data largely as MODFLOW 6 wants to see it: bare (numpy) arrays, or tables in the form of recarrays or pandas DataFrames. Part of what makes the imod-python features lag behind is that we try to provide data structures that are little nicer. Sometimes these data structures require developing a separate package like xugrid; the nice thing is that it can then be used beyond MODFLOW.

In general, the nicer data structures are:

• Xarray for DIS input
• Xugrid with 2D topologies for DISV input
• Xugrid with 1D topologies for SFR networks
• X, Y data for wells

A key element is that these all remain spatially associated (they are not stripped to cell indices). As a consequence, they can all be visualized relatively directly, we can apply spatial operations like clipping and regridding to very quickly build submodels, do local refinements, etc.

Another example: for SFR, flopy requires you to associate each reach with the appropriate cell index. imod-python can do this automatically, because the data structure (with a Ugrid1d network) is spatial and can be intersected with the DIS grid to do automatic assocation. With flopy, you'd have to write your own intersection logic, number the elements, then feed them to the Flopy SFR class.

But: we can always skip the "nicety" and just write the tables directly like MODFLOW 6 wants. You can subclass the imod.mf6.Package and add your own input which could basically take a pandas DataFrame. The number of methods required for a simple package isn't that large, so something like a mover package can be written very quickly. It might not support larger-than-memory datasets or support clipping, but it'll write the input. The water moved could be done in a very quick way, but then it's up to the user to make sure everything is connected in the right way. (It's an interesting question what a "user-friendly" and transparent water mover looks like: a GeoDataFrame with lines? Almost like a set of physical pipes?)

Working with imod-python and Flopy

Even if you find you need to write models with Flopy, you will likely develop the workflows with:

• pandas
• geopandas
• xarray
• xugrid, e.g. for the regridding

In practice, this means you can remain somewhat imod/FloPy agnostic. Ideally, your intermediate datasets could be almost directly ingested by the imod-python classes. To use FloPy, you would only need to strip the spatial context and pass the underlying NumPy arrays. This way, you can use imod-python in your “happy path,” but fall back to FloPy when needed -- without extra effort compared to starting with FloPy alone.

Bottom line

I would definitely send us an e-mail and schedule a Teams meeting! My colleague @JoerivanEngelen is currently on birth leave, but he'll be back shortly and he can definitely help clear things up. The iMOD User Day might also be worthwhile (if you didn't sign up already): https://oss.deltares.nl/web/imod/-/international-imod-user-day-2024-duplicate-2-2

In terms of MODFLOW 6 capabilities, there are unlikely to be immediate limitations. Building a large groundwater model from scratch always brings challenges, and an iterative approach is best: start simple (e.g. use Rivers before SFR, avoid convertible cells initially), validate performance, and only then add complexity.

If you communicate in advance which additional packages you will require, we can take this into account in our prioritization. By the time you are happy with earlier iterations, I wager that the odds are good that something like Particle Tracking will be up and running. And maybe this isn't clear: but if you start a project or collaboration with us, you can of course directly influence our development priorities. When we can align your needs with active development (especially in project form with some funding), we can accelerate work on features like the Water Mover or Particle Tracking packages, or even help with some of the workflows.

[phaehnel]

Yes. Also worth noting: the Flopy developers are looking into rewriting the MODFLOW 6 functionality, basing it around xarray as this makes netCDF input very straightforward

Cool, this sounds awesome! I really like working with xarrays and xugrids. They make it so much easier to handle spatial data than plain numpy arrays.

This is an excellent question. I think the answer is no: everything MODFLOW 6 offers should be possible

Thanks for the examples. I noticed, when switching from MODFLOW 2005 to MODFLOW 6m that there are often several ways now to achieve the same thing. Probably convinient for the user but difficult to deal with when developing pre- and post-processing tools.

But: we can always skip the "nicety" and just write the tables directly like MODFLOW 6 wants. You can subclass the imod.mf6.Package and add your own input which could basically take a pandas DataFrame.

That's interesting. I'll have to look into that at some point. Thanks for pointing it out.

To use FloPy, you would only need to strip the spatial context and pass the underlying NumPy arrays.

Yep. Exactly what I've been doing so far to compare flopy and imod-python. And it's also just sensible to work with xarray/xugrid anyways.

In terms of MODFLOW 6 capabilities, there are unlikely to be immediate limitations. Building a large groundwater model from scratch always brings challenges, and an iterative approach is best: start simple (e.g. use Rivers before SFR, avoid convertible cells initially), validate performance, and only then add complexity.

Yes, we are at the very beginning and we just try to think a bit ahead what tools might be relevant for the model beyond the first version. Just to omit going down a path were we can add more complex capabilities only with a lot of effort and rewriting of routines or even switching tools/software.

I would definitely send us an e-mail and schedule a Teams meeting! My colleague @JoerivanEngelen is currently on birth leave, but he'll be back shortly and he can definitely help clear things up. The iMOD User Day might also be worthwhile (if you didn't sign up already): https://oss.deltares.nl/web/imod/-/international-imod-user-day-2024-duplicate-2-2

If you communicate in advance which additional packages you will require, we can take this into account in our prioritization. By the time you are happy with earlier iterations, I wager that the odds are good that something like Particle Tracking will be up and running. And maybe this isn't clear: but if you start a project or collaboration with us, you can of course directly influence our development priorities. When we can align your needs with active development (especially in project form with some funding), we can accelerate work on features like the Water Mover or Particle Tracking packages, or even help with some of the workflows.

Thank you for offering a meeting! This would be really great and helpful for us. I'll forward our conversation to my colleagues and we will contact you then soon. We have been in contact from time to time with @JoerivanEngelen and also visited the iMOD User Day in Hamburg earlier this year. Some of my colleagues and I are also going to attend the User Day in Delft in November. But in general, a meeting to get to know each other and talk about what we are working on, about imod and future developments, or the NHI would be great. Our project is not as comprehensive and complex as the NHI at this stage but we are dealing with the same spatial scale and the challanges that come along with that.

[Huite]

Oh and I highly recommend writing binary input for MODFLOW 6! imod-python does so by default as you probably noticed, but it's substantially faster in writing.
For Fortran too, binary input is substantially faster, so the MODFLOW 6 model will initialize much faster.
I recommend writing ascii only for input checking if things are going wrong.

[phaehnel]

Yes, I am aware of that and noticed the speedup and decrease in storage requirement. As you mentioned, in model development, I like to be able to take a look at the input files in case MODFLOW or the results show strange behaviour. For a tested model, the binary approach is very handy indeed!