cross_section plot_pathline sorts trace by x coordinate

[Hugovdberg]

I'm plotting pathlines in a simple 2D model that initially move in the positive x-direction but then turn around to end up on the other boundary. When I plot it with plain matplotlib as follows I get this, which I would expect:

plt.plot(trace['x'], trace['z'])

However, the resultant plot looks similar to this, although the projection on the grid alters the x coordinates of the raw trace slightly:

sort_trace = np.sort(trace, axis=0)
plt.plot(sort_trace['x'], sort_trace['z'])

A comparison of two slices of the same trace using the actual PlotCrossSection:

Full code

import flopy
import matplotlib.pyplot as plt
import numpy as np

sim_name = "river_2D"
sim_ws = "river_2D"

sim = flopy.mf6.MFSimulation(sim_name=sim_name, sim_ws=sim_ws)

tdis = flopy.mf6.ModflowTdis(
    simulation=sim, time_units="DAYS"  # , start_date_time="2021-01-01 00:00"
)

ims = flopy.mf6.ModflowIms(simulation=sim, linear_acceleration="BICGSTAB")

gwf = flopy.mf6.ModflowGwf(simulation=sim, newtonoptions="NEWTON", save_flows=True)

nlay = 15
nrow = 1
ncol = 20
ifaces = np.linspace(10, -5, num=nlay + 1)
dis = flopy.mf6.ModflowGwfdis(
    model=gwf,
    length_units="METERS",
    nlay=nlay,
    nrow=nrow,
    ncol=ncol,
    delr=2,
    top=ifaces[0],
    botm=ifaces[1:],
    idomain=1,
)

ks = 0.5 * np.ones((nlay, nrow, ncol))
ks[:3, 0, 8:11] = 1000
ks[3, 0, 8:11] = 0.01
ks[9] = 0.0005
npf = flopy.mf6.ModflowGwfnpf(model=gwf, k=ks, k33=ks / 10,)

ic = flopy.mf6.ModflowGwfic(model=gwf, strt=8)

chd = flopy.mf6.ModflowGwfchd(
    model=gwf,
    stress_period_data=[
        ((lay, 0, col), 7 + col / 19) for lay in range(nlay) for col in [0, 19]
    ]
    + [((lay, 0, col), 9.5) for lay in range(3) for col in range(8, 11)],
)


hfb = flopy.mf6.ModflowGwfhfb(
    model=gwf,
    stress_period_data=[
        [(0, 0, 7), (0, 0, 8), 0.05],
        [(1, 0, 7), (1, 0, 8), 0.05],
        [(2, 0, 7), (2, 0, 8), 0.05],
        [(0, 0, 10), (0, 0, 11), 0.05],
        [(1, 0, 10), (1, 0, 11), 0.05],
        [(2, 0, 10), (2, 0, 11), 0.05],
    ],
)


wel = flopy.mf6.ModflowGwfwel(
    model=gwf,
    stress_period_data=[
        ((12, 0, 4), -0.1),
        ((13, 0, 4), -0.2),
        ((14, 0, 4), -0.1),
        ((12, 0, 15), 0.1),
        ((13, 0, 15), 0.2),
        ((14, 0, 15), 0.1),
    ],
)

oc = flopy.mf6.ModflowGwfoc(
    model=gwf,
    head_filerecord="model.hds",
    budget_filerecord="model.cbb",
    saverecord=[["HEAD", "ALL"], ["BUDGET", "ALL"]],
)

sim.write_simulation()
sim.run_simulation()

heads = flopy.utils.HeadFile("river_2D/model.hds")

final_heads = heads.get_data([0, 0])


# create forward particle group
sd = flopy.modpath.CellDataType(
    drape=0, columncelldivisions=1, rowcelldivisions=1, layercelldivisions=1
)
p = flopy.modpath.LRCParticleData(
    subdivisiondata=[sd], lrcregions=[np.array([[0, 0, 0, nlay - 1, 0, ncol - 1]])]
)
pg2 = flopy.modpath.ParticleGroupLRCTemplate(
    particlegroupname="PG3", particledata=p, filename="ex01a.pg3.sloc"
)

particlegroups = [pg2]
mp = flopy.modpath.Modpath7(
    modelname="model_mp",
    flowmodel=gwf,
    headfilename="model.hds",
    budgetfilename="model.cbb",
    model_ws=sim_ws,
)
mpbas = flopy.modpath.Modpath7Bas(mp, porosity=0.3)
mpsim = flopy.modpath.Modpath7Sim(
    mp,
    simulationtype="combined",
    trackingdirection="forward",
    weaksinkoption="pass_through",
    weaksourceoption="pass_through",
    budgetoutputoption="summary",
    stoptimeoption="extend",
    particlegroups=particlegroups,
)

# write modpath datasets
mp.write_input()

# run modpath
mp.run_model()

p = flopy.utils.PathlineFile("river_2D/model_mp.mppth")
trace = p.get_data(partid=50)

fig, (ax, ax1) = plt.subplots(nrows=2, figsize=(19, 8))
pxs = flopy.plot.PlotCrossSection(gwf, line={"row": 0}, ax=ax)
# c = pxs.plot_array(final_heads, cmap="viridis", head=final_heads)
pxs.plot_grid()
pxs.plot_bc("CHD", head=final_heads, color="none", edgecolor="red", hatch="\\")
pxs.plot_pathline(trace, method="cell", travel_time=65, colors="blue")
# pxs.plot_bc("WEL", head=final_heads, color="none", edgecolor="blue", hatch="||")
ax.set_aspect("equal")
ax.set_title("Trace with strictly increasing x coordinates (max travel time 65 days)")
# ax.grid("off")
# plt.colorbar(c)
pxs = flopy.plot.PlotCrossSection(gwf, line={"row": 0}, ax=ax1)
# c = pxs.plot_array(final_heads, cmap="viridis", head=final_heads)
pxs.plot_grid()
pxs.plot_bc("CHD", head=final_heads, color="none", edgecolor="red", hatch="\\")
pxs.plot_pathline(trace, method="cell", travel_time=85, colors="blue")
# pxs.plot_bc("WEL", head=final_heads, color="none", edgecolor="blue", hatch="||")
ax1.set_aspect("equal")
ax1.set_title(
    "Trace without strictly increasing x coordinates (max travel time 85 days)"
)
fig.savefig("cross_section.png", bbox_inches="tight")

[jlarsen-usgs]

@Hugovdberg I'll take a look at this.

[jlarsen-usgs]

@Hugovdberg

I've put together a fix for this. The code as it's written intersects and reprojects the particles on to the cross section. When it does this the pathline locations are stored in the cross section grid cell drawing order. This works well if the particles do not reverse direction over the course of the simulation, but as you've demonstrated if they do, the method does not draw the path correctly.

The new changes sort the stored values by travel time before plotting. With these changes, here is the result

.
Thanks for supplying a detailed issue description! I'll open a pull request and submit these changes today.

[Hugovdberg]

Thanks for the quick fix! And regarding the description, I know how frustrating it is to have an issue filed "it just doesn't work, sometimes".