Unable to load RCH package with inflow concentrations MF-USG model

[javgs-bd]

I'm trying to load a MF-USG model that has inflow concentrations on RCH Package. I can load the model by using forgive = True, check = False with flopy.modflow.Modflow.load. However the recharge package is not loaded.
When setting forgive = False, check = True the following error log indicates that the problem lies in the third line of the RCH package, where variable INCONC is specified.

Below the error log, I am showing some bits of the RCH package format.
I wonder if this has an easy fix, since with this setup the RCH package contains the array of rechrates for each SP followed by an array of rech concentrations.

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
C:\Users\JAVIER~1\AppData\Local\Temp/ipykernel_38636/732944418.py in <module>
     40 #--- https://stackoverflow.com/questions/51734703/extract-heads-from-modflow-usg-binary-output-using-flopy
     41 
---> 42 gwf = flopy.modflow.Modflow.load(
     43                                    model_name,
     44                                    version       = 'mfusg',

~\.conda\envs\flopy-develop\lib\site-packages\flopy\modflow\mf.py in load(cls, f, version, exe_name, verbose, model_ws, load_only, forgive, check)
    884                     else:
    885                         if "check" in package_load_args:
--> 886                             item.package.load(
    887                                 item.filehandle,
    888                                 ml,

~\.conda\envs\flopy-develop\lib\site-packages\flopy\modflow\mfrch.py in load(cls, f, model, nper, ext_unit_dict, check)
    493                             f"   loading rech stress period {iper + 1:3d}..."
    494                         )
--> 495                     t = Util2d.load(
    496                         f,
    497                         model,

~\.conda\envs\flopy-develop\lib\site-packages\flopy\utils\util_array.py in load(cls, f_handle, model, shape, dtype, name, ext_unit_dict, array_free_format, array_format)
   2838         #    array_format = model.array_format
   2839 
-> 2840         cr_dict = Util2d.parse_control_record(
   2841             f_handle.readline(),
   2842             current_unit=curr_unit,

~\.conda\envs\flopy-develop\lib\site-packages\flopy\utils\util_array.py in parse_control_record(line, current_unit, dtype, ext_unit_dict, array_format)
   3027                 npl, fmt, width, decimal = None, None, None, None
   3028         else:
-> 3029             locat = int(line[0:10].strip())
   3030             if isfloat:
   3031                 if len(line) >= 20:

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

RCH Package:

# MODFLOW-USGs Recharge Package
 3 50 CONC 
 1
 1 INCONC
INTERNAL  1.000000e+00  (FREE)  -1  RECHARGE
 1.700000e-05 1.700000e-05 1.700000e-05 1.700000e-05 1.700000e-05 1.700000e-05 1.700000e-05 1.700000e-05 1.700000e-05 
 1.700000e-05 1.700000e-05 1.700000e-05 1.700000e-05 1.700000e-05 1.700000e-05 1.700000e-05 1.700000e-05 1.700000e-05 
 1.700000e-05 1.700000e-05 1.700000e-05 1.700000e-05 1.700000e-05 1.700000e-05 1.700000e-05 1.700000e-05 1.700000e-05 

.
.
.
.
1.700000e-05 1.700000e-05 1.700000e-05 1.700000e-05 1.700000e-05 1.700000e-05 1.700000e-05 1.700000e-05 0.000000e+00 
1.700000e-05 1.700000e-05 1.700000e-05 1.700000e-05 1.700000e-05 1.700000e-05 1.700000e-05 1.700000e-05 1.700000e-05 
1.700000e-05 1.700000e-05 1.700000e-05 1.700000e-05
INTERNAL  1.000000e+00  (FREE)  -1  RECHARGE CONC
1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 
1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 
1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 
.
.
.

[cnicol-gwlogic]

Hi @javgs-bd, I plan on bringing some of this capability in to Flopy. There is a bit to do first to enable autotests for usg-transport however. But keep an eye out - I'll try to get this happening over the next couple of weeks.

[javgs-bd]

Many thanks @cnicol-gwlogic, really appreciate it.

[javgs-bd]

Hi @cnicol-gwlogic, just some notes on this topic that though might be good to pass along. I've noted that something similar is happening for my well package, but now the error is reported as a mismatch in data types.

AssertionError: MfList error: recarray dtype: (numpy.record, [('node', '<i4'), ('flux', '<f4'), ('c01', '<f4')]) doesn't match self dtype: [('node', '<i4'), ('flux', '<f4')]

Two things that might be worth noting is that i) previously, when using the old MODFLOW class and version 'mfusg' these packages were loaded correctly. I guess concentrations (or auxiliary variables) were being ignored (would have to check that), and now with the MfUsg class they are not. and ii) I don't know why I dont have troubles to load a CHD package I am using that does have concentrations, also defined as a c01 aux variable.

I'm leaving the whole traceback below.

In the meantime I am loading the model using the load_only attribute of the load method, and passing only the packages that do not conflict.

Regards!

AssertionError                            Traceback (most recent call last)
C:\Users\JAVIER~1\AppData\Local\Temp/ipykernel_20644/144936415.py in <module>
     40 #--- https://stackoverflow.com/questions/51734703/extract-heads-from-modflow-usg-binary-output-using-flopy
     41 
---> 42 gwf = flopy.mfusg.MfUsg.load(
     43                                    model_name,
     44                                    version       = 'mfusg',

~\.conda\envs\flopy-usg\lib\site-packages\flopy\mfusg\mfusg.py in load(cls, f, version, exe_name, verbose, model_ws, load_only, forgive, check)
    260         cls._set_mfpar_pval(model, ext_unit_dict, ext_pkg_d)
    261 
--> 262         files_successfully_loaded, files_not_loaded = cls._load_packages(
    263             model, ext_unit_dict, ext_pkg_d, load_only, forgive
    264         )

~\.conda\envs\flopy-usg\lib\site-packages\flopy\mfusg\mfusg.py in _load_packages(cls, model, ext_unit_dict, ext_pkg_d, load_only, forgive)
    339                     files_successfully_loaded,
    340                     files_not_loaded,
--> 341                 ) = cls._load_ext_unit_dict_paks(
    342                     model,
    343                     ext_unit_dict,

~\.conda\envs\flopy-usg\lib\site-packages\flopy\mfusg\mfusg.py in _load_ext_unit_dict_paks(cls, model, ext_unit_dict, load_only, item, forgive, files_successfully_loaded, files_not_loaded)
    411                     files_not_loaded.append(item.filename)
    412             else:
--> 413                 cls._ext_unit_d_load(model, ext_unit_dict, item)
    414                 files_successfully_loaded.append(item.filename)
    415                 if model.verbose:

~\.conda\envs\flopy-usg\lib\site-packages\flopy\mfusg\mfusg.py in _ext_unit_d_load(model, ext_unit_dict, ext_unit_d_item)
    456         package_load_args = getfullargspec(ext_unit_d_item.package.load)[0]
    457         if "check" in package_load_args:
--> 458             ext_unit_d_item.package.load(
    459                 ext_unit_d_item.filehandle,
    460                 model,

~\.conda\envs\flopy-usg\lib\site-packages\flopy\mfusg\mfusgwel.py in load(cls, f, model, nper, ext_unit_dict, check)
    462                 model.add_pop_key_list(iunitafr)
    463 
--> 464         wel = cls(
    465             model,
    466             ipakcb=ipakcb,

~\.conda\envs\flopy-usg\lib\site-packages\flopy\mfusg\mfusgwel.py in __init__(self, model, ipakcb, stress_period_data, cln_stress_period_data, dtype, extension, options, binary, unitnumber, filenames, add_package)
    186         filenames = self._prepare_filenames(filenames)
    187 
--> 188         super().__init__(
    189             model,
    190             ipakcb=ipakcb,

~\.conda\envs\flopy-usg\lib\site-packages\flopy\modflow\mfwel.py in __init__(self, model, ipakcb, stress_period_data, dtype, extension, options, binary, unitnumber, filenames, add_package)
    237 
    238         # initialize MfList
--> 239         self.stress_period_data = MfList(
    240             self, stress_period_data, binary=binary
    241         )

~\.conda\envs\flopy-usg\lib\site-packages\flopy\utils\util_list.py in __init__(self, package, data, dtype, model, list_free_format, binary)
     87         self.__data = {}
     88         if data is not None:
---> 89             self.__cast_data(data)
     90         self.__df = None
     91         if list_free_format is None:

~\.conda\envs\flopy-usg\lib\site-packages\flopy\utils\util_list.py in __cast_data(self, data)
    338 
    339                 if isinstance(d, np.recarray):
--> 340                     self.__cast_recarray(kper, d)
    341                 elif isinstance(d, np.ndarray):
    342                     self.__cast_ndarray(kper, d)

~\.conda\envs\flopy-usg\lib\site-packages\flopy\utils\util_list.py in __cast_recarray(self, kper, d)
    393 
    394     def __cast_recarray(self, kper, d):
--> 395         assert d.dtype == self.__dtype, (
    396             "MfList error: recarray dtype: {} doesn't match self dtype: "
    397             "{}".format(d.dtype, self.dtype)

AssertionError: MfList error: recarray dtype: (numpy.record, [('node', '<i4'), ('flux', '<f4'), ('c01', '<f4')]) doesn't match self dtype: [('node', '<i4'), ('flux', '<f4')]

[cnicol-gwlogic]

Thanks @javgs-bd - apologies, I have been too busy to get around to this yet. It is on my list, but it will take some time sorry.