How to create SFR input for a MF6 DISV grid in Flopy?

[nikobenho]

Hi,

I'm building a MF6 octree

I'm creating a MF6 disv-grid model in Flopy and want to represent a creek using the SFR package. Does anyone have experience adding SFR input to a MF6 DISV grid from a shapefile (or Shapely LineString)?

I've tried using the SFRmaker package. But unfortunately, the support for unstructured grids aren't fully implemented, and I haven't managed to find a suitable workaround.

Does anyone know of a way to create SFR input from a shapefile using Python & Flopy?

[Richard B. Winston]

This isn't exactly what you asked for but it could be a way forward for you.

1. Use ModelMuse to create your DISV grid.

2. Activate the SFR package in the Model|MODFLOW Packages and Programs dialog box.

3. Import the SFR data from a Shapefile using File|Import|Shapefile.

4. Export the model input files.

5. Import the model with Flopy and generate the rest of the model input using Flopy.

Incidentally, octtree refinement  would typically produce a DISU grid rather than a DISV grid.

[nikobenho]

Thank you for replying Richard,

It looks like this is the strategy I will choose. At first I wanted to keep my entire workflow in Python (reduce cognitive load), but grid setup is much more conveniently handled in ModelMuse (at least in my experience).

For anyone taking the same approach as i did (create SFR using ModelMuse, then import the model into Flopy) and encounters the problem where streambed bottom elevations increase in the downward direction, the following Python code will smoothen out the streambed elevation (similar to SFRmaker):

#Load SFR-package

sfr = m.get_package('sfr')

smoothed_sfr_elevation = []
current_min = 0
for count, element in enumerate(sfr.packagedata.array.rtp):
    if count == 0:
        smoothed_sfr_elevation.append(element)
        current_min = element
    else:
        if element > current_min:
            new_min = current_min - 0.0001
            smoothed_sfr_elevation.append(new_min)
            current_min = new_min
        else:
            smoothed_sfr_elevation.append(element)
            current_min = element

#Replace erroneous streambed elevations with smoothed streambed elevations

sfr.packagedata.array.rtp = np.array(smoothed_sfr_elevation)

I hope this can be of use to anyone in a similar situation.

/Nikolas