Issues with recharge/drying with multiple layers

[Selsey-sketch]

Hello!
Pretty new to Flopy but have experience with modeling (including MODFLOW). My model is using pinch outs and multiple layers. However, only my top layer is drying. Even if I set recharge = 0, only the top layer dries, which isn't realistic for my conceptual model and seems incorrect for how model layers should behave? Does anyone know more about this or has any suggestions? Thank you, I really appreciate it!! :)

[rmz-oz]

The most likely cause of this ("only layer 1 ever goes dry") is ICELLTYPE in the NPF package. Cells with icelltype=0 are treated as confined: their saturated thickness is held constant, so MODFLOW 6 never makes them dry, no matter how low the head gets. It just reports a head below the cell bottom. In flopy, ModflowGwfnpf uses icelltype=0 by default (mfgwfnpf.py). If only the top layer has icelltype=1 (for example icelltype=[1, 0, 0]), you get exactly what you describe.

From the MF6 input description (gwf-npf.dfn):

ICELLTYPE: flag for each cell that specifies how saturated thickness is treated. 0 means saturated thickness is held constant; >0 means saturated thickness varies with computed head when head is below the cell top

and in the source, the drying loop simply skips those cells (gwf-npf.f90):

! -- cycle if inactive or confined
if (this%ibound(n) == 0) cycle
if (this%icelltype(n) == 0) cycle

I checked it with a small cross section (3 layers, layer 2 pinched out with idomain=-1, bottoms 80/60/0 m, a constant head of 40 m in layer 3, no recharge). The water table should end up at about 40 m, so layers 1 and 2 should both be dry. Results with mf6 6.7.0 and flopy 3.11.0:

NPF setting	layer 1	layer 2
icelltype=[1, 0, 0]	all dry	none dry, head 40 m (below its 60 m bottom)
icelltype=0 (flopy default)	none dry, head 40 m (below its 80 m bottom)	none dry
icelltype=1 (all layers)	all dry	all dry
icelltype=1 + Newton	none flagged dry, heads below bottom	none flagged dry, heads below bottom

So the fix is to make the lower layers convertible (in STO, iconvert=1 as well if it's transient, see gwf-sto.dfn):

npf = flopy.mf6.ModflowGwfnpf(gwf, icelltype=1, k=k, k33=k33)  # or a per-layer list
sto = flopy.mf6.ModflowGwfsto(gwf, iconvert=1, ...)

A few related things that showed up in the same tests:

• With the Newton formulation (in my run newtonoptions="UNDER_RELAXATION", which flopy writes as NEWTON UNDER_RELAXATION in the name file), cells never go dry by design ("Cells will not dry when this option is used", gwf-nam.dfn). You won't see 1e30 values; look for head < cell bottom instead.
• Pinch-outs: with idomain=-1 the cell is a vertical pass-through and doesn't exist in the solution, so it can't be dry or wet (gwf-dis.dfn). In the head file those cells hold 1e30 while dry cells hold -1e30, so a check like abs(h) > 1e29 would count them as dry by mistake. Thin active cells behave like any other cell and follow their icelltype.
• Recharge: once the top cell is dry, RCHA moved the recharge to the first active cell below it in my runs. With fixed_cell=True it stayed on the dry layer 1 cell and the recharge rate was zero.
• If you're on MODFLOW-2005/NWT, the same idea is laytyp in LPF/UPW, also 0 by default in flopy. With LPF laytyp=[1, 0, 0] I got the same result: only layer 1 dry.

A quick check you can run on your model after it runs:

import numpy as np
ict = gwf.npf.icelltype.array
print("icelltype per layer:", [np.unique(ict[k]).tolist() for k in range(ict.shape[0])])
h = gwf.output.head().get_data()
idm = gwf.dis.idomain.array
for k in range(h.shape[0]):
    active = idm[k] > 0
    dry = (h[k] <= -1e29) & active
    below = (h[k] < gwf.dis.botm.array[k]) & ~dry & active
    print(f"layer {k+1}: dry={dry.sum()}  wet but head below bottom={below.sum()}")

If your lower layers already have icelltype=1 and you're not using Newton, can you share your NPF, STO and IMS settings and how you set up the pinch-outs (idomain -1 or thin cells)? Then it's a different cause and worth a closer look.