zonebudget for usg model

[hjia1005]

I am using bf.CellBudgetFile to read the cbb file from a usg model and do a zonebudget in python. The current zonbudusg.exe I have produces a weird result and I want to cross check.

Most of the functions in CellBudgetFile work as expected for the usg model except get_ts(). There is an index error but I don't have to use this function.

Another issue is, so far I can only get net flow from the binary file, which is fine for RCH, EVT, WEL but not enough for storage, ghb, riv and some other flow terms. For these I need both in and out. Any help will be appreciated.

[hjia1005]

Now figure out that the negative terms are OUT and positive terms are IN. Then did my own zone budget in flopy for the usg model and got a consistent result with the budget summary in the model list file. I also used the usgbud2smp.exe from PEST gw data utilities to process the cbb file and the results are the same as what I got in Python.

Here is the problem when I use zonbudusg.exe: we literally apply 100% of rainfall to a mine void but the result shows it receives 0 recharge throughout the model run. Also if I add the recharge in all zones together (which covers the entire model), the number doesn't match what is reported in the model list file.

When compared my result with the result from zonbudusg.exe (downloaded from USGS), all other terms are the same (storage, ghb, riv, drn, wel) except rch and evt. I realized that rch and evt have a different data structure in cbc file from other boundary package flow terms.

What I did was, the first array in rch is the nodes that receive recharge (basically the top cells for each column in the usg model, and their global node number are input to the model as IRCH) and the second array in rch are the flow into that node. But why the two arrays are not zipped as other packages such as ghb, riv and wel? Is there something wrong with this processing? If I am correct, then there could be something wrong with zonbudusg.exe...

Looking forward to someone from the development team clarifying the RCH and EVT packages in usg model. Also, if we use NRCHOP=3, then IRCH is omitted, then in what order we should put the RECH array into the model? @cnicol-gwlogic would you be able to look at my question?

[cnicol-gwlogic]

I don't really follow your question sorry. Mfusg Rch and evt are array-based data sets, whereas ghb etc are list based. That's why they are not stored in flopy in the same way.

If using an unstructured grid with different number of nodes per layer, you really must use nevtop and nrchop option values of 2 - i.e. you have to specify the node numbers you apply rch and evt to. Nrchop/nevtop values of 3 will not work as you expect.

I don't understand your screen shots - what is rch? In flopy you can access the irch and rech arrays using flopy.mfusg.MfusgRch.irch and .rech

Have you looked at your rch file and confirmed that it is in the correct format?

Can't help more than that sorry, as I don't really understand your issue.

[hjia1005]

Hi Chris,

Thank you very much for your reply. Sorry I was in a hurry and didn't include the full code. Variables rch and evt in the screenshot are the data extracted from cbb file using the following code:

CBB = bf.CellBudgetFile('model.cbb')
lst_kstpkper=CBB.get_kstpkper()
for kstpkper in lst_kstpkper:
rch=CBB.get_data(kstpkper=kstpkper,text='RECHARGE')
ghb=CBB.get_data(kstpkper=kstpkper,text='HEAD DEP')
! do the same for other flow terms
! then do something to get the zone budget done

Since I extracted the cell by cell flow, I did the zone budget and exported them to csv files. But the results for rch_in and evt_out are very different from what I got from zonbudusg.exe while other flow terms (storage, ghb, riv and wel), are the same. Of course exactly same zone file was used.

This model has 13 layers, 66120 rnodes and 401091 gnodes. RCH is applied on the top cell for each vertical column, yes we set both NRCHOP and NEVTOP as 2 and specified the node numbers (used gnode in model file and in the sequence of rnode from 1 to 66120), as shown in the first array of rch in the screenshot.

My result is consistent with the budget summary in the model list file. But zonbudusg gives a different result which worries me, especially it reports 0 recharge in one zone where we applied 100% rainfall as recharge. I really want to know which is wrong.

I hope it's clear now. Thanks again for your time.

[cnicol-gwlogic]

Re: "if we use NRCHOP=3, then IRCH is omitted, then in what order we should put the RECH array into the model" - if you look at the mfusg recharge package code, you'll see that for unstructured grids, NRCHOP = 3 means you apply recharge only to the top layer (i.e. layer 1). So the rech array is layer 1 nodes long. If layer 1 does not cover the entire model, I wouldn't use this option unless that is what you want to do (no recharge to layers other than layer 1). From memory, NRCHOP = 1 is the same. I always use 2.

Re: getting a recarray dataset from the recharge / evt cbb arrays for doing your zonebudget checks, you could try something along these lines (note irch is one-based, not zero-based; flopy cbb tools don't seem to convert these, but you can do that if you want):

import flopy.utils.binaryfile as bf
import numpy as np


cbb = bf.CellBudgetFile('dubv1tr01.cbc')
#print(cbb.list_records())
rec = cbb.get_data(text='RECHARGE')
#print(rec)
print(cbb.get_times())
print(cbb.get_kstpkper())
print('ntimes in cbc: ', len(rec))
print(rec[0][0].shape, rec[0][0].min(), rec[0][0].max(), rec[0][0].dtype)
print(rec[0][1].shape, rec[0][1].min(), rec[0][1].max(), rec[0][1].dtype)

# example recarray of kper1 recharge
kper1_data = np.rec.fromarrays(
    [rec[0][0].flatten(), rec[0][1].flatten()],
    dtype=[('irch','i'),('rech','f8')]
    )
print(kper1_data)
print(kper1_data['rech'])
print(kper1_data['irch'])

# example recarray of all kpers' recharge
totims = [
    np.array(
        [totim] * rec[kperkstp_index][1].shape[1])
    for kperkstp_index,totim in enumerate(cbb.get_times())
    ]
irch = [
    rec[kperkstp_index][0].flatten()
    for kperkstp_index,totim in enumerate(cbb.get_times())
    ]
rech = [
    rec[kperkstp_index][1].flatten()
    for kperkstp_index,totim in enumerate(cbb.get_times())
    ]

kper_data = np.rec.fromarrays(
    [np.concatenate(totims),
     np.concatenate(irch),
     np.concatenate(rech)],
    dtype=[('totim',float), ('irch',int),('rech',float)])

print(kper_data)
print(kper_data['totim'], kper_data['totim'].shape)
print(kper_data['irch'], kper_data['irch'].shape)
print(kper_data['rech'], kper_data['rech'].shape)

Re: your zonbudusg.exe differences - do you have any zones == 0? I seem to recall these being ignored (I could be wrong on that though). I doubt there is something randomly wrong like this with zonbudusg.exe. Maybe try a run with all nodes set to zone 1 or somethng and see if you get any differences.

[cnicol-gwlogic]

This is what I was thinking originally too. But I'm not so sure now, as those lines are inside a loop (passing over each array for each save time). The first entry in all usg models should I believe be "constant head" (despite no constant heads being in the model); this entry is n nodes (in the whole model) long, even for a compact budget form. And then if FLOW JA FACE is saved, BUFF/BUFFD increases in size again (38806 is the size of the recharge array, and 611856 is the size of buff/buffd == njag - which would always be way bigger than needed for recharge/evt):

I.e. what I thought was your issue is not an issue as far as I can see. But if you can put together a simple example proving there is an issue I am sure someone @ usgs would look at it.

[hjia1005]

@cnicol-gwlogic I read the code again, now understand and agree with what you just replied. NREAL is the size of full 3d of the entire model.

Now I think the issue could be with the repeated use of array BUFFD. Could it be possible that BUFFD(I)=DZERO reset some of the previously assigned rch/evt to 0? Another temporary array variable should have been used to read the ITYPE==3 arrays and then assign it to BUFFD. What do you think? I have to confirm it is truly an issue and not my mistake before I submit it to USGS :-)

DO 270 I=1,NREAL IF(I.NE.IBUFF(I)) THEN BUFFD(IBUFF(I))=BUFFD(I) BUFFD(I)=DZERO END I

[hjia1005]

@cnicol-gwlogic I can confirm the above, BUFFD(I)=DZERO is the issue. I manually checked my zones. If most or all of the gnode numbers in a specific zone are smaller than the size of evt arrays (66120 in my case) then the accumulated recharge or evaporation is very small or 0. Zone2 to zone4 in my model are like this and zone 4 got 0. If most or all of the gnode numbers in a specific zone are larger than the size of evt arrays (66120 in my case) they are well preserved because there is no chance for them to be re-assigned to 0 by BUFFD(I)=DZERO, such as zone 5 to zone 7 in my model. Now I can submit it to USGS.

@langevin-usgs where should I submit this issue related to zonbugusg.exe?

[cnicol-gwlogic]

I think you are spot on, yep, what you say makes sense. Need another array.

I think if you ordered your IRCH/IEVT arrays in node-ascending order you might not see this issue. But still, there is an issue there I think.

[hjia1005]

It's like a fight between gnode and rnode ~~~

[langevin-usgs]

I've been following this discussion a bit while on vacation. The USGS repository for MODFLOW-USG is located here. It sounds like you may have identified a bug in zonbudusg, which we should fix. If you had a simple flopy script that could set up a MODFLOW-USG simulation and/or if you could provide the MODFLOW-USG and zonbudusg input files as an attachment to a new issue on the MODFLOW-USG repository, we could work on a fix.

[hjia1005]

@langevin-usgs Thanks, the issue is submitted in MODFLOW-USG repository.
https://github.com/MODFLOW-USGS/mfusg/issues/50

[langevin-usgs]

Okay, we will take a look. Not sure of the time frame as it's a busy travel season. Hopefully you have found a workaround in the meantime.