Modifying Geometry (top and bottom) in Existing mfusg DISU models?

[javgs-bd]

Let's say I have an existing mfusg model which I can load using flopy.mfusg.MfUsg.load. This will allow me to get the disu package (object) of the model.

I also have a function to read a GSF file and create the arrays needed to build the model unstructured modelgrid object (verts, iverts, xc, yc, ncpl, top, botm).

Now, one would be able to modify the top and botm arrays of the modelgrid object using a plethora of methods (interpolation, constant values, etc, etc), and afterwards I can get all the arrays I would need to reconstruct the DISU object either from the original DISU, or the altered modelgrid object... all but one: the fahl array.

I understand that if I only change the top and botm arrays, then the fahl array should be the only other array that would change. However I cannot get this array from the modelgrid.

There is a function get_connection_area_fahl() but it is within a the UnstructuredModelCell class and I can't understand what it does from the source code (it's based in flopy.mf6.coordinates.modelgrid.ModelCell). There is also a get_unstruct_jagged_array_list() but seems to apply only to mf6 modelgrids.

¿Is there any way I'm missing for getting the updated fahl array without having to code it from scratch for mfusg models?

(btw: I would be glad to share the function that reads the gsf file, it's not that sophisticated but haven't found the time to go through the contributing formal process).

[cnicol-gwlogic]

Hi @javgs-bd - not that I know of sorry. I usually get Algomesh to redo this side of things for me, or I guess if you are using gridgen, a rebuild in that would do it too. I can add this to my list for flopy if you want but I wouldn't expect to get to it any time soon. If I did it, it would have to be very clear that it is only for layer thickness changes; the easiest way might be to track changes in thickness, and adjust fahl accordingly on the fly (rather than try to write something to work it out from cell geometry, which I am guessing is what that mf6 stuff does - better, but more work!).

Chris.

[cnicol-gwlogic]

Actually thinking a little more in this, it's not just Fahl it's also the cl12 array(s) that needs updating.

[javgs-bd]

Many thanks @cnicol-gwlogic. I'll try to code it then. If it works I'll try to add it to a PR.

[javgs-bd]

Hi @cnicol-gwlogic. I don't know if this is of interest, but I'm just reviewing this and am having serious doubts if my initial statement was correct. And am getting totally confused.

In the flopy Gridgen class docs, it states:

cl12 – A 1D vector of the cell connection distances, which are from the center of cell n to its shared face will cell m
fahl – A 1D vector of the cell connection information, which is flow area for a vertical connection and horizontal length for a horizontal connection

If that would be the case, none of these two arrays would change in a cell if the cell top and bottom change.

On the other hand, in the Gridgen report (ofr2014-1109.pdf) Table 9 reads:

Furthermore, the MFUSG (TM6-a45.pdf) report states (page 42) for the first, one-layer model example:

_Record 10 is the FAHL array, which contains the perpendicular flow area for all connections. Values in the FAHL array are either 0, 10,000 or 3,333.33. _

I get the 10,000 value because the bigger cells are 100 x 100. But the smaller ones are 33.33 x 33.33, and those lengths would not yield 3,333.33 in area. If the layer is 1 m thickness the interface area of bigger cells would be 100 and of smaller cells would be 33.33.

Finally, the input/output msfusg manual (mfusg_io_v_1_5) describes again:

FAHL—contains the area of the interface Anm between nodes n and m.

Which is consistent with Figure 7 in TM6-a45.

So I'm I wrong (most likely) or there are two wrong descriptions in the Gridgen flopy docs and in the first example of the MFUSG manual?

[javgs-bd]

Geez, just realized "cell dimensions of 100" also applies for the vertical direction." I guess then that, as it was most likely, the MFUSG manual was right. Still, I guess the description in the Gridgen docs is wrong... and fahl would be the only array to recalculate.

[langevin-usgs]

Sorry for the confusion on this. The FAHL input array for MODFLOW-USG is an area term for both horizontal and vertical connections. The GRIDGEN docstring in flopy is incorrect, as you note. Regarding this post, the FAHL array would need to be recalculated if tops and bottoms were to change.

This is different from MODFLOW 6, where instead of an area term for all connections, the user enters an HWVA (horizontal width / vertical area) array. In this case, the HWVA does not need to change if the tops and bottoms change, but it does have a mix of lengths and areas. Not sure which approach is more confusing!

[javgs-bd]

Many thanks @langevin-usgs!

[cnicol-gwlogic]

Mind bender! Yep, gridgen docstring needs adjustment. I could well be wrong, but I still think CL12 will need adjusting for vertical connections if you change cell tops and / or bottoms - along with FAHL for horizontal connections. For mfusg that is...

Maybe a function roughly along these lines is needed for updating fahl and cl12 for new top/bottom arrays vs old ones. Warning: I have not tested / thought much on it - check my logic if you want to use it.

def update_usg_disu_fahl_cl12(model, new_tops, new_bots):
    # assume disu tops and bots are "old ones" to be updated
    # Note: this assumes ivsd = 1 (so IVC exists) - 
    # otherwise you will need to flag vertical vs horizontal connections another way

    # calculate thickness change factors
    thickness_change_factors = (new_tops - new_bots) / (model.disu.top.array - model.disu.bottom.array)

    # expand old --> new change factors to align with fahl/cl12
    ja_thickness_change_factors = np.repeat(thickness_change_factors, model.disu.iac.array)

    # update disu tops and bots with new ones
    model.disu.top = new_tops   # <--these types of updates should now work for a flat array disu.nodes long
    model.disu.bottom = new_bots

    # update fahl for horizontal connections
    new_fahl = model.disu.fahl.array
    # note: I haven't thought about cln-cln / gwf-cln connections ivc==3/4 
    # (I think this is internal to usg anyway). Also check this logic for staggered / nested connections.
    new_fahl[model.disu.ivc.array != 1] *= ja_thickness_change_factors[model.disu.ivc.array != 1]    
    model.disu.fahl = new_fahl

    # update cl12 for vertical connections
    new_cl12 = model.disu.cl12.array
    new_cl12[model.disu.ivc.array == 1] *= ja_thickness_change_factors[model.disu.ivc.array == 1] 
    model.disu.cl12 = new_cl12

[langevin-usgs]

@cnicol-gwlogic, you are correct about CL12 also needing adjustment. Interestingly, for vertical connections, the CL12 information is redundant with top and bottom information. In both the mfusg and mf6 codes, we may use top and bottom information or CL12 information. As a side note a consistency check should probably be added in the MODFLOW codes to make sure CL12 and top and bottom aren't out of whack...

[cnicol-gwlogic]

Thanks @langevin-usgs. Yeah I also wondered about this redundancy issue / potential for error once I thought through the issue above. I wonder if it is best to override the provided vertical CL12 info with half thickness from the provided cell tops/bottoms (and tell the user that that is what happens) - this is probably safest given most users would be building / checking the rest of the model inputs around that top/bottom info (eg boundary conditions). The only issue is if someone is doing something weird on purpose (I think I have seen John Doherty use an offset vertical CL12 centre before for coal seam gas stuff). But maybe overriding of CL12 could be swtiched off using a keyword option for this unusual case.