Particle Tracking with MODPATH
This notebook shows how to create a particle tracking model using MODPATH.
To-Do
make the examples from a package and from a model layer faster
update toc
add cross section
import os
import sys
import flopy
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
import nlmod
nlmod.util.get_color_logger("INFO")
nlmod.show_versions()
Python version : 3.11.14
NumPy version : 2.4.4
Xarray version : 2026.4.0
Matplotlib version : 3.10.9
Flopy version : 3.10.0
nlmod version : 0.11.3dev
Groundwater Flow Model
We use the groundwater flow model from the 03_local_grid_refinement notebook. Make sure to run that notebook before you run this notebook.
# load lgr model dataset
model_ws = os.path.join("..", "examples", "03_local_grid_refinement")
model_name = "IJm_planeten"
ds = xr.open_dataset(os.path.join(model_ws, f"{model_name}.nc"))
# load simulation and groundwaterflow model
# set exe_name to point to mf6 version in nlmod bin directory
exe_name = os.path.join(os.path.dirname(nlmod.__file__), "bin", "mf6")
if sys.platform.startswith("win"):
exe_name += ".exe"
sim = flopy.mf6.MFSimulation.load("mfsim.nam", sim_ws=model_ws, exe_name=exe_name)
gwf = sim.get_model(model_name=model_name)
loading simulation...
loading simulation name file...
loading tdis package...
loading model gwf6...
loading package disv...
loading package npf...
loading package ic...
loading package oc...
loading package ghb...
loading package drn...
loading package chd...
loading package rch...
loading solution package ijm_planeten...
Modpath
Backward tracking
# list with xy coordinates to start particle tracking from
xy_start = [(101500, 496500), (101500, 499100)]
# create a modpath model
mpf = nlmod.modpath.mpf(gwf)
# create the basic modpath package
_mpfbas = nlmod.modpath.bas(mpf)
# find the nodes for given xy
nodes = nlmod.modpath.xy_to_nodes(xy_start, mpf, ds, layer=5)
# create a particle tracking group at the cell faces
pg = nlmod.modpath.pg_from_fdt(nodes)
# create the modpath simulation file
mpsim = nlmod.modpath.sim(mpf, pg, "backward", gwf=gwf)
adding Package: MPBAS
INFO:nlmod.modpath.modpath.pg_from_fdt:particle group with 9 particle per cell face, 54 particles per cell
adding Package: MPSIM
# run modpath model
nlmod.modpath.write_and_run(mpf, script_path="10_modpath.ipynb")
INFO:nlmod.modpath.modpath.write_and_run:write script 20260513_10_modpath.ipynb to modpath workspace
INFO:nlmod.modpath.modpath.write_and_run:write modpath files to model workspace
Writing packages:
Package: MPBAS
Package: MPSIM
INFO:nlmod.modpath.modpath.write_and_run:run modpath model
FloPy is using the following executable to run the model: ../../../../../../envs/latest/lib/python3.11/site-packages/nlmod/bin/mp7_2_002_provisional
MODPATH Version 7.2.002 PROVISIONAL
Program compiled Oct 16 2023 02:57:43 with IFORT compiler (ver. 20.21.7)
Run particle tracking simulation ...
Processing Time Step 1 Period 6. Time = 6.00000E+00 Steady-state flow
Processing Time Step 1 Period 5. Time = 5.00000E+00 Steady-state flow
Processing Time Step 1 Period 4. Time = 4.00000E+00 Steady-state flow
Processing Time Step 1 Period 3. Time = 3.00000E+00 Steady-state flow
Processing Time Step 1 Period 2. Time = 2.00000E+00 Steady-state flow
Processing Time Step 1 Period 1. Time = 1.00000E+00 Steady-state flow
Particle Summary:
0 particles are pending release.
0 particles remain active.
0 particles terminated at boundary faces.
0 particles terminated at weak sink cells.
0 particles terminated at weak source cells.
108 particles terminated at strong source/sink cells.
0 particles terminated in cells with a specified zone number.
0 particles were stranded in inactive or dry cells.
0 particles were unreleased.
0 particles have an unknown status.
Normal termination.
pdata = nlmod.modpath.load_pathline_data(mpf)
f, ax = plt.subplots(nrows=1, ncols=1, figsize=(10, 10))
ax.set_aspect("equal")
ax = nlmod.plot.modelgrid(ds, ax=ax)
for pid in np.unique(pdata["particleid"]):
pf = pdata[pdata["particleid"] == pid]
ax.plot(pf["x"], pf["y"], color="k", linewidth=0.5)
ax.plot(pf["x"], pf["y"], color="k", linewidth=0.5, label="pathline")
cids = [nlmod.grid.get_icell2d_from_xy(xy[0], xy[1], ds) for xy in xy_start]
ax.plot(
ds.x[cids],
ds.y[cids],
label="start of backwards tracking",
ls="",
marker="o",
color="red",
)
ax.set_title("pathlines")
ax.legend(loc="upper right");
f, ax = plt.subplots(nrows=1, ncols=1, figsize=(12, 4))
for _, pid in enumerate(np.unique(pdata["particleid"])):
pf = pdata[pdata["particleid"] == pid]
x0, y0, z0 = pf[["x", "y", "z"]][0]
distance = np.sqrt((pf["x"] - x0) ** 2 + (pf["y"] - y0) ** 2 + (pf["z"] - z0) ** 2)
ax.plot(pf["time"] / 365.25, distance, label=pid)
ax.set_ylabel("distance [m]")
ax.set_xlabel("time [year]")
ax.set_title("distance travelled per particle")
ax.grid()
Forward tracking
# list with xy coordinates to start particle tracking from
xy_start = [(101500, 496500), (101500, 499100)]
# create a modpath model
mpf = nlmod.modpath.mpf(gwf)
# create the basic modpath package
_mpfbas = nlmod.modpath.bas(mpf)
# find the nodes for given xy
nodes = nlmod.modpath.xy_to_nodes(xy_start, mpf, ds, layer=5)
# create a particle tracking group at the cell faces
# pg = nlmod.modpath.pg_from_pd(nodes, localx=0.5, localy=0.5, localz=1.0)
pg = nlmod.modpath.pg_from_fdt(nodes)
# create the modpath simulation file
mpsim = nlmod.modpath.sim(mpf, pg, "forward")
adding Package: MPBAS
INFO:nlmod.modpath.modpath.pg_from_fdt:particle group with 9 particle per cell face, 54 particles per cell
adding Package: MPSIM
# run modpath model
nlmod.modpath.write_and_run(mpf, script_path="10_modpath.ipynb")
INFO:nlmod.modpath.modpath.remove_output:removed old version of 'mp7_ijm_planeten.mppth'
INFO:nlmod.modpath.modpath.remove_output:removed old version of 'mp7_ijm_planeten.timeseries'
INFO:nlmod.modpath.modpath.remove_output:removed old version of 'mp7_ijm_planeten.mpend'
INFO:nlmod.modpath.modpath.write_and_run:write script 20260513_10_modpath.ipynb to modpath workspace
INFO:nlmod.modpath.modpath.write_and_run:write modpath files to model workspace
Writing packages:
Package: MPBAS
Package: MPSIM
INFO:nlmod.modpath.modpath.write_and_run:run modpath model
FloPy is using the following executable to run the model: ../../../../../../envs/latest/lib/python3.11/site-packages/nlmod/bin/mp7_2_002_provisional
MODPATH Version 7.2.002 PROVISIONAL
Program compiled Oct 16 2023 02:57:43 with IFORT compiler (ver. 20.21.7)
Run particle tracking simulation ...
Processing Time Step 1 Period 1. Time = 1.00000E+00 Steady-state flow
Processing Time Step 1 Period 2. Time = 2.00000E+00 Steady-state flow
Processing Time Step 1 Period 3. Time = 3.00000E+00 Steady-state flow
Processing Time Step 1 Period 4. Time = 4.00000E+00 Steady-state flow
Processing Time Step 1 Period 5. Time = 5.00000E+00 Steady-state flow
Processing Time Step 1 Period 6. Time = 6.00000E+00 Steady-state flow
Particle Summary:
0 particles are pending release.
0 particles remain active.
0 particles terminated at boundary faces.
0 particles terminated at weak sink cells.
0 particles terminated at weak source cells.
108 particles terminated at strong source/sink cells.
0 particles terminated in cells with a specified zone number.
0 particles were stranded in inactive or dry cells.
0 particles were unreleased.
0 particles have an unknown status.
Normal termination.
pdata = nlmod.modpath.load_pathline_data(mpf)
f, axl = plt.subplots(nrows=1, ncols=3, figsize=(30, 10))
for i, ax in enumerate(axl):
ax.set_aspect("equal")
ax = nlmod.plot.modelgrid(ds, ax=ax)
for pid in np.unique(pdata["particleid"]):
pf = pdata[pdata["particleid"] == pid]
ax.plot(pf["x"], pf["y"], color="k", linewidth=0.5)
ax.plot(pf["x"], pf["y"], color="k", linewidth=0.5, label="pathline")
cids = [nlmod.grid.get_icell2d_from_xy(xy[0], xy[1], ds) for xy in xy_start]
ax.plot(
ds.x[cids],
ds.y[cids],
label="start of forward tracking",
ls="",
marker="o",
color="red",
)
ax.set_title("pathlines")
ax.legend(loc="upper right")
if i == 1:
ax.set_xlim(101200, 101700)
ax.set_ylim(498700, 499300)
elif i == 2:
ax.set_xlim(101200, 101700)
ax.set_ylim(496300, 496700)
f, ax = plt.subplots(nrows=1, ncols=1, figsize=(12, 4))
for _, pid in enumerate(np.unique(pdata["particleid"])):
pf = pdata[pdata["particleid"] == pid]
x0, y0, z0 = pf[["x", "y", "z"]][0]
distance = np.sqrt((pf["x"] - x0) ** 2 + (pf["y"] - y0) ** 2 + (pf["z"] - z0) ** 2)
ax.plot(pf["time"] / 365.25, distance, label=pid)
ax.set_ylabel("distance [m]")
ax.set_xlabel("time [year]")
ax.set_title("distance travelled per particle")
ax.grid()
Backward tracking from general head boundaries
# create a modpath model
mpf = nlmod.modpath.mpf(gwf)
# create the basic modpath package
_mpfbas = nlmod.modpath.bas(mpf)
# get the nodes from a package
nodes = nlmod.modpath.package_to_nodes(gwf, "GHB", ibound=mpf.ib)
# create a particle tracking group from cell centers
pg = nlmod.modpath.pg_from_pd(nodes, localx=0.5, localy=0.5, localz=0.5)
# create the modpath simulation file
mpsim = nlmod.modpath.sim(mpf, pg, "backward", gwf=gwf)
adding Package: MPBAS
adding Package: MPSIM
# run modpath model
nlmod.modpath.write_and_run(mpf, script_path="10_modpath.ipynb")
INFO:nlmod.modpath.modpath.remove_output:removed old version of 'mp7_ijm_planeten.mppth'
INFO:nlmod.modpath.modpath.remove_output:removed old version of 'mp7_ijm_planeten.timeseries'
INFO:nlmod.modpath.modpath.remove_output:removed old version of 'mp7_ijm_planeten.mpend'
INFO:nlmod.modpath.modpath.write_and_run:write script 20260513_10_modpath.ipynb to modpath workspace
INFO:nlmod.modpath.modpath.write_and_run:write modpath files to model workspace
Writing packages:
Package: MPBAS
Package: MPSIM
INFO:nlmod.modpath.modpath.write_and_run:run modpath model
FloPy is using the following executable to run the model: ../../../../../../envs/latest/lib/python3.11/site-packages/nlmod/bin/mp7_2_002_provisional
MODPATH Version 7.2.002 PROVISIONAL
Program compiled Oct 16 2023 02:57:43 with IFORT compiler (ver. 20.21.7)
Run particle tracking simulation ...
Processing Time Step 1 Period 6. Time = 6.00000E+00 Steady-state flow
Processing Time Step 1 Period 5. Time = 5.00000E+00 Steady-state flow
Processing Time Step 1 Period 4. Time = 4.00000E+00 Steady-state flow
Processing Time Step 1 Period 3. Time = 3.00000E+00 Steady-state flow
Processing Time Step 1 Period 2. Time = 2.00000E+00 Steady-state flow
Processing Time Step 1 Period 1. Time = 1.00000E+00 Steady-state flow
Particle Summary:
0 particles are pending release.
0 particles remain active.
0 particles terminated at boundary faces.
0 particles terminated at weak sink cells.
0 particles terminated at weak source cells.
2922 particles terminated at strong source/sink cells.
0 particles terminated in cells with a specified zone number.
0 particles were stranded in inactive or dry cells.
0 particles were unreleased.
0 particles have an unknown status.
Normal termination.
pdata = nlmod.modpath.load_pathline_data(mpf)
f, axl = plt.subplots(nrows=1, ncols=3, figsize=(30, 10))
for i, ax in enumerate(axl):
ax.set_aspect("equal")
ax = nlmod.plot.modelgrid(ds, ax=ax)
for pid in np.unique(pdata["particleid"]):
pf = pdata[pdata["particleid"] == pid]
ax.plot(pf["x"], pf["y"], color="k", linewidth=0.5)
ax.plot(pf["x"], pf["y"], color="k", linewidth=0.5, label="pathline")
if i > 0:
cids = np.where((ds["rws_oppwater_cond"] != 0).values)[0]
ax.plot(
ds.x[cids],
ds.y[cids],
label="start of backwards tracking",
ls="",
marker="o",
color="red",
)
ax.set_title("pathlines")
ax.legend(loc="upper right")
if i == 1:
ax.set_xlim(101000, 102000)
ax.set_ylim(498300, 499300)
elif i == 2:
ax.set_xlim(101000, 102000)
ax.set_ylim(496300, 497300)
f, ax = plt.subplots(nrows=1, ncols=1, figsize=(12, 4))
for i, pid in enumerate(np.unique(pdata["particleid"])):
pf = pdata[pdata["particleid"] == pid]
x0, y0, z0 = pf[["x", "y", "z"]][0]
distance = np.sqrt((pf["x"] - x0) ** 2 + (pf["y"] - y0) ** 2 + (pf["z"] - z0) ** 2)
ax.plot(pf["time"] / 365.25, distance, label=pid)
ax.set_xlim(0, 5000)
ax.set_ylabel("distance [m]")
ax.set_xlabel("time [year]")
ax.set_title("distance travelled per particle")
ax.grid()
Forward tracking from each cell in the top layer
Stop after 10 years.
# create a modpath model
mpf = nlmod.modpath.mpf(gwf)
# create the basic modpath package
_mpfbas = nlmod.modpath.bas(mpf)
# get nodes of all cells in the top modellayer
nodes = nlmod.modpath.layer_to_nodes(mpf, 0)
# create a particle tracking group from cell centers
pg = nlmod.modpath.pg_from_pd(nodes, localx=0.5, localy=0.5, localz=0.5)
# create the modpath simulation file
mpsim = nlmod.modpath.sim(mpf, pg, "forward", gwf=gwf, stoptime=10 * 365)
adding Package: MPBAS
adding Package: MPSIM
# run modpath model
nlmod.modpath.write_and_run(mpf, script_path="10_modpath.ipynb")
INFO:nlmod.modpath.modpath.remove_output:removed old version of 'mp7_ijm_planeten.mppth'
INFO:nlmod.modpath.modpath.remove_output:removed old version of 'mp7_ijm_planeten.timeseries'
INFO:nlmod.modpath.modpath.remove_output:removed old version of 'mp7_ijm_planeten.mpend'
INFO:nlmod.modpath.modpath.write_and_run:write script 20260513_10_modpath.ipynb to modpath workspace
INFO:nlmod.modpath.modpath.write_and_run:write modpath files to model workspace
Writing packages:
Package: MPBAS
Package: MPSIM
INFO:nlmod.modpath.modpath.write_and_run:run modpath model
FloPy is using the following executable to run the model: ../../../../../../envs/latest/lib/python3.11/site-packages/nlmod/bin/mp7_2_002_provisional
MODPATH Version 7.2.002 PROVISIONAL
Program compiled Oct 16 2023 02:57:43 with IFORT compiler (ver. 20.21.7)
Run particle tracking simulation ...
Processing Time Step 1 Period 1. Time = 1.00000E+00 Steady-state flow
Processing Time Step 1 Period 2. Time = 2.00000E+00 Steady-state flow
Processing Time Step 1 Period 3. Time = 3.00000E+00 Steady-state flow
Processing Time Step 1 Period 4. Time = 4.00000E+00 Steady-state flow
Processing Time Step 1 Period 5. Time = 5.00000E+00 Steady-state flow
Processing Time Step 1 Period 6. Time = 6.00000E+00 Steady-state flow
Particle Summary:
0 particles are pending release.
1483 particles remain active.
0 particles terminated at boundary faces.
0 particles terminated at weak sink cells.
0 particles terminated at weak source cells.
565 particles terminated at strong source/sink cells.
0 particles terminated in cells with a specified zone number.
0 particles were stranded in inactive or dry cells.
0 particles were unreleased.
0 particles have an unknown status.
Normal termination.
pdata = nlmod.modpath.load_pathline_data(mpf)
f, axl = plt.subplots(nrows=1, ncols=3, figsize=(30, 10))
for i, ax in enumerate(axl):
ax.set_aspect("equal")
ax = nlmod.plot.modelgrid(ds, ax=ax)
for pid in np.unique(pdata["particleid"]):
pf = pdata[pdata["particleid"] == pid]
ax.plot(pf["x"], pf["y"], color="k", linewidth=0.5)
ax.plot(pf["x"], pf["y"], color="k", linewidth=0.5, label="pathline")
if i > 0:
ax.plot(
ds.x.values,
ds.y.values,
label="start of forward tracking",
ls="",
marker="o",
color="red",
)
ax.set_title("pathlines")
ax.legend(loc="upper right")
if i == 1:
ax.set_xlim(101000, 102000)
ax.set_ylim(498300, 499300)
elif i == 2:
ax.set_xlim(101000, 102000)
ax.set_ylim(496300, 497300)
f, ax = plt.subplots(nrows=1, ncols=1, figsize=(12, 4))
for i, pid in enumerate(np.unique(pdata["particleid"])):
pf = pdata[pdata["particleid"] == pid]
x0, y0, z0 = pf[["x", "y", "z"]][0]
distance = np.sqrt((pf["x"] - x0) ** 2 + (pf["y"] - y0) ** 2 + (pf["z"] - z0) ** 2)
ax.plot(pf["time"] / 365.25, distance, label=pid)
ax.set_xlim(0, 11)
ax.set_ylabel("distance [m]")
ax.set_xlabel("time [year]")
ax.set_title("distance travelled per particle")
ax.grid()
