How to read COEFF AI2MD from binary file in the MDT module in MODFLOW-USG

[MartiBayer]

Dear all,
I am working on a simple 1D tansport model with MODFLOW-USG, and I am trying to read the values of COEFF AI2MD from the MDT module which are written in a binary file.

I can read concentration values using this code:

import flopy.utils.binaryfile as bf
import os
folder_path = r"C:\test"
os.chdir(folder_path)
file1 = 'test_CONC001.con'
obj1 = bf.HeadFile(file1, text='CONC')
print('\n\nlist_records() returns:\n')
obj1.list_records()
times = obj1.get_times()
t_test=36525
j=0
for t in times:
if t<t_test:
j+=1
conc_vector_previous_time = obj1.get_data(totim=times[j-1])
print('\n\nread vector of concentration for simulation time =', t_test, '\n')
print(conc_vector_previous_time)

But when I try to read COEFF AI2MD values, with this code:

import flopy.utils.binaryfile as bf
import os
folder_path = r"C:\test"
os.chdir(folder_path)
file2 = 'test_COEFF_AI2MD.dat'
obj2 = bf.HeadFile(file2, text='COEFF AI2MD')

I get this error:
ValueError: Could not determine the precision of the headfile test_COEFF_AI2MD.dat

Any ideas on how can I solve this? attached both binary files (produced by MODFLOW-USG). I have tried precision='single' or precision='double' but it did not work

Thank you very much!!
Martí

PS: I am new here. Could it be an option to upload the files so you can see what happens? Thanks again, any help will be much appreciated!!! I have discovered this very usefull on-line binary file viewer https://iamkate.com/code/binary-file-viewer/ and I can see that the data is in the files I am trying to read.

[langevin-usgs]

@cnicol-gwlogic, any thoughts on this? I don't know the structure of the AI2MD file offhand, but presumably a custom reader could be created if we don't have one yet. @MartiBayer, have you tried the CellBudgetFile reader?

[cnicol-gwlogic]

Not me sorry - I've never used mdt. Looks like it is written to disk via ULASAV, so I can't see why there should be an issue. Might be one for Sorab?

[MartiBayer]

I have not tried the CellButgetFile reader. Is there a version for python?
I have found this
https://rogiersbart.github.io/RMODFLOW/reference/rmf_read_cbc.html
but I was looking for something to implement in the python code I am using.

What I want to do is to compute the mass of contaminant diffused into the matrix, I can do that reading the .out ASCII file, but I only get total values for the entire model, not values on a cel-by-cell basis. I have put the two binary files I menthin in the first post here:
https://1drv.ms/f/c/24b83e0090bdddd5/ErVOLH0L3eBIpgacS7goMAcBuSOhbyF-wpaQHWAiNxMZmg?e=MbdFeC
in case it helps seeing what the issue may be.

Thank you very much for your time!!

[langevin-usgs]

Hey @MartiBayer, CellBudgetFile is a flopy utility. It's worth a shot. If it doesn't work, drop your mdt file (hopefully a small one) here into a github message. You can drop files into these messages.

    precision = "double"
    budobj = flopy.utils.binaryfile.CellBudgetFile(fpth, precision=precision)

Also try precision="single" if double doesn't work.

[langevin-usgs]

I was able to grab your file and try it with CellBudgetFile and didn't work. I'll take a look.

[sorab]

The subroutine that writes the MDT files is almost the same as that which writes the CON files (for concentration). It is in single precision. Could it be the header (it is Character*16 and blank spaces may be different)

[langevin-usgs]

hi @MartiBayer, would it be possible for you to zip up the model input datasets and post them here? Like you, I can also read the concentration file with flopy but I can't read the mdt file.

[MartiBayer]

base_case_4feb (2).zip

Here the model files, let me know if you can open them, I think I put them all in there.
MDT output is written in unit 88

thank you so much for your help!!

[sorab]

Hi @MartiBayer please post the files in native (non-GMS) format.

[MartiBayer]

base_case_4feb_native.zip

Here the files in native text (I hope this is the right format!)

[sorab]

Ok. Here is the scoop. The file on unit 88 was not opened in the mfn (NAME) file. So it did not know form and access for the file. If you open it as in the zip file here it should work.

mfn.zip

[MartiBayer]

I did not get this working.

If I add the file with unit 88 in the name file, then modflow creates the file but it is empty after the simulation (nothing is written).
If I do not add the file with unit 88 in the name file, then modflow creates a file with name "fort.88" and the file has data on it. If I look at it using the online binary file viewer I can see that there is data with header COEFF AI2MD.

As you have seen, I use GMS to generate the files and then run modflow either from GMS of from python with Flopy.

Thanks again

[MartiBayer]

Next things I will try:

1. Use a the USGS compiled version of modflow with the native modflow files adding the unit 88 file in the mfn file as Sorab suggested (so far, I only used the .exe files provided with GMS).
2. Modify the modflow source code to write the MDT file in ASCII format, either using the GMS source code ot the original USGS source code (I guess the later is always the best option).

If you have any other idea that may work, let me know!

Thank you very much guys for your help and inputs! It was a pleasure to interact with you in the first post I write.

I hope I can keep active in this forum and help others too!

Cheers,
Martí

[langevin-usgs]

@MartiBayer, I am happy to help on the flopy side with reading this file; we just need to know the structure. @sorab, can you verify and maybe upload a binary .ai1 file and I will check our readers. One question I have is whether this .ai1 file will be an unstructured or structured file. The .con file appears to be structured.

[sorab]

ai1.zip Here is the file. It is structured since the problem is structured.

[langevin-usgs]

It turns out flopy can't read this file because there are two different types of data in the file ("COEFF AI1MD" and "COEFF AI2MD"). However, this code below shows how you can read the mdt file with some relatively simple python code. Hopefully you can work with this. I tested it with the mdt file that @sorab provided above.

import warnings
import numpy as np
import flopy


def binaryread(file, vartype, shape=(1,), charlen=16):
    """
    Read character bytes, scalar or array values from a binary file.

    Parameters
    ----------
    file : file object
        is an open file object
    vartype : type
        is the return variable type: bytes, numpy.int32,
        numpy.float32, or numpy.float64. Using str is deprecated since
        bytes is preferred.
    shape : tuple, default (1,)
        is the shape of the returned array (shape(1, ) returns a single
        value) for example, shape = (nlay, nrow, ncol)
    charlen : int, default 16
        is the length character bytes.  Note that arrays of bytes
        cannot be returned, only multi-character bytes.  Shape has no
        affect on bytes.

    Raises
    ------
    EOFError
    """

    if vartype == str:
        # handle a hang-over from python2
        warnings.warn(
            "vartype=str is deprecated; use vartype=bytes instead.",
            DeprecationWarning,
        )
        vartype = bytes
    if vartype == bytes:
        # read character bytes of length charlen
        result = file.read(charlen)
        if len(result) < charlen:
            raise EOFError
    else:
        # find the number of values
        nval = np.prod(shape)
        result = np.fromfile(file, vartype, nval)
        if result.size < nval:
            raise EOFError
        if nval != 1:
            result = np.reshape(result, shape)
    return result


if __name__ == "__main__":

    mdt_file = "base_case_4feb.ai1"
    inttype = np.int32
    floattype = np.float32
    dt = [
        ("kstp", inttype),
        ("kper", inttype),
        ("pertim", floattype),
        ("totim", floattype),
        ("text", "S16"),
        ("ncol", inttype),
        ("nrow", inttype),
        ("ilay", inttype),
    ]

    with open (mdt_file, "rb") as f:
        while True:
            result = {}
            for varname, vartype in dt:
                r = binaryread(f, vartype)[0]
                if "text" in varname:
                    r = r.decode("ascii", "strict").strip()
                result[varname] = r
            nrow = result["nrow"]
            ncol = result["ncol"]
            result["data"] = binaryread(f, np.float32, shape=(nrow, ncol))
            print (result)

This is what the program outputs for each record:

{'kstp': np.int32(33), 'kper': np.int32(2), 'pertim': np.float32(36525.0), 'totim': np.float32(60266.25), 'text': 'COEFF AI1MD', 'ncol': np.int32(400), 'nrow': np.int32(1), 'ilay': np.int32(1), 'data': array([[0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
        0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
        0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
        0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
        0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
        0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
        0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
        0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
        0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
        0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
        0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
        0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
        0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
        0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
        0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
        0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
        0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
        0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
        0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
        0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
        0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
        0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
        0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
        0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
        0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]],
      dtype=float32)}

[MartiBayer]

Thank you so much!! I could get this working form the files generated in GMS by writting the MDT file in the same unit as the concentration file (base_case_4feb001.con). Now I can read both CONC and COEFF AI2MD in just one binary file. GREAT!

Now I realize that computing mass in matrix from AI2MD values seems not to be that easy, I thought that the average concentration in the matrix could be computed just as AIOLD2MD/DIFFLENMD (since AIOLD2MD is the integral from 0 to DIFFLENMD, dividing AIOLD2MD by DIFFLENMD should result in the average concentration, but maybe I am missing something somewhere)