Modflow-MetaSWAP coupler mappings should use grid-agnostic wells or well adapter

[JoerivanEngelen]

The msw.CouplerMapping , msw.Sprinkling, and couplers.MetaMod (in primod) currently use the deprecated mf6.WellDisStructured object. We should start supporting the new architecture.

This will also require updates to:

• primod because coupler logic is now placed there
• The example in the iMOD Documentation

We might be able to use the Mf6Well adapter directly here:

• Less work to update existing _render method
• No mf6 data needs to be provided to render method. With grid agnostic packages, this would be (idomain, bottom, top, k)
• However, this makes it harder to clip and regrid this package object this would require the regrid_like to update the indices

The alternative is providing the grid agnostic Well package:

• Easy to clip & regrid
• This requires extra information to the _render method, all contained in the GroundwaterFlowModel object.

The third option is creating a GridWell package: A grid of well rates, comparable to other boundary conditions.

• Behaves similair to other packages
• Quite a change in public API: Would users like this?

Also: Since MODFLOW6 models have an idomain correction after regridding (some NPF cells may become inactive at edges, where idomain cells were not inactive), we always have to regrid the GroundwaterFlowModel first, before regridding MetaSWAP model.

[JoerivanEngelen]

Refinement

• CouplerMapping takes a StructuredDiscretization object. It has to be initialized again each time the Modflow model is clipped, regridded, or masked. So for this object using an Mf6Well object is sufficient.
• Sprinkling is more complicated, crude sketch:

if type(well) == GridAgnosticWell
    sprink = Sprinkling(well)
    sprink.regrid_like(new_idomain)

    gwf_regridded = gwf_model.regrid_like(new_idomain)
    idomain_regridded = gwf_regridded["dis"]["idomain"]
    top_regridded = gwf_regridded["dis"]["top"]
    bottom_regridded = gwf_regridded["dis"]["bottom"]
    k_regridded = gwf_regridded["npf"]["k"]

    sprink._render(file, index, svat, idomain_regridded, top_regridded, bottom_regridded, k_regridded)
    # In _render, self.well.to_mf6_pkg is called

elif type(well) == Mf6Wel
    mf6_wel = well.to_mf6_pkg(idomain, top, bottom, k)
    sprink = Sprinkling(mf6_wel)

    gwf_regridded = gwf_model.regrid_like(new_idomain)
    idomain_regridded = gwf_regridded["dis"]["idomain"]
    top_regridded = gwf_regridded["dis"]["top"]
    bottom_regridded = gwf_regridded["dis"]["bottom"]
    k_regridded = gwf_regridded["npf"]["k"]
    mf6_wel_regrid = well.to_mf6_pkg(idomain_regridded, top_regridded, bottom_regridded, k_regridded)
    sprink.regrid_like(new_idomain)
    sprink.update_wel(mf6_wel_regrid)

    sprink._render(file, index, svat)

[JoerivanEngelen]

For this specific issue, running coupled simulations is not a hard requirement, as the rendered textfiles need to remain the same. However, for the future solving https://github.com/Deltares/imod-python/issues/1154 would be useful

[JoerivanEngelen]

Doing some more thinking on this: It is probably better to refactor the MetaSWAP code base in such a way, that the MODFLOW objects are only necessary upon writing, instead of already upon initialization. For example for the Sprinkling package:

Current state:

class Sprinkling:
    def __init__(
        self,
        max_abstraction_groundwater: xr.DataArray,
        max_abstraction_surfacewater: xr.DataArray,
        modflow_wel: WellDisStructured
     ):
        self.well = modflow_wel
...

    def write(
        self,
        directory: Union[str, Path],
        index: np.ndarray,
        svat: xr.DataArray,
    ):
...

Initial idea:

class Sprinkling:
    def __init__(
        self,
        max_abstraction_groundwater: xr.DataArray,
        max_abstraction_surfacewater: xr.DataArray,
        modflow_wel: Well
     ):
        self.well = modflow_wel
...

    def write(
        self,
        directory: Union[str, Path],
        index: np.ndarray,
        svat: xr.DataArray,
    ):
...

Proposed:

class Sprinkling:
    def __init__(
        self,
        max_abstraction_groundwater: xr.DataArray,
        max_abstraction_surfacewater: xr.DataArray,
        mf6_wellname: str,
     ):
        self.mf6_wellname = mf6_wellname
...

    def write(
        self,
        directory: Union[str, Path],
        index: np.ndarray,
        svat: xr.DataArray,
        modflow_wel: Mf6Wel
    ):
...

This propsed approach has the following advantages over the initial idea:

1. Modflow objects are only used when really necessary. This has the advantage that we can directly use the more low-level Mf6Wel object instead of the grid agnostic Well and LayeredWell packages, as the latter still have to be assigned to cells. We avoid having to conduct this twice.
2. Mutations of Modflow data, for example regrid_like, do not have to be done twice. In the initial idea, regridding would have to be called once for the MetaSWAP model and once for the Modflow6Simulation, as we do not check wether the Modflow package assigned to the two different models is a copy or the same package for each model. The proposed approach also means no calls to Modflow6Package.regrid_like are done outside the mf6 module, reducing clutter a bit.

The same approach can be taken for the CouplerMapping object in iMOD Python, and the NodeSvatMapping, RechargeSvatMapping, WellSvatMapping objects in primod.

The MetaMod.write method would be the place fetch the Modflow packages and pass them on through to MetaSwapModel.write

[JoerivanEngelen]

The proposed approach means quite a serious refactoring exercise. Best approach here is probably first trying to make this:

class Sprinkling:
    def __init__(
        self,
        max_abstraction_groundwater: xr.DataArray,
        max_abstraction_surfacewater: xr.DataArray,
        modflow_wel: Mf6Wel
     ):
        self.well = modflow_wel
...

    def write(
        self,
        directory: Union[str, Path],
        index: np.ndarray,
        svat: xr.DataArray,
    ):
...

With this we can test if no changes in derived mappings occur when coupling to Mf6Wel, instead of WellDisStructured. It means the Sprinkling package still cannot be regridded, but this is acceptable for now, as this is already the present state.

I will create a separate issue with the proposed refactoring.