bug: get_structured_faceflows outputting wrong flows

[vcantarella]

Describe the bug

Hi, I created a 2D vertical modflow6 model which has only one row, thus flow should occur only on the x and z directions. When I apply the posprocessing get_structured_faceflows I would expect fff (front face flow) to be 0 in all cells, however at the output I get flf (lower face flow) as a 0 array, and my fff is non zero.

To Reproduce

I am providing an example of a 2D vertical slice model below with the structured face flow output for inspection. Thank you.

import flopy
import os
import matplotlib.pyplot as plt
import numpy as np
#from particle_track import particle_track, prepare_arrays

# Model creation:
name = 'mf6_verif_2'
model_directory = "modflow_verification_2"
sim = flopy.mf6.MFSimulation(
    sim_name=name, exe_name="mf6", version="mf6", sim_ws=model_directory
)

# Simulation time:
tdis = flopy.mf6.ModflowTdis(
    sim, pname="tdis", time_units="DAYS", nper=1, perioddata=[(1.0, 1, 1.0)]
)

# Nam file
model_nam_file = "{}.nam".format(name)

# Groundwater flow object:
gwf = flopy.mf6.ModflowGwf(
    sim,
    modelname=name,
    model_nam_file=model_nam_file,
    save_flows=True,
)

# Grid properties:
Lx = 2000 #problem lenght [m]
Ly = 1 #problem width [m]
H = 100  #aquifer height [m]
delx = 1 #block size x direction
dely = 1 #block size y direction
delz = 1 #block size z direction

nlay = 100

ncol = int(Lx/delx) # number of columns

nrow = int(Ly/dely) # number of layers

# Flopy Discretizetion Objects (DIS)
bottom_array = np.ones((nlay,nrow,ncol))
bottom_range = np.arange(99,-1,-1)
bottom_range = bottom_range[:,np.newaxis,np.newaxis]
bottom_array = bottom_array * bottom_range

dis = flopy.mf6.ModflowGwfdis(
    gwf,
    xorigin = 0.,
    yorigin = 0.,
    nlay=nlay,
    nrow=nrow,
    ncol=ncol,
    delr=dely,
    delc=delx,
    top=100.,
    botm=bottom_array,
)

# Flopy initial Conditions

h0 = 200
start = h0 * np.ones((nlay, nrow, ncol))
ic = flopy.mf6.ModflowGwfic(gwf, pname="ic", strt=start)

# Node property flow

k = 1e-4 * np.ones((nlay, nrow, ncol)) # Model conductivity in m/s
k[9:30,:,300:1701] = 1e-8
npf = flopy.mf6.ModflowGwfnpf(
    gwf,
    icelltype=0, #This we define the model as convertible (water table aquifer)
    k=k,
)

#boundary conditions:

## Constant head    
chd_rec = []

h = np.linspace(101, 110, ncol)
i = 0
for col in range(0, ncol):
    #((layer,row,col),head,iface)
    chd_rec.append(((0, 0, col), h[i], 6))
    i+=1

chd = flopy.mf6.ModflowGwfchd(
    gwf,
    auxiliary=[('iface',)],
    stress_period_data=chd_rec,
    print_input=True,
    print_flows=True,
    save_flows=True,
)

#Output control and Solver
headfile = "{}.hds".format(name)
head_filerecord = [headfile]
budgetfile = "{}.cbb".format(name)
budget_filerecord = [budgetfile]
saverecord = [("HEAD", "ALL"), ("BUDGET", "ALL")]
printrecord = [("HEAD", "LAST")]
oc = flopy.mf6.ModflowGwfoc(
    gwf,
    saverecord=saverecord,
    head_filerecord=head_filerecord,
    budget_filerecord=budget_filerecord,
    printrecord=printrecord,
)
ims = flopy.mf6.ModflowIms(
    sim,
    pname="ims",
    complexity="SIMPLE",
    #linear_acceleration="BICGSTAB",
    outer_maximum = 10,
    inner_maximum = 1500,
    inner_dvclose=1e-3,
    rcloserecord = [0.01, 'STRICT']
)

# Solving
sim.write_simulation()
sim.check()
success, buff = sim.run_simulation()
if not success:
    raise Exception("MODFLOW 6 did not terminate normally.")

## Reading output files

head = gwf.output.head()
head_array = head.get_data()

# Printing output:
fig, ax = plt.subplots(1, 1, figsize=(24, 6), constrained_layout=True)
# first subplot
contour_intervals = np.arange(101, 110, 0.5)

#plotting cross section of vertical model:

ax.set_title("Head Results")
modelmap = flopy.plot.PlotCrossSection(model=gwf, ax=ax, line = {"row": 0})
pa = modelmap.plot_array(np.log(k))
quadmesh = modelmap.plot_bc("CHD")
linecollection = modelmap.plot_grid(lw=0.05, color="0.5")
contours = modelmap.contour_array(
    head_array,
    levels=contour_intervals,
    colors="black",
)
ax.clabel(contours, fmt="%2.1f")
cb = plt.colorbar(pa, shrink=0.5, ax=ax)

'''
Structure Face Flows
'''
budget = gwf.output.budget()
flow_ja_face = budget.get_data(text = 'FLOW-JA-FACE')[0]
frf, fff, flf = flopy.mf6.utils.postprocessing.get_structured_faceflows(flow_ja_face, grb_file=os.path.join(model_directory,gwf.name+'.dis.grb'), verbose=True)

print(frf)
print(fff)
print(flf)

Expected behavior flf should be a non zero array and fff should be a zero array.

Screenshots

Desktop (please complete the following information):

• OS: [e.g. iOS] Windows
• Version : flopy: 3.3.7

Thanks !

[vcantarella]

Looking at the source code of the function I think the problem is in this snippet:

    # fill flow terms
    vmult = [-1.0, -1.0, -1.0]
    flows = [frf, fff, flf]
    for n in range(grb.nodes):
        i0, i1 = ia[n] + 1, ia[n + 1]
        for j in range(i0, i1):
            jcol = ja[j]
            if jcol > n:
                if jcol == n + 1:
                    ipos = 0
                elif jcol == n + grb.ncol: # <======HERE
                    ipos = 1
                else:
                    ipos = 2
                flows[ipos][n] = vmult[ipos] * flowja[j]
    return frf.reshape(shape), fff.reshape(shape), flf.reshape(shape)

It seems like it already assumes that jcol == n+grb.ncol means it is in the next row, but in this case it should jump straight to the next layer.

[langevin-usgs]

Thanks for picking this up @w-bonelli. I looked at this a bit and it should be easy to fix, but it is trickier than it seems. Curious to see what you come up with.

[jdhughes-usgs]

@w-bonelli we have run into edge case issues with function before. It seems we should have better testing for this. Probably need cases (1, nrow, 1), (1, 1, ncol), (1, nrow, ncol), (2, nrow, 1), (2, 1, ncol), (2, nrow, ncol), (3, nrow, 1), (3, 1, ncol), and (3, nrow, ncol) to cover all bases.

[wpbonelli]

Changes were needed to the way we determine which face each flow belongs to, while iterating through the intercell flow matrix to build right/front/lower face arrays.

In the given model, flows with m - n == ncol go to lower face as noted by @vcantarella. Previously the logic was

• 0 (right) if m - n == 1
• 1 (front) if m - n == ncol
• 2 (lower) otherwise

With #1968 it becomes

• 0 (right) if m - n == 1
• 2 (lower) if m - n == ncpl
• 1 (front) otherwise

But this is not generic, it only works if the grid is extended (has multiple cells) in all 3 dimensions. Special treatment is added for "degenerate" cases among those listed above by @jdhughes-usgs