Modflow 6 external files not supported?

[tomvansteijn]

I have tried to link to an external file for a Modflow 6 recharge package, without success. It seems the package is trying to read the supplied filename as a list. Is it true that reading from external files is currently not supported?

I could try to write a test if the functionality is there but needs to be fixed.

[langevin-usgs]

@tomvansteijn, it would be good if you included the syntax for what you were trying to do. Without that, it's hard to help.

External files do work with MODFLOW 6 for flopy. Here is an example:

# rch
recharge = {0: {'filename': 'myrecharge.dat', 'data': np.ones((nrow, ncol)) * 0.005}}
rch = flopy.mf6.ModflowGwfrcha(gwf, recharge=recharge)

This will create a new file called myrecharge.dat that has all 0.005 values in it. If you leave out the 'data' key, then the recharge package will just point to that file, which you must have created some other way. The code above will create a recharge package that looks like this:

BEGIN options
  READASARRAYS
END options

BEGIN period  1
  recharge
    OPEN/CLOSE  'myrecharge.dat'  FACTOR  1.0
END period  1

It takes some time getting used to the syntax, but once you do, I think you'll find that nearly all of the MODFLOW 6 capabilities are supported.

[tomvansteijn]

Thanks for the explanation. This looks great. I tried the syntax for Modflow 2005 as in the example notebook for external file handling:

https://github.com/modflowpy/flopy/blob/develop/examples/Notebooks/flopy3_external_file_handling.ipynb

Like this:

rechargefile = r"ref/recharge.ref"  # existing file
rch = flopy.mf6.ModflowGwfrcha(gwf, pname='rch', recharge=rechargefile)

Which obviously could not have worked.

Would it be helpful if I reworked the notebook "flopy3_external_file_handling.ipynb" for Modflow 6 as an addition to the example notebooks? Or is there already an example somewhere of external file input for Modflow 6?

[langevin-usgs]

Much of this is described here: https://github.com/modflowpy/flopy/blob/spaulins-develop-grid/examples/Notebooks/flopy3_mf6_tutorial.ipynb. Though we are clearly going to need a better way to document all of these new MODFLOW 6 capabilities. Right now, we are trying to label all the MODFLOW 6 notebooks to start with flopy3_mf6_xxx so that they can be found more easily among the smorgasbord of notebooks.

Closing this issue. If you want to add another notebook that does a better job describing this, please do so with the flopy3_mf6_xxx naming and submit as a PR. Make sure the notebook follows the syntax at the top of the other notebooks or it will fail on Travis.

[tomvansteijn]

Thanks for the link! Just for the sake of completeness, is there a notebook describing the use of a binary format for external files? I use np.savetxt for the regular external files, but have not found a similar function for a binary format. There is ndarray.tofile:

https://docs.scipy.org/doc/numpy/reference/generated/numpy.ndarray.tofile.html#numpy.ndarray.tofile

But this does not create the format that Modflow expects. Perhaps it is missing a header. In the docs I only find the specification for binary grid, output, budget and LAK, MAW, and SFR packages. However I can add the flag 'binary' to the input dictionary for external files, e.g.:

recharge_ext = [{'filename': rechargedatfile, 'data': recharge, 'binary': True}]
rch = flopy.mf6.ModflowGwfrcha(gwf, pname='rch', recharge=recharge_ext)

And the flag appears in the package file:

BEGIN period  1
  recharge  LAYERED
    OPEN/CLOSE  'data\recharge.dat'  FACTOR  1.0  (BINARY)
END period  1

But the data format (as saved by flopy) is the same as when binary=False.

[langevin-usgs]

There is an example of writing a binary array file here: https://github.com/modflowpy/flopy/blob/develop/examples/Notebooks/flopy3_save_binary_data_file.ipynb. But MODFLOW 6 has a bug in it right now with reading these files. We expect to release version 6.0.4 by the end of the week.

[tomvansteijn]

Thanks for the update! Using version 6.0.4 I was able to use binary files as external files for Modflow packages. Apparently the Modflow 6 executable can only read double precision files. I see this is also specified in the Modflow 6 Input and Output documentation, so this is an expected limitation.

[langevin-usgs]

I looked into this a little more and am hoping that @spaulins-usgs can take a look at this. The following is a simple code that works with MODFLOW 6.0.4, which has a fix for reading binary files.

def write_binary_file(fname, arr):
    text = 'head'
    pertim = np.float64(1.0)
    nrow = arr.shape[0]
    ncol = arr.shape[1]
    precision = 'double'
    header = flopy.utils.BinaryHeader.create(bintype=text, precision=precision,
                                             text=text, nrow=nrow, ncol=ncol,
                                             ilay=1, pertim=pertim,
                                             totim=pertim, kstp=1, kper=1)
    flopy.utils.Util2d.write_bin(arr.shape, fname, arr, header_data=header)
    return
    
sim = flopy.mf6.MFSimulation()
tdis = flopy.mf6.modflow.mftdis.ModflowTdis(sim)
ims = flopy.mf6.modflow.mfims.ModflowIms(sim)
gwf = flopy.mf6.ModflowGwf(sim)

# set up grid information
nlay, nrow, ncol = 3, 12, 10
top = 0.
dz = 1.
botm = [top - (k + 1) * dz for k in range(nlay)]
dis = flopy.mf6.ModflowGwfdis(gwf, nlay=nlay, nrow=nrow, ncol=ncol,
                              top=top, botm=botm)

ic = flopy.mf6.ModflowGwfic(gwf)

k1 = np.random.random((nrow, ncol))

# Layer 2 K from binary file
k2a = np.random.random((nrow, ncol))
fname = 'model.npf.k1.bin'
write_binary_file(os.path.join(model_ws, fname), k2a)
k2 = {'filename': fname, 'binary': True, 'iprn':1}

k3 = np.random.random((nrow, ncol))
npf = flopy.mf6.ModflowGwfnpf(gwf, k=[k1, k2, k3])

chd = flopy.mf6.ModflowGwfchd(gwf, stress_period_data=[[(0, 0, 0), 0.],
                                                       [(0, nrow - 1, ncol - 1), 1.],])
sim.write_simulation()

But this could use a couple tweaks.

1. If you try to access npf.k.array, it doesn't work, presumable because it can't read the binary file. But flopy has code for reading it, so it would be nice to get that to work.

2. It looks like npf.k[1] is returning the values for npf.k[0], which is wrong. This is probably the result of it not being able to load the array in the binary file.

3. It would be nice if we could specify the following

k2 = {'filename': fname, 'binary': True, 'iprn':1, 'data':k2a}

and have flopy automatically do what the write_binary_file function shown here does. This would make it much easier to use binary files.

[spaulins-usgs]

Reading and writing binary array and list files is now supported. The tutorial notebook has an example and t505_test.py has some test cases.