Struggling to use unsaturated flow package (UFZ) Please help!

[AH-CSM]

Hello!
I am trying to add a UFZ package to my modflow 6 model since I am trying to have layer depths above my water table. I understand that there are example problems available in the master folders but they are all for a single layer. I am having trouble developing the modflow 6 UFZ package language lines in a way that I can :

a) feed in a set of vertical k values for multiple layers
b) create and feed in arrays of one values to use for the variables that I will not need since I don't want ET simulated just yet.
c) create an idomain file that includes all layers since it seems that I need this for the UFZ?

I basically want to use the UFZ package in as basic a form as possible just to allow my layer depths to be above my water table.
Any advice would be greatly appreciated. I can share the code I have so far if needed.
Thank you in advance!

[emorway-usgs]

Hello @AH-CSM, I wonder if one of the MODFLOW6 autotests would help provide the example FloPy script you're seeking. Perhaps this example/autotest, which uses FloPy to generate a MODFLOW 6 model with UZF will demonstrate some of what you are looking for?

[AH-CSM]

It would help if my k's were all the same as those examples look to be. However, I am trying to change the k for each layer and build an array with a for loop since my grid size is 150x150 with 5 layers. The masters example helped as I am adding and removing things from that example code to try and make this work for multi layers. An example of what I tried for the pakdata is below:

uzf_pakdata = []
iuzf = 0
for n in range(Nlay):

    if n == 0:
        uzf_kv = 50
        for i in range(nrow):
            for j in range(ncol):
                cellid = [n, i, j]
                uzfcell = [
                    iuzf,
                    0,
                    cellid,
                    1,
                    -1,
                    surf_dep,
                    uzf_kv,
                    thtr * extwc_mult,
                    thts,
                    thti * extwc_mult,
                    eps,
                    "surfrate",
                    ]
                uzf_pakdata.append(uzfcell)
    if n == 1:
        uzf_kv = 1e-1
        for i in range(nrow):
            for j in range(ncol):
                cellid = [n, i, j]
                uzfcell = [
                    iuzf,
                    0,
                    cellid,
                    1,
                    -1,
                    surf_dep,
                    uzf_kv,
                    thtr * extwc_mult,
                    thts,
                    thti * extwc_mult,
                    eps,
                    "surfrate",
                    ]
                uzf_pakdata.append(uzfcell)
##repeats through nlay=4##
iuzf += 1

I may be approaching the package data wrong, but if I can figure this out for the package data, I can do the same loop for the stress period data as well I assume.
Just an FYI I am a SUPER newbie to Python AND Flopy and am struggling to learn all of them at once and I am very behind on how some of these things work and Googling only helps so much. Thank you for understanding and helping!

[emorway-usgs]

I think you're pretty close...what about adding something like this just after the first for loop:

for n in range(Nlay):
    if n==0:
        uzf_kv = 2.2   # or whatever value you want for layer 1
        landflg = 1
    if n==1:
        uzf_kv = 3.2   
        landflg = 0   # will remain 0 for the rest of the layers
    if n==2:
        uzf_kv = 3.5   # or whatever value you want
    if n==3:
        uzf_kv = 4.0   # or whatever value you want
    if n==4:
        uzf_kv = 5.5.0   # or whatever value you want
    
    # Note: I didn't look over the code below too closely, but see at least one item you might need to address 
    # (see commented code below)
    for i in range(nrow):
        for j in range(ncol):
            cellid = (n, i, j)  # switched this to a tuple
            if n < 4:  # I'm assuming all cells in your model domain are active, you will need to adjust this if that's not the case
                vertconn = iuzf + (nrow*ncol)  # connects the iuzf ID that will be immediately below
            elif n== 4:
                vertconn = -1

            uzfcell = [
                iuzf,
                cellid, # made a suggestion here
                landflg,
                vertconn,
                surf_dep,
                uzf_kv,
                thtr * extwc_mult,
                thts,
                thti * extwc_mult,
                eps,
                "surfrate",
                ]
            uzf_pakdata.append(uzfcell)

In your example code, you had uzf_kv = 50 for layer 1 and reduced uzf_kv in layer 2 to 0.1. This may contribute to perching which could make the solver work pretty hard. I mention that in case you notice lengthy runtimes.

[AH-CSM]

Thank you so much for your suggestions. The model will run now and I just have to deal with convergence issues. Would you recommend utilizing an RCH and the UFZ in the same simulation?

[emorway-usgs]

Sorry, @AH-CSM, I'm not willing to make a recommendation on how you set up your model to solve the problem you're working on. However, that's the fun part, so I'll leave it to you.