handling of 'None' reaches mf6 SFR package

[aleaf]

tl;dr: flopy won't load or (correctly) write an SFR package like this one (can provide full working example if needed):

BEGIN options
END options

BEGIN packagedata
  1  1 1 10      35.50000000       1.76000000       0.44800000     293.20000000       1.00000000       1.00000000  0.037  1       1.00000000  0
  2  NONE  35.499134 1.7631012109518447 0.44803908 293.21634 1.0 1.0 0.03700000047683716 1 1.0 0
  3  NONE  35.149265 3.28178281913793 0.001 293.21634 1.0 1.0 0.03700000047683716 1 1.0 0
  4  1 1 9      35.50000000       1.76000000       0.44800000     293.20000000       1.00000000       1.00000000  0.037  1       1.00000000  0
END packagedata

BEGIN connectiondata
  1  2
  2  3
  3  4
END connectiondata

In the MODFLOW-6 SFR package input, 'NONE' must be specified for SFR reaches that don't have a connection to the groundwater system (for example, those located in cells with IDOMAIN < 1). If I make this model:

sim_name = 'testsim'
model_name = 'testmodel'
exe_name = 'mf6'
out_dir = 'tmp'

# set up simulation
tdis_name = '{}.tdis'.format(sim_name)
sim = mf6.MFSimulation(sim_name=sim_name,
                   version='mf6', exe_name=exe_name,
                   sim_ws=out_dir)
tdis_rc = [(6.0, 2, 1.0), (6.0, 3, 1.0)]
tdis = mf6.ModflowTdis(sim, time_units='DAYS', nper=2,
                          perioddata=tdis_rc)
ims_package = mf6.ModflowIms(sim, print_option='ALL',
                               complexity='SIMPLE', outer_hclose=0.00001,
                               outer_maximum=50, under_relaxation='NONE',
                               inner_maximum=30,
                               inner_hclose=0.00001,
                               linear_acceleration='CG',
                               preconditioner_levels=7,
                               preconditioner_drop_tolerance=0.01,
                               number_orthogonalizations=2)
sim.register_ims_package(ims_package, [model_name])
# create model instance
model = mf6.ModflowGwf(sim, modelname=model_name,
                         model_nam_file='{}.nam'.format(model_name))

dis = mf6.ModflowGwfdis(model, length_units='FEET', nlay=1,
                                     nrow=1, ncol=10, delr=500.0,
                                     delc=500.0,
                                     top=100.0, botm=50.0,
                                     filename='{}.dis'.format(model_name))
sfr = mf6.ModflowGwfsfr(model, packagedata=[
    (0, (0, 0, 9),  35.5, 1.76, 0.448, 293.2, 1.0, 1.0, 0.037, 1, 1.0, 0),
    (1, None,  35.5, 1.76, 0.448, 293.2, 1.0, 1.0, 0.037, 1, 1.0, 0),
    (2, None,  35.5, 1.76, 0.448, 293.2, 1.0, 1.0, 0.037, 1, 1.0, 0),
    (3, (0, 0, 8),  35.5, 1.76, 0.448, 293.2, 1.0, 1.0, 0.037, 1, 1.0, 0)],
                        connectiondata=[(0, 1),
                                        (1, 2),
                                        (2, 3),
                                       ]
)
sim.write_simulation()

it writes ok, but the SFR package looks like this:

BEGIN options
END options

BEGIN packagedata
  1  1 1 10      35.50000000       1.76000000       0.44800000     293.20000000       1.00000000       1.00000000  0.037  1       1.00000000  0
  2
  3
  4  1 1 9      35.50000000       1.76000000       0.44800000     293.20000000       1.00000000       1.00000000  0.037  1       1.00000000  0
END packagedata

BEGIN connectiondata
  1  2
  2  3
  3  4
END connectiondata

I'm not sure if empty reaches like 2 and 3 above are valid MODFLOW-6 input, but I can't load the above file with mf6.MFSimulation.load('testsim', sim_ws='tmp')

Changing reaches 2 and 3 to

  2  NONE  35.499134 1.7631012109518447 0.44803908 293.21634 1.0 1.0 0.03700000047683716 1 1.0 0
  3  NONE  35.149265 3.28178281913793 0.001 293.21634 1.0 1.0 0.03700000047683716 1 1.0 0

is valid for MODFLOW6, but when I try loading with mf6.MFSimulation.load('testsim', sim_ws='tmp'), I get this error

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-64-7283ca140f7b> in <module>
----> 1 sim2 = mf6.MFSimulation.load('testsim', sim_ws='tmp')

~/Documents/GitHub/flopy/flopy/mf6/modflow/mfsimulation.py in load(cls, sim_name, version, exe_name, sim_ws, strict, verbosity_level, load_only, verify_data)
    692             if verbosity_level.value >= VerbosityLevel.normal.value:
    693                 print('  loading model {}...'.format(item[0].lower()))
--> 694             instance._models[item[2]] = model_obj.load(
    695                 instance,
    696                 instance.structure.model_struct_objs[item[0].lower()], item[2],

~/Documents/GitHub/flopy/flopy/mf6/modflow/mfgwf.py in load(cls, simulation, structure, modelname, model_nam_file, version, exe_name, strict, model_rel_path, load_only)
     99              exe_name='mf6.exe', strict=True, model_rel_path='.',
    100              load_only=None):
--> 101         return mfmodel.MFModel.load_base(simulation, structure, modelname,
    102                                          model_nam_file, 'gwf', version,
    103                                          exe_name, strict, model_rel_path,

~/Documents/GitHub/flopy/flopy/mf6/mfmodel.py in load_base(cls, simulation, structure, modelname, model_nam_file, mtype, version, exe_name, strict, model_rel_path, load_only)
    626                     print('    loading package {}...'.format(ftype))
    627                 # load package
--> 628                 instance.load_package(ftype, fname, pname, strict, None)
    629 
    630         # load referenced packages

~/Documents/GitHub/flopy/flopy/mf6/mfmodel.py in load_package(self, ftype, fname, pname, strict, ref_path, dict_package_name, parent_package)
   1116                               parent_file=parent_package)
   1117         try:
-> 1118             package.load(strict)
   1119         except ReadAsArraysException:
   1120             #  create ReadAsArrays package and load it instead

~/Documents/GitHub/flopy/flopy/mf6/mfpackage.py in load(self, strict)
   1530 
   1531         try:
-> 1532             self._load_blocks(fd_input_file, strict)
   1533         except ReadAsArraysException as err:
   1534             fd_input_file.close()

~/Documents/GitHub/flopy/flopy/mf6/mfpackage.py in _load_blocks(self, fd_input_file, strict, max_blocks)
   1652                                 self.blocks[block_key].structure.name))
   1653 
-> 1654                         self.blocks[block_key].load(block_header_info,
   1655                                                     fd_input_file, strict)
   1656                         self._simulation_data.mfdata[self.blocks[block_key].

~/Documents/GitHub/flopy/flopy/mf6/mfpackage.py in load(self, block_header, fd, strict)
    639                         print('        loading data {}..'
    640                               '.'.format(dataset.structure.name))
--> 641                     next_line = dataset.load(line, fd_block,
    642                                              self.block_headers[-1],
    643                                              initial_comment)

~/Documents/GitHub/flopy/flopy/mf6/data/mfdatalist.py in load(self, first_line, file_handle, block_header, pre_data_comments)
    697                                         self._simulation_data, self._path,
    698                                         self._current_key)
--> 699         return file_access.load_from_package(
    700             first_line, file_handle, self._get_storage_obj(), pre_data_comments)
    701 

~/Documents/GitHub/flopy/flopy/mf6/data/mffileaccess.py in load_from_package(self, first_line, file_handle, storage, pre_data_comments)
    802         else:
    803             have_newrec_line, newrec_line, self._data_line =\
--> 804                 self.read_list_data_from_file(file_handle, storage,
    805                                               self._current_key,
    806                                               current_line, self._data_line)

~/Documents/GitHub/flopy/flopy/mf6/data/mffileaccess.py in read_list_data_from_file(self, file_handle, storage, current_key, current_line, data_line, store_internal)
    976                                     # is a cellid
    977                                     try:
--> 978                                         cellid_tuple += \
    979                                             (int(arr_line[sub_entry[0]]) - 1,)
    980                                     except:

ValueError: invalid literal for int() with base 10: 'NONE'

[spaulins-usgs]

@aleaf, confirmed that this is a problem. You can find a quick fix that is working for me on this fork:

https://github.com/spaulins-usgs/flopy

I will revisit this problem next week and see if I can come up with a more elegant solution.

[aleaf]

thanks @spaulins-usgs

[aleaf]

@spaulins-usgs, thanks, I can supply None now for sfr reaches without connections, but flopy is still writing these reaches as empty lines in packagedata (except for the reach number). For example, if I add a reach 3 at the end of the package data input in the above example with the same values as reach 2, except with None for the cellid, I get this:

BEGIN packagedata
  1  1 1 10      35.50000000       1.76000000       0.44800000     293.20000000       1.00000000       1.00000000  0.037  1       1.00000000  0
  2  1 1 9      35.50000000       1.76000000       0.44800000     293.20000000       1.00000000       1.00000000  0.037  1       1.00000000  0
  3
END packagedata

[spaulins-usgs]

@aleaf-usgs Sorry I forgot to mention that in addition to the problem I fixed with the flopy code there is also a problem with your example. Your example passes in the python value None for a cellid. Flopy is expecting the text 'none', not the python value none ('none' is just text that happens to be the same name as the python value).

When I address your next issue (#835) I will add some additional cellid validation code to flopy so that the python value None returns an 'invalid cellid' error.

[jdhughes-usgs]

@spaulins-usgs can we close this?

[spaulins-usgs]

@jdhughes-usgs yes closing this issue.