diff --git a/notebooks/ex-prt-mp7-p03.ipynb b/notebooks/ex-prt-mp7-p03.ipynb
new file mode 100644
index 00000000..637b3b9a
--- /dev/null
+++ b/notebooks/ex-prt-mp7-p03.ipynb
@@ -0,0 +1,1196 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "id": "ed827fd3",
+ "metadata": {},
+ "source": [
+ "## Particle tracking through a transient flow system\n",
+ "\n",
+ "Application of a MODFLOW 6 particle-tracking (PRT) model to solve example 3 from the MODPATH 7 documentation.\n",
+ "\n",
+ "This is a transient MODFLOW 6 simulation with a flow system similar to examples 1 and 2.\n",
+ "Starting at 90,000 days, a 25 particles are released from 4 cells in the grid's top left corner.\n",
+ "This occurs every 20 days for 200 days. 1000 particles are released in total.\n",
+ "The model's 2 wells begin to pump at a constant rate 100,000 days into the simulation,\n",
+ "and continue to pump at a constant rate for the remainder of the simulation.\n",
+ "\n",
+ "The original problem specifies three stress periods:\n",
+ "\n",
+ "| Stress period | Type | Time steps | Length (days) | Multiplier |\n",
+ "|:--------------|:-------------|:-----------|:--------------|:-----------|\n",
+ "| 1 | steady-state | 1 | 100000 | 1 |\n",
+ "| 2 | transient | 10 | 36500 | 1.5 |\n",
+ "| 3 | steady-state | 1 | 100000 | 1 |\n",
+ "\n",
+ "The original solution achieves particle release timing by setting MODPATH 7's reference time option to 2 and value to 0.9. Since PRT does not yet support fractional time steps, we modify the problems' time discretization to 5 stress periods, with a finer time discretization in the 2nd for particle release:\n",
+ "\n",
+ "| Stress period | Type | Time steps | Length (days) | Multiplier |\n",
+ "|:--------------|:-------------|:-----------|:--------------|:-----------|\n",
+ "| 1 | steady-state | 1 | 90000 | 1 |\n",
+ "| 2 | steady-state | 10 | 200 | 1 |\n",
+ "| 3 | steady-state | 1 | 9800 | 1 |\n",
+ "| 4 | transient | 10 | 36500 | 1.5 |\n",
+ "| 5 | steady-state | 1 | 100000 | 1 |\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "id": "0b23bcc9",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "perioddata = [\n",
+ " # perlen, nstp, tsmult\n",
+ " (90000, 1, 1),\n",
+ " (200, 10, 1),\n",
+ " (9800, 1, 1),\n",
+ " (36500, 10, 1.5),\n",
+ " (100000, 1, 1),\n",
+ "]"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "630777b9",
+ "metadata": {},
+ "source": [
+ "\n",
+ "First import dependencies."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "id": "64aee853",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import os\n",
+ "import sys\n",
+ "import matplotlib.pyplot as plt\n",
+ "import flopy\n",
+ "import numpy as np\n",
+ "import pandas as pd\n",
+ "from pathlib import Path\n",
+ "\n",
+ "sys.path.append(os.path.join(\"..\", \"common\"))\n",
+ "import config"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "ba2958a5",
+ "metadata": {},
+ "source": [
+ "Setup model name and workspace variables."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "id": "0969a70c",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "ws = config.base_ws\n",
+ "sim_name = \"mp7-p03\"\n",
+ "example_name = \"ex-prt-\" + sim_name\n",
+ "sim_ws = Path(ws) / example_name\n",
+ "\n",
+ "nm_mf6 = sim_name\n",
+ "nm_prt = sim_name + \"_prt\"\n",
+ "\n",
+ "headfile = f\"{sim_name}.hds\"\n",
+ "budgetfile = f\"{sim_name}.cbb\"\n",
+ "budgetfile_prt = f\"{nm_prt}.cbb\"\n",
+ "trackcsvfile_prt = f\"{nm_prt}.trk.csv\"\n",
+ "pathlinefile_mp7 = f\"{sim_name}_mp.mppth\"\n",
+ "endpointfile_mp7 = f\"{sim_name}_mp.mpend\""
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "19845f3a",
+ "metadata": {},
+ "source": [
+ "Define some shared variables, starting with discretization and other flow model data."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "id": "c55171b8",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "nlay, nrow, ncol = 3, 21, 20\n",
+ "delr = delc = 500.0\n",
+ "top = 350.0\n",
+ "botm = [220.0, 200.0, 0.0]\n",
+ "laytyp = [1, 0, 0]\n",
+ "kh = [50.0, 0.01, 200.0]\n",
+ "kv = [10.0, 0.01, 20.0]\n",
+ "rchv = 0.005\n",
+ "riv_h = 320.0\n",
+ "riv_z = 317.0\n",
+ "riv_c = 1.0e5"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "60c05f6b",
+ "metadata": {},
+ "source": [
+ "Define well data. Although this notebook will refer to layer/row/column indices starting at 1, indices in FloPy (and more generally in Python) are zero-based. A negative discharge indicates pumping, while a positive value indicates injection."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "id": "4cbb2742",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "wells = [\n",
+ " # layer, row, col, discharge\n",
+ " (0, 10, 9, -75000),\n",
+ " (2, 12, 4, -100000),\n",
+ "]"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "8c4019e7",
+ "metadata": {},
+ "source": [
+ "Define the drain location."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "id": "aa7fb4c2",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "drain = (0, 14, (9, 20))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "0e756027",
+ "metadata": {},
+ "source": [
+ "Configure locations for particle tracking to terminate. We have three explicitly defined termination zones:\n",
+ "\n",
+ "- 2: the well in layer 1, at row 11, column 10\n",
+ "- 3: the well in layer 3, at row 13, column 5\n",
+ "- 4: the drain in layer 1, running through row 15 from column 10-20\n",
+ "- 5: the river in layer 1, running through column 20\n",
+ "\n",
+ "MODFLOW 6 reserves zone number 1 to indicate that particles may move freely within the zone."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "id": "6ff855b9",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "zone_maps = []\n",
+ "\n",
+ "# zone 1 is the default (non-terminating)\n",
+ "def fill_zone_1():\n",
+ " return np.ones((nrow, ncol), dtype=np.int32)\n",
+ "\n",
+ "# zone map for layer 1\n",
+ "za = fill_zone_1()\n",
+ "za[wells[0][1:3]] = 2 # well\n",
+ "za[drain[1], drain[2][0] : drain[2][1]] = 4 # drain\n",
+ "za[:, ncol - 1] = 5 # river\n",
+ "zone_maps.append(za)\n",
+ "\n",
+ "# zone map for layer 2 (use constant to apply zone 1 to all cells)\n",
+ "zone_maps.append(1)\n",
+ "\n",
+ "# zone map for layer 3\n",
+ "za = fill_zone_1()\n",
+ "za[wells[1][1:3]] = 3 # well\n",
+ "zone_maps.append(za)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "5883261f",
+ "metadata": {},
+ "source": [
+ "TODO: plot zone map\n",
+ "\n",
+ "Define particles to track. Particles are released from the top of a 2x2 square of cells in the upper left of the midel grid's top layer."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "id": "3a07c5c0",
+ "metadata": {
+ "lines_to_end_of_cell_marker": 0,
+ "lines_to_next_cell": 1
+ },
+ "outputs": [],
+ "source": [
+ "\n",
+ "rel_minl = rel_maxl = 0\n",
+ "rel_minr = 2\n",
+ "rel_maxr = 3\n",
+ "rel_minc = 2\n",
+ "rel_maxc = 3\n",
+ "celldata = flopy.modpath.CellDataType(\n",
+ " drape=0,\n",
+ " rowcelldivisions=5,\n",
+ " columncelldivisions=5,\n",
+ " layercelldivisions=1,\n",
+ ")\n",
+ "lrcregions = [[rel_minl, rel_minr, rel_minc, rel_maxl, rel_maxr, rel_maxc]]\n",
+ "lrcpd = flopy.modpath.LRCParticleData(\n",
+ " subdivisiondata=celldata,\n",
+ " lrcregions=lrcregions,\n",
+ ")\n",
+ "pg = flopy.modpath.ParticleGroupLRCTemplate(\n",
+ " particlegroupname=\"PG1\",\n",
+ " particledata=lrcpd,\n",
+ " filename=f\"{sim_name}.pg1.sloc\",\n",
+ " releasedata=(10, 0, 20)\n",
+ ")\n",
+ "pgs = [pg]\n",
+ "defaultiface = {\"RECHARGE\": 6, \"ET\": 6}"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "6f0b870b",
+ "metadata": {},
+ "source": [
+ "Define a function to create the MODFLOW 6 simulation, which will include a groundwater flow (GWF) model and a particle tracking (PRT) model."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "id": "5c6edfa2",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def build_mf6_model():\n",
+ " print(\"Building MODFLOW 6 model\")\n",
+ "\n",
+ " # simulation\n",
+ " sim = flopy.mf6.MFSimulation(\n",
+ " sim_name=sim_name,\n",
+ " exe_name=\"mf6\",\n",
+ " version=\"mf6\",\n",
+ " sim_ws=sim_ws\n",
+ " )\n",
+ "\n",
+ " # temporal discretization\n",
+ " tdis = flopy.mf6.modflow.mftdis.ModflowTdis(\n",
+ " sim,\n",
+ " pname=\"tdis\",\n",
+ " time_units=\"DAYS\",\n",
+ " nper=len(perioddata),\n",
+ " perioddata=perioddata\n",
+ " )\n",
+ "\n",
+ " # groundwater flow (gwf) model\n",
+ " model_nam_file = f\"{sim_name}.nam\"\n",
+ " gwf = flopy.mf6.ModflowGwf(\n",
+ " sim, modelname=sim_name, model_nam_file=model_nam_file, save_flows=True\n",
+ " )\n",
+ "\n",
+ " # iterative model solver (ims) package\n",
+ " ims = flopy.mf6.modflow.mfims.ModflowIms(\n",
+ " sim,\n",
+ " pname=\"ims\",\n",
+ " complexity=\"SIMPLE\",\n",
+ " )\n",
+ "\n",
+ " # grid discretization\n",
+ " dis = flopy.mf6.modflow.mfgwfdis.ModflowGwfdis(\n",
+ " gwf,\n",
+ " pname=\"dis\",\n",
+ " nlay=nlay,\n",
+ " nrow=nrow,\n",
+ " ncol=ncol,\n",
+ " length_units=\"FEET\",\n",
+ " delr=delr,\n",
+ " delc=delc,\n",
+ " top=top,\n",
+ " botm=botm,\n",
+ " )\n",
+ "\n",
+ " # initial conditions\n",
+ " ic = flopy.mf6.modflow.mfgwfic.ModflowGwfic(gwf, pname=\"ic\", strt=320)\n",
+ "\n",
+ " # node property flow\n",
+ " npf = flopy.mf6.modflow.mfgwfnpf.ModflowGwfnpf(\n",
+ " gwf, pname=\"npf\", icelltype=laytyp, k=kh, k33=kv, save_flows=True, save_specific_discharge=True, save_saturation=True\n",
+ " )\n",
+ "\n",
+ " # recharge\n",
+ " rch = flopy.mf6.modflow.mfgwfrcha.ModflowGwfrcha(gwf, recharge=rchv)\n",
+ "\n",
+ " # storage\n",
+ " sto = flopy.mf6.modflow.mfgwfsto.ModflowGwfsto(\n",
+ " gwf, save_flows=True, iconvert=1, ss=0.0001, sy=0.1,\n",
+ " steady_state={0: True, 4: True},\n",
+ " transient={3: True},\n",
+ " )\n",
+ "\n",
+ " # wells\n",
+ " def no_flow(w):\n",
+ " return w[0], w[1], w[2], 0\n",
+ "\n",
+ " nf_wells = [no_flow(w) for w in wells] \n",
+ " wel = flopy.mf6.modflow.mfgwfwel.ModflowGwfwel(\n",
+ " gwf,\n",
+ " maxbound=2,\n",
+ " stress_period_data={0: nf_wells, 1: nf_wells, 2: nf_wells, 3: wells, 4: wells},\n",
+ " )\n",
+ "\n",
+ " # river\n",
+ " riv_iface = 6\n",
+ " riv_iflowface = -1\n",
+ " rd = [[(0, i, ncol - 1), riv_h, riv_c, riv_z, riv_iface, riv_iflowface] for i in range(nrow)]\n",
+ " flopy.mf6.modflow.mfgwfriv.ModflowGwfriv( \n",
+ " gwf,\n",
+ " auxiliary=[\"iface\", \"iflowface\"],\n",
+ " stress_period_data={0: rd}\n",
+ " )\n",
+ "\n",
+ " # drain (set auxiliary IFACE var to 6 for top of cell)\n",
+ " drn_iface = 6\n",
+ " drn_iflowface = -1\n",
+ " dd = [\n",
+ " [drain[0], drain[1], i + drain[2][0], 322.5, 100000.0, drn_iface, drn_iflowface]\n",
+ " for i in range(drain[2][1] - drain[2][0])\n",
+ " ]\n",
+ " drn = flopy.mf6.modflow.mfgwfdrn.ModflowGwfdrn(\n",
+ " gwf,\n",
+ " auxiliary=[\"iface\", \"iflowface\"],\n",
+ " stress_period_data={0: dd})\n",
+ "\n",
+ " # output control\n",
+ " oc = flopy.mf6.modflow.mfgwfoc.ModflowGwfoc(\n",
+ " gwf,\n",
+ " pname=\"oc\",\n",
+ " saverecord=[(\"HEAD\", \"ALL\"), (\"BUDGET\", \"ALL\")],\n",
+ " head_filerecord=[headfile],\n",
+ " budget_filerecord=[budgetfile],\n",
+ " )\n",
+ "\n",
+ " # Instantiate the MODFLOW 6 prt model\n",
+ " prt = flopy.mf6.ModflowPrt(\n",
+ " sim, modelname=nm_prt, model_nam_file=\"{}.nam\".format(nm_prt)\n",
+ " )\n",
+ "\n",
+ " # Instantiate the MODFLOW 6 prt discretization package\n",
+ " flopy.mf6.modflow.mfgwfdis.ModflowGwfdis(\n",
+ " prt, pname=\"dis\",\n",
+ " nlay=nlay, nrow=nrow, ncol=ncol,\n",
+ " length_units=\"FEET\",\n",
+ " delr=delr, delc=delc,\n",
+ " top=top, botm=botm,\n",
+ " )\n",
+ "\n",
+ " # Instantiate the MODFLOW 6 prt model input package\n",
+ " porosity = 0.1\n",
+ " flopy.mf6.ModflowPrtmip(prt, pname=\"mip\", porosity=porosity)\n",
+ "\n",
+ " releasepts = list(lrcpd.to_prp(gwf.modelgrid))\n",
+ " flopy.mf6.ModflowPrtprp(\n",
+ " prt, pname=\"prp1\", filename=\"{}_1.prp\".format(nm_prt),\n",
+ " nreleasepts=len(releasepts), packagedata=releasepts,\n",
+ " perioddata={\n",
+ " 0: [],\n",
+ " 1: [\"STEPS\"] + [str(i + 1) for i in range(10)],\n",
+ " 2: [],\n",
+ " 3: [],\n",
+ " 4: []\n",
+ " },\n",
+ " )\n",
+ "\n",
+ " # Instantiate the MODFLOW 6 prt output control package\n",
+ " flopy.mf6.ModflowPrtoc(\n",
+ " prt,\n",
+ " pname=\"oc\",\n",
+ " budget_filerecord=[budgetfile_prt],\n",
+ " trackcsv_filerecord=[trackcsvfile_prt],\n",
+ " saverecord=[(\"BUDGET\", \"ALL\")],\n",
+ " )\n",
+ "\n",
+ " # Instantiate the MODFLOW 6 prt flow model interface\n",
+ " flopy.mf6.ModflowPrtfmi(prt, packagedata=[\n",
+ " (\"GWFHEAD\", headfile),\n",
+ " (\"GWFBUDGET\", budgetfile)\n",
+ " ])\n",
+ "\n",
+ " # Create the MODFLOW 6 gwf-prt model exchange\n",
+ " flopy.mf6.ModflowGwfprt(\n",
+ " sim, exgtype=\"GWF6-PRT6\",\n",
+ " exgmnamea=nm_mf6, exgmnameb=nm_prt,\n",
+ " filename=\"{}.gwfprt\".format(nm_mf6),\n",
+ " )\n",
+ "\n",
+ " # Create an explicit model solution (EMS) for the MODFLOW 6 prt model\n",
+ " ems = flopy.mf6.ModflowEms(\n",
+ " sim, pname=\"ems\",\n",
+ " filename=f\"{nm_prt}.ems\",\n",
+ " )\n",
+ " sim.register_solution_package(ems, [prt.name])\n",
+ "\n",
+ " return sim"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "d7334d85",
+ "metadata": {},
+ "source": [
+ "Define a function create the MODPATH 7 model & simulation."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 10,
+ "id": "e576d2a8",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def build_mp7_model(gwf):\n",
+ " print(\"Building MODPATH 7 model...\")\n",
+ "\n",
+ " mp = flopy.modpath.Modpath7(\n",
+ " modelname=f\"{sim_name}_mp\",\n",
+ " flowmodel=gwf,\n",
+ " exe_name=\"mp7\",\n",
+ " model_ws=sim_ws,\n",
+ " )\n",
+ " mpbas = flopy.modpath.Modpath7Bas(mp, porosity=0.1, defaultiface=defaultiface)\n",
+ " mpsim = flopy.modpath.Modpath7Sim(\n",
+ " mp,\n",
+ " simulationtype=\"combined\",\n",
+ " trackingdirection=\"forward\",\n",
+ " weaksinkoption=\"pass_through\",\n",
+ " weaksourceoption=\"pass_through\",\n",
+ " budgetoutputoption=\"summary\",\n",
+ " referencetime=[1, 0, 0],\n",
+ " timepointdata=[30, 2000.0],\n",
+ " zonedataoption=\"on\",\n",
+ " zones=zone_maps,\n",
+ " particlegroups=pgs,\n",
+ " )\n",
+ "\n",
+ " return mp"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "90cdd160",
+ "metadata": {},
+ "source": [
+ "Define a function to build both simulations."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 11,
+ "id": "04051816",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def build_models():\n",
+ " mf6sim = build_mf6_model()\n",
+ " gwf = mf6sim.get_model(nm_mf6)\n",
+ " mp7sim = build_mp7_model(gwf)\n",
+ " return mf6sim, mp7sim"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "dfca56d1",
+ "metadata": {},
+ "source": [
+ "Define a function to run the MODFLOW 6 and MOPATH 7 models/simulations."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 12,
+ "id": "b0933a31",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def run_models(mf6sim, mp7sim):\n",
+ " mf6sim.write_simulation()\n",
+ " success, buff = mf6sim.run_simulation(silent=True, report=True)\n",
+ " for line in buff:\n",
+ " print(line)\n",
+ " assert success, \"Failed to run MODFLOW 6\"\n",
+ "\n",
+ " mp7sim.write_input()\n",
+ " success, buff = mp7sim.run_model(silent=True, report=True)\n",
+ " for line in buff:\n",
+ " print(line)\n",
+ " assert success, \"Failed to run MODPATH 7\""
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "7b46090d",
+ "metadata": {},
+ "source": [
+ "Define a function to determine termination zones."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 13,
+ "id": "c07eae3a",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def get_term_zone_nns(grid):\n",
+ " wel_locs = [w[0:3] for w in wells]\n",
+ " riv_locs = [(0, i, 19) for i in range(20)]\n",
+ " drn_locs = [(drain[0], drain[1], d) for d in range(drain[2][0], drain[2][1])]\n",
+ " wel_nids = grid.get_node(wel_locs)\n",
+ " riv_nids = grid.get_node(riv_locs)\n",
+ " drn_nids = grid.get_node(drn_locs)\n",
+ " return wel_nids, drn_nids, riv_nids"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "2c650d87",
+ "metadata": {},
+ "source": [
+ "Define functions to load pathline data from the budget file created by the PRT model."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 14,
+ "id": "6ef054ad",
+ "metadata": {
+ "lines_to_end_of_cell_marker": 0,
+ "lines_to_next_cell": 1
+ },
+ "outputs": [],
+ "source": [
+ "from flopy.plot.plotutil import to_mp7_pathlines, to_mp7_endpoints\n",
+ "\n",
+ "def load_mf6_pathlines(p):\n",
+ " # read pathlines into dataframe from csv\n",
+ " pls = pd.read_csv(p)\n",
+ "\n",
+ " # assign a unique particle index column incrementing an integer\n",
+ " # for each unique combination of irpt, iprp, imdl, and trelease\n",
+ " id_cols = [\"imdl\", \"iprp\", \"irpt\", \"trelease\"]\n",
+ " pls = pls.sort_values(id_cols)\n",
+ " particles = pls.groupby(id_cols)\n",
+ " pls[\"particleid\"] = particles.ngroup()\n",
+ "\n",
+ " # select endpoints\n",
+ " eps = pls.sort_values(\"t\").groupby(id_cols).tail(1)\n",
+ "\n",
+ " # add a terminating zone column to pathlines dataframe\n",
+ " def set_term_zone(row):\n",
+ " id = row[\"particleid\"]\n",
+ " ep = eps[eps[\"particleid\"] == id].iloc[0]\n",
+ " return ep[\"izone\"]\n",
+ " pls[\"izone_term\"] = pls.apply(lambda r: set_term_zone(r), axis=1)\n",
+ "\n",
+ " # return a dict keyed on capture zone\n",
+ " return {\n",
+ " \"well\": to_mp7_pathlines(pls[(pls[\"izone_term\"] == 2) | (pls[\"izone_term\"] == 3)]),\n",
+ " \"drain\": to_mp7_pathlines(pls[pls[\"izone_term\"] == 4]),\n",
+ " \"river\": to_mp7_pathlines(pls[pls[\"izone_term\"] == 5]),\n",
+ " }\n",
+ "\n",
+ "def load_mf6_endpoints(p):\n",
+ " # read pathlines into dataframe from csv\n",
+ " pls = pd.read_csv(p)\n",
+ "\n",
+ " # assign a unique particle index column incrementing an integer\n",
+ " # for each unique combination of irpt, iprp, imdl, and trelease\n",
+ " id_cols = [\"imdl\", \"iprp\", \"irpt\", \"trelease\"]\n",
+ " pls = pls.sort_values(id_cols)\n",
+ " particles = pls.groupby(id_cols)\n",
+ " pls[\"particleid\"] = particles.ngroup()\n",
+ "\n",
+ " # select endpoints\n",
+ " eps = pls.sort_values(\"t\").groupby(id_cols).tail(1)\n",
+ "\n",
+ " # add a terminating zone column to pathlines dataframe\n",
+ " def set_term_zone(row):\n",
+ " id = row[\"particleid\"]\n",
+ " ep = eps[eps[\"particleid\"] == id].iloc[0]\n",
+ " return ep[\"izone\"]\n",
+ " pls[\"izone_term\"] = pls.apply(lambda r: set_term_zone(r), axis=1)\n",
+ "\n",
+ " # return a dict keyed on capture zone\n",
+ " return {\n",
+ " \"all\": to_mp7_endpoints(pls),\n",
+ " \"well\": to_mp7_endpoints(pls[(pls[\"izone_term\"] == 2) | (pls[\"izone_term\"] == 3)]),\n",
+ " \"drain\": to_mp7_endpoints(pls[pls[\"izone_term\"] == 4]),\n",
+ " \"river\": to_mp7_endpoints(pls[pls[\"izone_term\"] == 5]),\n",
+ " }"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "da89757b",
+ "metadata": {},
+ "source": [
+ "Define functions to load pathline and endpoint data from the MODPATH 7 models' output files."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 15,
+ "id": "3ffaa04c",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from flopy.utils.modpathfile import EndpointFile\n",
+ "\n",
+ "def load_mp7_pathlines(pth, grid):\n",
+ " wel_nids, drn_nids, riv_nids = get_term_zone_nns(grid)\n",
+ " mp7pathlines = {}\n",
+ " p = flopy.utils.PathlineFile(pth)\n",
+ " mp7pathlines[\"well\"] = p.get_destination_pathline_data(wel_nids, to_recarray=True)\n",
+ " mp7pathlines[\"drain\"] = p.get_destination_pathline_data(drn_nids, to_recarray=True)\n",
+ " mp7pathlines[\"river\"] = p.get_destination_pathline_data(riv_nids, to_recarray=True)\n",
+ " return mp7pathlines\n",
+ "\n",
+ "def load_mp7_endpoints(pth, grid):\n",
+ " wel_nids, drn_nids, riv_nids = get_term_zone_nns(grid)\n",
+ " mp7endpoints = {}\n",
+ " p = EndpointFile(pth)\n",
+ " # endpts = pthobj.get_alldata()\n",
+ " mp7endpoints[\"well\"] = p.get_destination_endpoint_data(wel_nids)\n",
+ " mp7endpoints[\"drain\"] = p.get_destination_endpoint_data(drn_nids)\n",
+ " mp7endpoints[\"river\"] = p.get_destination_endpoint_data(riv_nids)\n",
+ " return mp7endpoints"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "29367ff5",
+ "metadata": {},
+ "source": [
+ "Define functions to plot the heads and pathline results, with pathlines colored by capture zone."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 16,
+ "id": "dce13e18",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from flopy.export.vtk import Vtk\n",
+ "import pyvista as pv\n",
+ "\n",
+ "def plot_startpoints(ax, gwf, sp, title):\n",
+ " ax.set_aspect(\"equal\")\n",
+ " ax.set_title(title)\n",
+ " mv = flopy.plot.PlotMapView(model=gwf, ax=ax)\n",
+ " mv.plot_grid(lw=0.5)\n",
+ " mv.plot_bc(\"WEL\", plotAll=True)\n",
+ " mv.plot_bc(\"DRN\")\n",
+ " mv.plot_bc(\"RIV\")\n",
+ " mv.plot_endpoint(sp, label=\"start points\", color=\"yellow\")\n",
+ " ax.set_xlim([1000, 2000])\n",
+ " ax.set_ylim([9000, 10000])\n",
+ "\n",
+ "def plot_endpoints(ax, gwf, ep, title):\n",
+ " ax.set_aspect(\"equal\")\n",
+ " ax.set_title(title)\n",
+ "\n",
+ " mv = flopy.plot.PlotMapView(model=gwf, ax=ax)\n",
+ " mv.plot_grid(lw=0.5)\n",
+ " mv.plot_bc(\"DRN\", alpha=0.2)\n",
+ " mv.plot_bc(\"RIV\", alpha=0.2)\n",
+ " mv.plot_bc(\"WEL\", alpha=0.2, plotAll=True)\n",
+ "\n",
+ " if ep[\"well\"] is not None:\n",
+ " mv.plot_endpoint(ep[\"well\"], label=\"end points\", color=\"red\")\n",
+ " if ep[\"drain\"] is not None:\n",
+ " mv.plot_endpoint(ep[\"drain\"], label=\"end points\", color=\"green\")\n",
+ " if ep[\"river\"] is not None:\n",
+ " mv.plot_endpoint(ep[\"river\"], label=\"end points\", color=\"blue\")\n",
+ "\n",
+ " # inset axes, zoom into particle endpoint clusters, first well..\n",
+ " axins1 = ax.inset_axes([-0.82, 0.25, 0.5, 0.5])\n",
+ " mvins1 = flopy.plot.PlotMapView(model=gwf, ax=axins1)\n",
+ " mvins1.plot_grid(lw=0.5)\n",
+ "\n",
+ " if ep[\"well\"] is not None:\n",
+ " mvins1.plot_endpoint(ep[\"well\"], label=\"end points\", color=\"red\")\n",
+ " if ep[\"drain\"] is not None:\n",
+ " mvins1.plot_endpoint(ep[\"drain\"], label=\"end points\", color=\"green\")\n",
+ " if ep[\"river\"] is not None:\n",
+ " mvins1.plot_endpoint(ep[\"river\"], label=\"end points\", color=\"blue\")\n",
+ "\n",
+ " axins1.set_xlim([2000, 3000])\n",
+ " axins1.set_ylim([3900, 4600])\n",
+ " ax.indicate_inset_zoom(axins1)\n",
+ "\n",
+ " # ..then river\n",
+ " axins2 = ax.inset_axes([1.12, 0.25, 0.5, 0.5])\n",
+ " mvins2 = flopy.plot.PlotMapView(model=gwf, ax=axins2)\n",
+ " mvins2.plot_grid(lw=0.5)\n",
+ "\n",
+ " if ep[\"well\"] is not None:\n",
+ " mvins2.plot_endpoint(ep[\"well\"], label=\"end points\", color=\"red\")\n",
+ " if ep[\"drain\"] is not None:\n",
+ " mvins2.plot_endpoint(ep[\"drain\"], label=\"end points\", color=\"green\")\n",
+ " if ep[\"river\"] is not None:\n",
+ " mvins2.plot_endpoint(ep[\"river\"], label=\"end points\", color=\"blue\")\n",
+ "\n",
+ " axins2.set_xlim([9000, 10000])\n",
+ " axins2.set_ylim([2500, 3900])\n",
+ " ax.indicate_inset_zoom(axins2)\n",
+ "\n",
+ "def plot_pathlines(ax, gwf, hd, pl, title):\n",
+ " ax.set_aspect(\"equal\")\n",
+ " ax.set_title(title)\n",
+ " mv = flopy.plot.PlotMapView(model=gwf, ax=ax)\n",
+ " mv.plot_grid(lw=0.5)\n",
+ " mv.plot_bc(\"WEL\", plotAll=True)\n",
+ " mv.plot_bc(\"DRN\")\n",
+ " mv.plot_bc(\"RIV\")\n",
+ " hd = mv.plot_array(hd, alpha=0.2)\n",
+ " cb = plt.colorbar(hd, shrink=0.5, ax=ax)\n",
+ " cb.set_label(\"Head\")\n",
+ " mv.plot_pathline(\n",
+ " pl['well'],\n",
+ " layer=\"all\",\n",
+ " colors=[\"red\"],\n",
+ " label=\"captured by well\",\n",
+ " )\n",
+ " mv.plot_pathline(\n",
+ " pl['drain'],\n",
+ " layer=\"all\",\n",
+ " colors=[\"green\"],\n",
+ " label=\"captured by drain\",\n",
+ " )\n",
+ " mv.plot_pathline(\n",
+ " pl['river'],\n",
+ " layer=\"all\",\n",
+ " colors=[\"blue\"],\n",
+ " label=\"captured by river\",\n",
+ " )\n",
+ " mv.ax.legend()\n",
+ "\n",
+ " # inset axes, zoom into particle release locations....\n",
+ " axins = ax.inset_axes([-0.82, 0.5, 0.5, 0.5])\n",
+ " mvins = flopy.plot.PlotMapView(model=gwf, ax=axins)\n",
+ " mvins.plot_grid(lw=0.5)\n",
+ " mvins.plot_pathline(\n",
+ " pl['well'],\n",
+ " layer=\"all\",\n",
+ " colors=[\"red\"],\n",
+ " label=\"captured by well\",\n",
+ " )\n",
+ " mvins.plot_pathline(\n",
+ " pl['drain'],\n",
+ " layer=\"all\",\n",
+ " colors=[\"green\"],\n",
+ " label=\"captured by drain\",\n",
+ " )\n",
+ " mvins.plot_pathline(\n",
+ " pl['river'],\n",
+ " layer=\"all\",\n",
+ " colors=[\"blue\"],\n",
+ " label=\"captured by river\",\n",
+ " )\n",
+ " axins.set_xlim([1000, 2000])\n",
+ " axins.set_ylim([8500, 9500])\n",
+ " ax.indicate_inset_zoom(axins)\n",
+ "\n",
+ "def plot_pathlines_xc(ax, gwf, hd, pl, line, title):\n",
+ " # ax.set_aspect(\"equal\")\n",
+ " ax.set_title(title)\n",
+ " mv = flopy.plot.PlotCrossSection(model=gwf, ax=ax, line=line)\n",
+ " mv.plot_grid(lw=0.5)\n",
+ " mv.plot_bc(\"WEL\")\n",
+ " mv.plot_bc(\"DRN\")\n",
+ " mv.plot_bc(\"RIV\")\n",
+ " hd = mv.plot_array(hd, alpha=0.2)\n",
+ " cb = plt.colorbar(hd, shrink=0.5, ax=ax)\n",
+ " cb.set_label(\"Head\")\n",
+ " mv.plot_pathline(\n",
+ " pl['well'],\n",
+ " colors=[\"red\"],\n",
+ " label=\"captured by well\",\n",
+ " )\n",
+ " mv.plot_pathline(\n",
+ " pl['drain'],\n",
+ " colors=[\"green\"],\n",
+ " label=\"captured by drain\",\n",
+ " )\n",
+ " mv.plot_pathline(\n",
+ " pl['river'],\n",
+ " colors=[\"blue\"],\n",
+ " label=\"captured by river\",\n",
+ " )\n",
+ " mv.ax.legend()\n",
+ "\n",
+ "def plot_results_2d(gwf, heads, ep1, ep2, pl1, pl2, title1=\"MODFLOW 6 PRT\", title2=\"MODPATH 7\"):\n",
+ " fig, axes = plt.subplots(ncols=2, nrows=5, figsize=(18, 18))\n",
+ " plt.suptitle(\n",
+ " t=\"Example 3: Structured transient forward-tracking model with multiple release times\",\n",
+ " fontsize=14,\n",
+ " y=0.92,\n",
+ " )\n",
+ "\n",
+ " # map view pathlines\n",
+ " plot_pathlines(axes[0][0], gwf, heads, pl1, f\"{title1} map view\")\n",
+ " plot_pathlines(axes[0][1], gwf, heads, pl2, f\"{title2} map view\")\n",
+ "\n",
+ " # map view end points\n",
+ " plot_endpoints(axes[1][0], gwf, ep1, f\"{title1} end points\")\n",
+ " plot_endpoints(axes[1][1], gwf, ep2, f\"{title2} end points\")\n",
+ "\n",
+ " # cross sections\n",
+ " xc_row = 11\n",
+ " xc_line = {\"row\": xc_row}\n",
+ " plot_pathlines_xc(axes[2][0], gwf, heads, pl1, xc_line, f\"{title1} cross section, row {xc_row}\")\n",
+ " plot_pathlines_xc(axes[2][1], gwf, heads, pl2, xc_line, f\"{title2} cross section, row {xc_row}\")\n",
+ " xc_row = 12\n",
+ " xc_line = {\"row\": xc_row}\n",
+ " plot_pathlines_xc(axes[3][0], gwf, heads, pl1, xc_line, f\"{title1} cross section, row {xc_row}\")\n",
+ " plot_pathlines_xc(axes[3][1], gwf, heads, pl2, xc_line, f\"{title2} cross section, row {xc_row}\")\n",
+ " xc_row = 14\n",
+ " xc_line = {\"row\": xc_row}\n",
+ " plot_pathlines_xc(axes[4][0], gwf, heads, pl1, xc_line, f\"{title1} cross section, row {xc_row}\")\n",
+ " plot_pathlines_xc(axes[4][1], gwf, heads, pl2, xc_line, f\"{title2} cross section, row {xc_row}\")\n",
+ " \n",
+ " plt.show()\n",
+ "\n",
+ "def plot_results_3d(gwf, heads, path1, path2, title1=\"MODFLOW 6 PRT\", title2=\"MODPATH 7\"):\n",
+ " # read the VTK files into PyVista meshes\n",
+ " grid = pv.read(path1.parent / f\"{path1.stem}_000000.vtk\")\n",
+ " mf6_pls = pv.read(path1.parent / f\"{path1.stem}_pathline.vtk\")\n",
+ " mp7_pls = pv.read(path2.parent / f\"{path2.stem}_pathline.vtk\")\n",
+ "\n",
+ " # rotate and scale the plot\n",
+ " axes = pv.Axes(show_actor=True, actor_scale=2.0, line_width=5)\n",
+ " grid.rotate_z(160, point=axes.origin, inplace=True)\n",
+ " mf6_pls.rotate_z(160, point=axes.origin, inplace=True)\n",
+ " mp7_pls.rotate_z(160, point=axes.origin, inplace=True)\n",
+ "\n",
+ " # check grid properties\n",
+ " assert grid.n_cells == gwf.modelgrid.nnodes\n",
+ " print(\"Grid has\", grid.n_cells, \"cells\")\n",
+ " print(\"Grid has\", grid.n_arrays, \"arrays\")\n",
+ "\n",
+ " # set PyVista theme and backend\n",
+ " pv.set_plot_theme(\"document\")\n",
+ " pv.set_jupyter_backend('trame')\n",
+ " # pv.set_jupyter_backend('static')\n",
+ "\n",
+ " # create the plot and add the grid and pathline meshes\n",
+ " p = pv.Plotter(shape=(1, 2))\n",
+ "\n",
+ " # todo: when PRT updated to output capture zone,\n",
+ " # add scalars and use them to color pathline pts\n",
+ "\n",
+ " # add subplot for MODFLOW 6 PRT\n",
+ " p.subplot(0, 0)\n",
+ " p.add_text(title1, font_size=20)\n",
+ " p.add_mesh(grid, opacity=0.05)\n",
+ " p.add_mesh(mf6_pls)\n",
+ "\n",
+ " # add subplot for MODPATH 7\n",
+ " p.subplot(0, 1)\n",
+ " p.add_text(title2, font_size=20)\n",
+ " p.add_mesh(grid, opacity=0.05)\n",
+ " p.add_mesh(mp7_pls)\n",
+ "\n",
+ " # zoom in and show the plot\n",
+ " p.camera.zoom(2.4)\n",
+ " p.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "256ac3d3",
+ "metadata": {},
+ "source": [
+ "Define a function to wrap the entire scenario."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 19,
+ "id": "4ed7763d",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import flopy.utils.binaryfile as bf\n",
+ "\n",
+ "def scenario():\n",
+ " # build models\n",
+ " mf6sim, mp7sim = build_models()\n",
+ " gwf = mf6sim.get_model(nm_mf6)\n",
+ " grid = gwf.modelgrid\n",
+ "\n",
+ " # run models\n",
+ " run_models(mf6sim, mp7sim)\n",
+ "\n",
+ " # load results\n",
+ " hds = flopy.utils.HeadFile(sim_ws / headfile).get_data()\n",
+ " mf6_pl = load_mf6_pathlines(sim_ws / trackcsvfile_prt)\n",
+ " mf6_ep = load_mf6_endpoints(sim_ws / trackcsvfile_prt)\n",
+ " mp7_pl = load_mp7_pathlines(sim_ws / pathlinefile_mp7, grid)\n",
+ " mp7_ep = load_mp7_endpoints(sim_ws / endpointfile_mp7, grid)\n",
+ "\n",
+ " # export pathline results to VTK files\n",
+ " vtk = Vtk(model=gwf, binary=True, vertical_exageration=True, smooth=False)\n",
+ " vtk.add_model(gwf)\n",
+ " vtk.add_pathline_points(mf6_pl[\"well\"] + mf6_pl[\"river\"])\n",
+ " mf6_vtk_path = mf6sim.sim_path / f\"{sim_name}_mf6.vtk\"\n",
+ " vtk.write(mf6_vtk_path)\n",
+ " vtk = Vtk(model=gwf, binary=True, vertical_exageration=True, smooth=False)\n",
+ " vtk.add_model(gwf)\n",
+ " vtk.add_pathline_points([mp7_pl[\"well\"], mp7_pl[\"river\"]])\n",
+ " mp7_vtk_path = mf6sim.sim_path / f\"{sim_name}_mp7.vtk\"\n",
+ " vtk.write(mp7_vtk_path)\n",
+ "\n",
+ " # plot results\n",
+ " plot_results_2d(gwf, hds, mf6_ep, mp7_ep, mf6_pl, mp7_pl)\n",
+ " plot_results_3d(gwf, hds, mf6_vtk_path, mp7_vtk_path)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "831dc470",
+ "metadata": {},
+ "source": [
+ "Run the MODPATH 7 example problem 3 scenario."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 20,
+ "id": "d2c0d39b",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Building MODFLOW 6 model\n",
+ "Building MODPATH 7 model...\n",
+ "writing simulation...\n",
+ " writing simulation name file...\n",
+ " writing simulation tdis package...\n",
+ " writing solution package ims...\n",
+ " writing solution package ems...\n",
+ " writing package mp7-p03.gwfprt...\n",
+ " writing model mp7-p03...\n",
+ " writing model name file...\n",
+ " writing package dis...\n",
+ " writing package ic...\n",
+ " writing package npf...\n",
+ " writing package rcha_0...\n",
+ " writing package sto...\n",
+ " writing package wel_0...\n",
+ " writing package riv_0...\n",
+ "INFORMATION: maxbound in ('gwf6', 'riv', 'dimensions') changed to 21 based on size of stress_period_data\n",
+ " writing package drn_0...\n",
+ "INFORMATION: maxbound in ('gwf6', 'drn', 'dimensions') changed to 11 based on size of stress_period_data\n",
+ " writing package oc...\n",
+ " writing model mp7-p03_prt...\n",
+ " writing model name file...\n",
+ " writing package dis...\n",
+ " writing package mip...\n",
+ " writing package prp1...\n",
+ " writing package oc...\n",
+ " writing package fmi...\n",
+ " MODFLOW 6\n",
+ " U.S. GEOLOGICAL SURVEY MODULAR HYDROLOGIC MODEL\n",
+ " VERSION 6.5.0.dev0+prt (preliminary) 08/18/2023\n",
+ " ***DEVELOP MODE***\n",
+ "\n",
+ " MODFLOW 6 compiled Aug 18 2023 16:34:51 with Intel(R) Fortran Intel(R) 64\n",
+ " Compiler Classic for applications running on Intel(R) 64, Version 2021.7.0\n",
+ " Build 20220726_000000\n",
+ "\n",
+ "This software is preliminary or provisional and is subject to \n",
+ "revision. It is being provided to meet the need for timely best \n",
+ "science. The software has not received final approval by the U.S. \n",
+ "Geological Survey (USGS). No warranty, expressed or implied, is made \n",
+ "by the USGS or the U.S. Government as to the functionality of the \n",
+ "software and related material nor shall the fact of release \n",
+ "constitute any such warranty. The software is provided on the \n",
+ "condition that neither the USGS nor the U.S. Government shall be held \n",
+ "liable for any damages resulting from the authorized or unauthorized \n",
+ "use of the software.\n",
+ "\n",
+ " \n",
+ " Run start date and time (yyyy/mm/dd hh:mm:ss): 2023/09/15 20:27:10\n",
+ " \n",
+ " Writing simulation list file: mfsim.lst\n",
+ " Using Simulation name file: mfsim.nam\n",
+ " \n",
+ " Solving: Stress period: 1 Time step: 1\n",
+ " Solving: Stress period: 2 Time step: 1\n",
+ " Solving: Stress period: 2 Time step: 2\n",
+ " Solving: Stress period: 2 Time step: 3\n",
+ " Solving: Stress period: 2 Time step: 4\n",
+ " Solving: Stress period: 2 Time step: 5\n",
+ " Solving: Stress period: 2 Time step: 6\n",
+ " Solving: Stress period: 2 Time step: 7\n",
+ " Solving: Stress period: 2 Time step: 8\n",
+ " Solving: Stress period: 2 Time step: 9\n",
+ " Solving: Stress period: 2 Time step: 10\n",
+ " Solving: Stress period: 3 Time step: 1\n",
+ " Solving: Stress period: 4 Time step: 1\n",
+ " Solving: Stress period: 4 Time step: 2\n",
+ " Solving: Stress period: 4 Time step: 3\n",
+ " Solving: Stress period: 4 Time step: 4\n",
+ " Solving: Stress period: 4 Time step: 5\n",
+ " Solving: Stress period: 4 Time step: 6\n",
+ " Solving: Stress period: 4 Time step: 7\n",
+ " Solving: Stress period: 4 Time step: 8\n",
+ " Solving: Stress period: 4 Time step: 9\n",
+ " Solving: Stress period: 4 Time step: 10\n",
+ " Solving: Stress period: 5 Time step: 1\n",
+ " \n",
+ " Run end date and time (yyyy/mm/dd hh:mm:ss): 2023/09/15 20:27:11\n",
+ " Elapsed run time: 0.601 Seconds\n",
+ " \n",
+ " Normal termination of simulation.\n",
+ "\n",
+ "MODPATH Version 7.2.001 \n",
+ "Program compiled Jul 05 2023 20:35:59 with IFORT compiler (ver. 20.21.7) \n",
+ " \n",
+ " \n",
+ "Run particle tracking simulation ...\n",
+ "Processing Time Step 1 Period 1. Time = 9.00000E+04 Steady-state flow \n",
+ "Processing Time Step 1 Period 2. Time = 9.00200E+04 Steady-state flow \n",
+ "Processing Time Step 2 Period 2. Time = 9.00400E+04 Steady-state flow \n",
+ "Processing Time Step 3 Period 2. Time = 9.00600E+04 Steady-state flow \n",
+ "Processing Time Step 4 Period 2. Time = 9.00800E+04 Steady-state flow \n",
+ "Processing Time Step 5 Period 2. Time = 9.01000E+04 Steady-state flow \n",
+ "Processing Time Step 6 Period 2. Time = 9.01200E+04 Steady-state flow \n",
+ "Processing Time Step 7 Period 2. Time = 9.01400E+04 Steady-state flow \n",
+ "Processing Time Step 8 Period 2. Time = 9.01600E+04 Steady-state flow \n",
+ "Processing Time Step 9 Period 2. Time = 9.01800E+04 Steady-state flow \n",
+ "Processing Time Step 10 Period 2. Time = 9.02000E+04 Steady-state flow \n",
+ "Processing Time Step 1 Period 3. Time = 1.00000E+05 Steady-state flow \n",
+ "Processing Time Step 1 Period 4. Time = 1.00322E+05 Transient flow \n",
+ "Processing Time Step 2 Period 4. Time = 1.00805E+05 Transient flow \n",
+ "Processing Time Step 3 Period 4. Time = 1.01530E+05 Transient flow \n",
+ "Processing Time Step 4 Period 4. Time = 1.02617E+05 Transient flow \n",
+ "Processing Time Step 5 Period 4. Time = 1.04247E+05 Transient flow \n",
+ "Processing Time Step 6 Period 4. Time = 1.06693E+05 Transient flow \n",
+ "Processing Time Step 7 Period 4. Time = 1.10362E+05 Transient flow \n",
+ "Processing Time Step 8 Period 4. Time = 1.15864E+05 Transient flow \n",
+ "Processing Time Step 9 Period 4. Time = 1.24119E+05 Transient flow \n",
+ "Processing Time Step 10 Period 4. Time = 1.36500E+05 Transient flow \n",
+ "Processing Time Step 1 Period 5. Time = 2.36500E+05 Steady-state flow \n",
+ "\n",
+ "Particle Summary:\n",
+ " 0 particles are pending release.\n",
+ " 0 particles remain active.\n",
+ " 627 particles terminated at boundary faces.\n",
+ " 0 particles terminated at weak sink cells.\n",
+ " 0 particles terminated at weak source cells.\n",
+ " 373 particles terminated at strong source/sink cells.\n",
+ " 0 particles terminated in cells with a specified zone number.\n",
+ " 0 particles were stranded in inactive or dry cells.\n",
+ " 0 particles were unreleased.\n",
+ " 0 particles have an unknown status.\n",
+ " \n",
+ "Normal termination. \n"
+ ]
+ },
+ {
+ "ename": "ValueError",
+ "evalue": "select with an empty condition list is not possible",
+ "output_type": "error",
+ "traceback": [
+ "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
+ "\u001b[1;31mValueError\u001b[0m Traceback (most recent call last)",
+ "\u001b[1;32mc:\\Users\\wbonelli\\dev\\modflow6-examples\\notebooks\\ex-prt-mp7-p03.ipynb Cell 36\u001b[0m line \u001b[0;36m1\n\u001b[1;32m----> 1\u001b[0m scenario()\n",
+ "\u001b[1;32mc:\\Users\\wbonelli\\dev\\modflow6-examples\\notebooks\\ex-prt-mp7-p03.ipynb Cell 36\u001b[0m line \u001b[0;36m1\n\u001b[0;32m 13\u001b[0m hds \u001b[39m=\u001b[39m flopy\u001b[39m.\u001b[39mutils\u001b[39m.\u001b[39mHeadFile(sim_ws \u001b[39m/\u001b[39m headfile)\u001b[39m.\u001b[39mget_data()\n\u001b[0;32m 14\u001b[0m mf6_pl \u001b[39m=\u001b[39m load_mf6_pathlines(sim_ws \u001b[39m/\u001b[39m trackcsvfile_prt)\n\u001b[1;32m---> 15\u001b[0m mf6_ep \u001b[39m=\u001b[39m load_mf6_endpoints(sim_ws \u001b[39m/\u001b[39;49m trackcsvfile_prt)\n\u001b[0;32m 16\u001b[0m mp7_pl \u001b[39m=\u001b[39m load_mp7_pathlines(sim_ws \u001b[39m/\u001b[39m pathlinefile_mp7, grid)\n\u001b[0;32m 17\u001b[0m mp7_ep \u001b[39m=\u001b[39m load_mp7_endpoints(sim_ws \u001b[39m/\u001b[39m endpointfile_mp7, grid)\n",
+ "\u001b[1;32mc:\\Users\\wbonelli\\dev\\modflow6-examples\\notebooks\\ex-prt-mp7-p03.ipynb Cell 36\u001b[0m line \u001b[0;36m5\n\u001b[0;32m 50\u001b[0m pls[\u001b[39m\"\u001b[39m\u001b[39mizone_term\u001b[39m\u001b[39m\"\u001b[39m] \u001b[39m=\u001b[39m pls\u001b[39m.\u001b[39mapply(\u001b[39mlambda\u001b[39;00m r: set_term_zone(r), axis\u001b[39m=\u001b[39m\u001b[39m1\u001b[39m)\n\u001b[0;32m 52\u001b[0m \u001b[39m# return a dict keyed on capture zone\u001b[39;00m\n\u001b[0;32m 53\u001b[0m \u001b[39mreturn\u001b[39;00m {\n\u001b[0;32m 54\u001b[0m \u001b[39m\"\u001b[39m\u001b[39mall\u001b[39m\u001b[39m\"\u001b[39m: to_mp7_endpoints(pls),\n\u001b[1;32m---> 55\u001b[0m \u001b[39m\"\u001b[39m\u001b[39mwell\u001b[39m\u001b[39m\"\u001b[39m: to_mp7_endpoints(pls[(pls[\u001b[39m\"\u001b[39;49m\u001b[39mizone_term\u001b[39;49m\u001b[39m\"\u001b[39;49m] \u001b[39m==\u001b[39;49m \u001b[39m2\u001b[39;49m) \u001b[39m|\u001b[39;49m (pls[\u001b[39m\"\u001b[39;49m\u001b[39mizone_term\u001b[39;49m\u001b[39m\"\u001b[39;49m] \u001b[39m==\u001b[39;49m \u001b[39m3\u001b[39;49m)]),\n\u001b[0;32m 56\u001b[0m \u001b[39m\"\u001b[39m\u001b[39mdrain\u001b[39m\u001b[39m\"\u001b[39m: to_mp7_endpoints(pls[pls[\u001b[39m\"\u001b[39m\u001b[39mizone_term\u001b[39m\u001b[39m\"\u001b[39m] \u001b[39m==\u001b[39m \u001b[39m4\u001b[39m]),\n\u001b[0;32m 57\u001b[0m \u001b[39m\"\u001b[39m\u001b[39mriver\u001b[39m\u001b[39m\"\u001b[39m: to_mp7_endpoints(pls[pls[\u001b[39m\"\u001b[39m\u001b[39mizone_term\u001b[39m\u001b[39m\"\u001b[39m] \u001b[39m==\u001b[39m \u001b[39m5\u001b[39m]),\n\u001b[0;32m 58\u001b[0m }\n",
+ "File \u001b[1;32mc:\\Users\\wbonelli\\AppData\\Local\\miniconda3\\envs\\modflow6-examples\\lib\\site-packages\\flopy\\plot\\plotutil.py:2726\u001b[0m, in \u001b[0;36mto_mp7_endpoints\u001b[1;34m(data)\u001b[0m\n\u001b[0;32m 2722\u001b[0m startlayers \u001b[39m=\u001b[39m startpts[\u001b[39m\"\u001b[39m\u001b[39milay\u001b[39m\u001b[39m\"\u001b[39m]\u001b[39m.\u001b[39mto_numpy()\n\u001b[0;32m 2723\u001b[0m conditions \u001b[39m=\u001b[39m [\n\u001b[0;32m 2724\u001b[0m startpts[seqn_key]\u001b[39m.\u001b[39meq(row[seqn_key]) \u001b[39mfor\u001b[39;00m i, row \u001b[39min\u001b[39;00m startpts\u001b[39m.\u001b[39miterrows()\n\u001b[0;32m 2725\u001b[0m ]\n\u001b[1;32m-> 2726\u001b[0m endpts[\u001b[39m\"\u001b[39m\u001b[39mx0\u001b[39m\u001b[39m\"\u001b[39m] \u001b[39m=\u001b[39m np\u001b[39m.\u001b[39;49mselect(conditions, startxs)\n\u001b[0;32m 2727\u001b[0m endpts[\u001b[39m\"\u001b[39m\u001b[39my0\u001b[39m\u001b[39m\"\u001b[39m] \u001b[39m=\u001b[39m np\u001b[39m.\u001b[39mselect(conditions, startys)\n\u001b[0;32m 2728\u001b[0m endpts[\u001b[39m\"\u001b[39m\u001b[39mz0\u001b[39m\u001b[39m\"\u001b[39m] \u001b[39m=\u001b[39m np\u001b[39m.\u001b[39mselect(conditions, startzs)\n",
+ "File \u001b[1;32mc:\\Users\\wbonelli\\AppData\\Local\\miniconda3\\envs\\modflow6-examples\\lib\\site-packages\\numpy\\lib\\function_base.py:818\u001b[0m, in \u001b[0;36mselect\u001b[1;34m(condlist, choicelist, default)\u001b[0m\n\u001b[0;32m 816\u001b[0m \u001b[39m# Now that the dtype is known, handle the deprecated select([], []) case\u001b[39;00m\n\u001b[0;32m 817\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39mlen\u001b[39m(condlist) \u001b[39m==\u001b[39m \u001b[39m0\u001b[39m:\n\u001b[1;32m--> 818\u001b[0m \u001b[39mraise\u001b[39;00m \u001b[39mValueError\u001b[39;00m(\u001b[39m\"\u001b[39m\u001b[39mselect with an empty condition list is not possible\u001b[39m\u001b[39m\"\u001b[39m)\n\u001b[0;32m 820\u001b[0m choicelist \u001b[39m=\u001b[39m [np\u001b[39m.\u001b[39masarray(choice) \u001b[39mfor\u001b[39;00m choice \u001b[39min\u001b[39;00m choicelist]\n\u001b[0;32m 822\u001b[0m \u001b[39mtry\u001b[39;00m:\n",
+ "\u001b[1;31mValueError\u001b[0m: select with an empty condition list is not possible"
+ ]
+ }
+ ],
+ "source": [
+ "scenario()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "4c399341",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# wdir = Path(\"/Users/wes/dev/modpath-v7/examples/ex03_mf6\")\n",
+ "# \n",
+ "# mf6sim, mp7sim = build_models()\n",
+ "# gwf = mf6sim.get_model(nm_mf6)\n",
+ "# grid = gwf.modelgrid\n",
+ "# \n",
+ "# # load results\n",
+ "# hds = flopy.utils.HeadFile(sim_ws / headfile).get_data()\n",
+ "# cbb = bf.CellBudgetFile(sim_ws / budgetfile_prt)\n",
+ "# \n",
+ "# mf6_ep = load_mf6_endpoints(cbb, gwf)\n",
+ "# mp7_ep = load_mp7_endpoints(sim_ws / endpointfile_mp7, grid)\n",
+ "# mf6_pl = load_mf6_pathlines(cbb, gwf)\n",
+ "# mp7_pl = load_mp7_pathlines(sim_ws / pathlinefile_mp7, grid)\n",
+ "# \n",
+ "# ref_ep = load_mp7_endpoints(wdir / \"ex03a_mf6.endpoint7\", grid)\n",
+ "# ref_pl = load_mp7_pathlines(wdir / \"ex03a_mf6.pathline7\", grid)\n",
+ "# \n",
+ "# plot_results_2d(gwf, hds, ref_ep, mp7_ep, ref_pl, mp7_pl, \"MODPATH 7 (reference)\", \"MODPATH 7\")"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "venv",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.10.12"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/scripts/ex-prt-mp7-p03.py b/scripts/ex-prt-mp7-p03.py
new file mode 100644
index 00000000..642d84cb
--- /dev/null
+++ b/scripts/ex-prt-mp7-p03.py
@@ -0,0 +1,800 @@
+# ---
+# jupyter:
+# jupytext:
+# text_representation:
+# extension: .py
+# format_name: light
+# format_version: '1.5'
+# jupytext_version: 1.14.7
+# kernelspec:
+# display_name: venv
+# language: python
+# name: python3
+# ---
+
+# ## Particle tracking through a transient flow system
+#
+# Application of a MODFLOW 6 particle-tracking (PRT) model to solve example 3 from the MODPATH 7 documentation.
+#
+# This is a transient MODFLOW 6 simulation with a flow system similar to examples 1 and 2.
+# Starting at 90,000 days, a 25 particles are released from 4 cells in the grid's top left corner.
+# This occurs every 20 days for 200 days. 1000 particles are released in total.
+# The model's 2 wells begin to pump at a constant rate 100,000 days into the simulation,
+# and continue to pump at a constant rate for the remainder of the simulation.
+#
+# The original problem specifies three stress periods:
+#
+# | Stress period | Type | Time steps | Length (days) | Multiplier |
+# |:--------------|:-------------|:-----------|:--------------|:-----------|
+# | 1 | steady-state | 1 | 100000 | 1 |
+# | 2 | transient | 10 | 36500 | 1.5 |
+# | 3 | steady-state | 1 | 100000 | 1 |
+#
+# The original solution achieves particle release timing by setting MODPATH 7's reference time option to 2 and value to 0.9. Since PRT does not yet support fractional time steps, we modify the problems' time discretization to 5 stress periods, with a finer time discretization in the 2nd for particle release:
+#
+# | Stress period | Type | Time steps | Length (days) | Multiplier |
+# |:--------------|:-------------|:-----------|:--------------|:-----------|
+# | 1 | steady-state | 1 | 90000 | 1 |
+# | 2 | steady-state | 10 | 200 | 1 |
+# | 3 | steady-state | 1 | 9800 | 1 |
+# | 4 | transient | 10 | 36500 | 1.5 |
+# | 5 | steady-state | 1 | 100000 | 1 |
+#
+
+perioddata = [
+ # perlen, nstp, tsmult
+ (90000, 1, 1),
+ (200, 10, 1),
+ (9800, 1, 1),
+ (36500, 10, 1.5),
+ (100000, 1, 1),
+]
+
+#
+# First import dependencies.
+
+# +
+import os
+import sys
+import matplotlib.pyplot as plt
+import flopy
+import numpy as np
+import pandas as pd
+from pathlib import Path
+
+sys.path.append(os.path.join("..", "common"))
+import config
+# -
+
+# Setup model name and workspace variables.
+
+# +
+ws = config.base_ws
+sim_name = "mp7-p03"
+example_name = "ex-prt-" + sim_name
+sim_ws = Path(ws) / example_name
+
+nm_mf6 = sim_name
+nm_prt = sim_name + "_prt"
+
+headfile = f"{sim_name}.hds"
+budgetfile = f"{sim_name}.cbb"
+budgetfile_prt = f"{nm_prt}.cbb"
+trackcsvfile_prt = f"{nm_prt}.trk.csv"
+pathlinefile_mp7 = f"{sim_name}_mp.mppth"
+endpointfile_mp7 = f"{sim_name}_mp.mpend"
+# -
+
+# Define some shared variables, starting with discretization and other flow model data.
+
+nlay, nrow, ncol = 3, 21, 20
+delr = delc = 500.0
+top = 350.0
+botm = [220.0, 200.0, 0.0]
+laytyp = [1, 0, 0]
+kh = [50.0, 0.01, 200.0]
+kv = [10.0, 0.01, 20.0]
+rchv = 0.005
+riv_h = 320.0
+riv_z = 317.0
+riv_c = 1.0e5
+
+# Define well data. Although this notebook will refer to layer/row/column indices starting at 1, indices in FloPy (and more generally in Python) are zero-based. A negative discharge indicates pumping, while a positive value indicates injection.
+
+wells = [
+ # layer, row, col, discharge
+ (0, 10, 9, -75000),
+ (2, 12, 4, -100000),
+]
+
+# Define the drain location.
+
+drain = (0, 14, (9, 20))
+
+# Configure locations for particle tracking to terminate. We have three explicitly defined termination zones:
+#
+# - 2: the well in layer 1, at row 11, column 10
+# - 3: the well in layer 3, at row 13, column 5
+# - 4: the drain in layer 1, running through row 15 from column 10-20
+# - 5: the river in layer 1, running through column 20
+#
+# MODFLOW 6 reserves zone number 1 to indicate that particles may move freely within the zone.
+
+# +
+zone_maps = []
+
+# zone 1 is the default (non-terminating)
+def fill_zone_1():
+ return np.ones((nrow, ncol), dtype=np.int32)
+
+# zone map for layer 1
+za = fill_zone_1()
+za[wells[0][1:3]] = 2 # well
+za[drain[1], drain[2][0] : drain[2][1]] = 4 # drain
+za[:, ncol - 1] = 5 # river
+zone_maps.append(za)
+
+# zone map for layer 2 (use constant to apply zone 1 to all cells)
+zone_maps.append(1)
+
+# zone map for layer 3
+za = fill_zone_1()
+za[wells[1][1:3]] = 3 # well
+zone_maps.append(za)
+# -
+
+# TODO: plot zone map
+#
+# Define particles to track. Particles are released from the top of a 2x2 square of cells in the upper left of the midel grid's top layer.
+
+# +
+
+rel_minl = rel_maxl = 0
+rel_minr = 2
+rel_maxr = 3
+rel_minc = 2
+rel_maxc = 3
+celldata = flopy.modpath.CellDataType(
+ drape=0,
+ rowcelldivisions=5,
+ columncelldivisions=5,
+ layercelldivisions=1,
+)
+lrcregions = [[rel_minl, rel_minr, rel_minc, rel_maxl, rel_maxr, rel_maxc]]
+lrcpd = flopy.modpath.LRCParticleData(
+ subdivisiondata=celldata,
+ lrcregions=lrcregions,
+)
+pg = flopy.modpath.ParticleGroupLRCTemplate(
+ particlegroupname="PG1",
+ particledata=lrcpd,
+ filename=f"{sim_name}.pg1.sloc",
+ releasedata=(10, 0, 20)
+)
+pgs = [pg]
+defaultiface = {"RECHARGE": 6, "ET": 6}
+# -
+
+# Define a function to create the MODFLOW 6 simulation, which will include a groundwater flow (GWF) model and a particle tracking (PRT) model.
+
+def build_mf6_model():
+ print("Building MODFLOW 6 model")
+
+ # simulation
+ sim = flopy.mf6.MFSimulation(
+ sim_name=sim_name,
+ exe_name="mf6",
+ version="mf6",
+ sim_ws=sim_ws
+ )
+
+ # temporal discretization
+ tdis = flopy.mf6.modflow.mftdis.ModflowTdis(
+ sim,
+ pname="tdis",
+ time_units="DAYS",
+ nper=len(perioddata),
+ perioddata=perioddata
+ )
+
+ # groundwater flow (gwf) model
+ model_nam_file = f"{sim_name}.nam"
+ gwf = flopy.mf6.ModflowGwf(
+ sim, modelname=sim_name, model_nam_file=model_nam_file, save_flows=True
+ )
+
+ # iterative model solver (ims) package
+ ims = flopy.mf6.modflow.mfims.ModflowIms(
+ sim,
+ pname="ims",
+ complexity="SIMPLE",
+ )
+
+ # grid discretization
+ dis = flopy.mf6.modflow.mfgwfdis.ModflowGwfdis(
+ gwf,
+ pname="dis",
+ nlay=nlay,
+ nrow=nrow,
+ ncol=ncol,
+ length_units="FEET",
+ delr=delr,
+ delc=delc,
+ top=top,
+ botm=botm,
+ )
+
+ # initial conditions
+ ic = flopy.mf6.modflow.mfgwfic.ModflowGwfic(gwf, pname="ic", strt=320)
+
+ # node property flow
+ npf = flopy.mf6.modflow.mfgwfnpf.ModflowGwfnpf(
+ gwf, pname="npf", icelltype=laytyp, k=kh, k33=kv, save_flows=True, save_specific_discharge=True, save_saturation=True
+ )
+
+ # recharge
+ rch = flopy.mf6.modflow.mfgwfrcha.ModflowGwfrcha(gwf, recharge=rchv)
+
+ # storage
+ sto = flopy.mf6.modflow.mfgwfsto.ModflowGwfsto(
+ gwf, save_flows=True, iconvert=1, ss=0.0001, sy=0.1,
+ steady_state={0: True, 4: True},
+ transient={3: True},
+ )
+
+ # wells
+ def no_flow(w):
+ return w[0], w[1], w[2], 0
+
+ nf_wells = [no_flow(w) for w in wells]
+ wel = flopy.mf6.modflow.mfgwfwel.ModflowGwfwel(
+ gwf,
+ maxbound=2,
+ stress_period_data={0: nf_wells, 1: nf_wells, 2: nf_wells, 3: wells, 4: wells},
+ )
+
+ # river
+ riv_iface = 6
+ riv_iflowface = -1
+ rd = [[(0, i, ncol - 1), riv_h, riv_c, riv_z, riv_iface, riv_iflowface] for i in range(nrow)]
+ flopy.mf6.modflow.mfgwfriv.ModflowGwfriv(
+ gwf,
+ auxiliary=["iface", "iflowface"],
+ stress_period_data={0: rd}
+ )
+
+ # drain (set auxiliary IFACE var to 6 for top of cell)
+ drn_iface = 6
+ drn_iflowface = -1
+ dd = [
+ [drain[0], drain[1], i + drain[2][0], 322.5, 100000.0, drn_iface, drn_iflowface]
+ for i in range(drain[2][1] - drain[2][0])
+ ]
+ drn = flopy.mf6.modflow.mfgwfdrn.ModflowGwfdrn(
+ gwf,
+ auxiliary=["iface", "iflowface"],
+ stress_period_data={0: dd})
+
+ # output control
+ oc = flopy.mf6.modflow.mfgwfoc.ModflowGwfoc(
+ gwf,
+ pname="oc",
+ saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")],
+ head_filerecord=[headfile],
+ budget_filerecord=[budgetfile],
+ )
+
+ # Instantiate the MODFLOW 6 prt model
+ prt = flopy.mf6.ModflowPrt(
+ sim, modelname=nm_prt, model_nam_file="{}.nam".format(nm_prt)
+ )
+
+ # Instantiate the MODFLOW 6 prt discretization package
+ flopy.mf6.modflow.mfgwfdis.ModflowGwfdis(
+ prt, pname="dis",
+ nlay=nlay, nrow=nrow, ncol=ncol,
+ length_units="FEET",
+ delr=delr, delc=delc,
+ top=top, botm=botm,
+ )
+
+ # Instantiate the MODFLOW 6 prt model input package
+ porosity = 0.1
+ flopy.mf6.ModflowPrtmip(prt, pname="mip", porosity=porosity)
+
+ releasepts = list(lrcpd.to_prp(gwf.modelgrid))
+ flopy.mf6.ModflowPrtprp(
+ prt, pname="prp1", filename="{}_1.prp".format(nm_prt),
+ nreleasepts=len(releasepts), packagedata=releasepts,
+ perioddata={
+ 0: [],
+ 1: ["STEPS"] + [str(i + 1) for i in range(10)],
+ 2: [],
+ 3: [],
+ 4: []
+ },
+ )
+
+ # Instantiate the MODFLOW 6 prt output control package
+ flopy.mf6.ModflowPrtoc(
+ prt,
+ pname="oc",
+ budget_filerecord=[budgetfile_prt],
+ trackcsv_filerecord=[trackcsvfile_prt],
+ saverecord=[("BUDGET", "ALL")],
+ )
+
+ # Instantiate the MODFLOW 6 prt flow model interface
+ flopy.mf6.ModflowPrtfmi(prt, packagedata=[
+ ("GWFHEAD", headfile),
+ ("GWFBUDGET", budgetfile)
+ ])
+
+ # Create the MODFLOW 6 gwf-prt model exchange
+ flopy.mf6.ModflowGwfprt(
+ sim, exgtype="GWF6-PRT6",
+ exgmnamea=nm_mf6, exgmnameb=nm_prt,
+ filename="{}.gwfprt".format(nm_mf6),
+ )
+
+ # Create an explicit model solution (EMS) for the MODFLOW 6 prt model
+ ems = flopy.mf6.ModflowEms(
+ sim, pname="ems",
+ filename=f"{nm_prt}.ems",
+ )
+ sim.register_solution_package(ems, [prt.name])
+
+ return sim
+
+
+# Define a function create the MODPATH 7 model & simulation.
+
+def build_mp7_model(gwf):
+ print("Building MODPATH 7 model...")
+
+ mp = flopy.modpath.Modpath7(
+ modelname=f"{sim_name}_mp",
+ flowmodel=gwf,
+ exe_name="mp7",
+ model_ws=sim_ws,
+ )
+ mpbas = flopy.modpath.Modpath7Bas(mp, porosity=0.1, defaultiface=defaultiface)
+ mpsim = flopy.modpath.Modpath7Sim(
+ mp,
+ simulationtype="combined",
+ trackingdirection="forward",
+ weaksinkoption="pass_through",
+ weaksourceoption="pass_through",
+ budgetoutputoption="summary",
+ referencetime=[1, 0, 0],
+ timepointdata=[30, 2000.0],
+ zonedataoption="on",
+ zones=zone_maps,
+ particlegroups=pgs,
+ )
+
+ return mp
+
+
+# Define a function to build both simulations.
+
+def build_models():
+ mf6sim = build_mf6_model()
+ gwf = mf6sim.get_model(nm_mf6)
+ mp7sim = build_mp7_model(gwf)
+ return mf6sim, mp7sim
+
+
+# Define a function to run the MODFLOW 6 and MOPATH 7 models/simulations.
+
+def run_models(mf6sim, mp7sim):
+ mf6sim.write_simulation()
+ success, buff = mf6sim.run_simulation(silent=True, report=True)
+ for line in buff:
+ print(line)
+ assert success, "Failed to run MODFLOW 6"
+
+ mp7sim.write_input()
+ success, buff = mp7sim.run_model(silent=True, report=True)
+ for line in buff:
+ print(line)
+ assert success, "Failed to run MODPATH 7"
+
+
+# Define a function to determine termination zones.
+
+def get_term_zone_nns(grid):
+ wel_locs = [w[0:3] for w in wells]
+ riv_locs = [(0, i, 19) for i in range(20)]
+ drn_locs = [(drain[0], drain[1], d) for d in range(drain[2][0], drain[2][1])]
+ wel_nids = grid.get_node(wel_locs)
+ riv_nids = grid.get_node(riv_locs)
+ drn_nids = grid.get_node(drn_locs)
+ return wel_nids, drn_nids, riv_nids
+
+
+# Define functions to load pathline data from the budget file created by the PRT model.
+
+# +
+from flopy.plot.plotutil import to_mp7_pathlines, to_mp7_endpoints
+
+def load_mf6_pathlines(p):
+ # read pathlines into dataframe from csv
+ pls = pd.read_csv(p)
+
+ # assign a unique particle index column incrementing an integer
+ # for each unique combination of irpt, iprp, imdl, and trelease
+ id_cols = ["imdl", "iprp", "irpt", "trelease"]
+ pls = pls.sort_values(id_cols)
+ particles = pls.groupby(id_cols)
+ pls["particleid"] = particles.ngroup()
+
+ # select endpoints
+ eps = pls.sort_values("t").groupby(id_cols).tail(1)
+
+ # add a terminating zone column to pathlines dataframe
+ def set_term_zone(row):
+ id = row["particleid"]
+ ep = eps[eps["particleid"] == id].iloc[0]
+ return ep["izone"]
+ pls["izone_term"] = pls.apply(lambda r: set_term_zone(r), axis=1)
+
+ # return a dict keyed on capture zone
+ return {
+ "well": to_mp7_pathlines(pls[(pls["izone_term"] == 2) | (pls["izone_term"] == 3)]),
+ "drain": to_mp7_pathlines(pls[pls["izone_term"] == 4]),
+ "river": to_mp7_pathlines(pls[pls["izone_term"] == 5]),
+ }
+
+def load_mf6_endpoints(p):
+ # read pathlines into dataframe from csv
+ pls = pd.read_csv(p)
+
+ # assign a unique particle index column incrementing an integer
+ # for each unique combination of irpt, iprp, imdl, and trelease
+ id_cols = ["imdl", "iprp", "irpt", "trelease"]
+ pls = pls.sort_values(id_cols)
+ particles = pls.groupby(id_cols)
+ pls["particleid"] = particles.ngroup()
+
+ # select endpoints
+ eps = pls.sort_values("t").groupby(id_cols).tail(1)
+
+ # add a terminating zone column to pathlines dataframe
+ def set_term_zone(row):
+ id = row["particleid"]
+ ep = eps[eps["particleid"] == id].iloc[0]
+ return ep["izone"]
+ pls["izone_term"] = pls.apply(lambda r: set_term_zone(r), axis=1)
+
+ # return a dict keyed on capture zone
+ return {
+ "all": to_mp7_endpoints(pls),
+ "well": to_mp7_endpoints(pls[(pls["izone_term"] == 2) | (pls["izone_term"] == 3)]),
+ "drain": to_mp7_endpoints(pls[pls["izone_term"] == 4]),
+ "river": to_mp7_endpoints(pls[pls["izone_term"] == 5]),
+ }
+# -
+
+# Define functions to load pathline and endpoint data from the MODPATH 7 models' output files.
+
+# +
+from flopy.utils.modpathfile import EndpointFile
+
+def load_mp7_pathlines(pth, grid):
+ wel_nids, drn_nids, riv_nids = get_term_zone_nns(grid)
+ mp7pathlines = {}
+ p = flopy.utils.PathlineFile(pth)
+ mp7pathlines["well"] = p.get_destination_pathline_data(wel_nids, to_recarray=True)
+ mp7pathlines["drain"] = p.get_destination_pathline_data(drn_nids, to_recarray=True)
+ mp7pathlines["river"] = p.get_destination_pathline_data(riv_nids, to_recarray=True)
+ return mp7pathlines
+
+def load_mp7_endpoints(pth, grid):
+ wel_nids, drn_nids, riv_nids = get_term_zone_nns(grid)
+ mp7endpoints = {}
+ p = EndpointFile(pth)
+ # endpts = pthobj.get_alldata()
+ mp7endpoints["well"] = p.get_destination_endpoint_data(wel_nids)
+ mp7endpoints["drain"] = p.get_destination_endpoint_data(drn_nids)
+ mp7endpoints["river"] = p.get_destination_endpoint_data(riv_nids)
+ return mp7endpoints
+
+
+# -
+
+# Define functions to plot the heads and pathline results, with pathlines colored by capture zone.
+
+# +
+from flopy.export.vtk import Vtk
+import pyvista as pv
+
+def plot_startpoints(ax, gwf, sp, title):
+ ax.set_aspect("equal")
+ ax.set_title(title)
+ mv = flopy.plot.PlotMapView(model=gwf, ax=ax)
+ mv.plot_grid(lw=0.5)
+ mv.plot_bc("WEL", plotAll=True)
+ mv.plot_bc("DRN")
+ mv.plot_bc("RIV")
+ mv.plot_endpoint(sp, label="start points", color="yellow")
+ ax.set_xlim([1000, 2000])
+ ax.set_ylim([9000, 10000])
+
+def plot_endpoints(ax, gwf, ep, title):
+ ax.set_aspect("equal")
+ ax.set_title(title)
+
+ mv = flopy.plot.PlotMapView(model=gwf, ax=ax)
+ mv.plot_grid(lw=0.5)
+ mv.plot_bc("DRN", alpha=0.2)
+ mv.plot_bc("RIV", alpha=0.2)
+ mv.plot_bc("WEL", alpha=0.2, plotAll=True)
+
+ if ep["well"] is not None:
+ mv.plot_endpoint(ep["well"], label="end points", color="red")
+ if ep["drain"] is not None:
+ mv.plot_endpoint(ep["drain"], label="end points", color="green")
+ if ep["river"] is not None:
+ mv.plot_endpoint(ep["river"], label="end points", color="blue")
+
+ # inset axes, zoom into particle endpoint clusters, first well..
+ axins1 = ax.inset_axes([-0.82, 0.25, 0.5, 0.5])
+ mvins1 = flopy.plot.PlotMapView(model=gwf, ax=axins1)
+ mvins1.plot_grid(lw=0.5)
+
+ if ep["well"] is not None:
+ mvins1.plot_endpoint(ep["well"], label="end points", color="red")
+ if ep["drain"] is not None:
+ mvins1.plot_endpoint(ep["drain"], label="end points", color="green")
+ if ep["river"] is not None:
+ mvins1.plot_endpoint(ep["river"], label="end points", color="blue")
+
+ axins1.set_xlim([2000, 3000])
+ axins1.set_ylim([3900, 4600])
+ ax.indicate_inset_zoom(axins1)
+
+ # ..then river
+ axins2 = ax.inset_axes([1.12, 0.25, 0.5, 0.5])
+ mvins2 = flopy.plot.PlotMapView(model=gwf, ax=axins2)
+ mvins2.plot_grid(lw=0.5)
+
+ if ep["well"] is not None:
+ mvins2.plot_endpoint(ep["well"], label="end points", color="red")
+ if ep["drain"] is not None:
+ mvins2.plot_endpoint(ep["drain"], label="end points", color="green")
+ if ep["river"] is not None:
+ mvins2.plot_endpoint(ep["river"], label="end points", color="blue")
+
+ axins2.set_xlim([9000, 10000])
+ axins2.set_ylim([2500, 3900])
+ ax.indicate_inset_zoom(axins2)
+
+def plot_pathlines(ax, gwf, hd, pl, title):
+ ax.set_aspect("equal")
+ ax.set_title(title)
+ mv = flopy.plot.PlotMapView(model=gwf, ax=ax)
+ mv.plot_grid(lw=0.5)
+ mv.plot_bc("WEL", plotAll=True)
+ mv.plot_bc("DRN")
+ mv.plot_bc("RIV")
+ hd = mv.plot_array(hd, alpha=0.2)
+ cb = plt.colorbar(hd, shrink=0.5, ax=ax)
+ cb.set_label("Head")
+ mv.plot_pathline(
+ pl['well'],
+ layer="all",
+ colors=["red"],
+ label="captured by well",
+ )
+ mv.plot_pathline(
+ pl['drain'],
+ layer="all",
+ colors=["green"],
+ label="captured by drain",
+ )
+ mv.plot_pathline(
+ pl['river'],
+ layer="all",
+ colors=["blue"],
+ label="captured by river",
+ )
+ mv.ax.legend()
+
+ # inset axes, zoom into particle release locations....
+ axins = ax.inset_axes([-0.82, 0.5, 0.5, 0.5])
+ mvins = flopy.plot.PlotMapView(model=gwf, ax=axins)
+ mvins.plot_grid(lw=0.5)
+ mvins.plot_pathline(
+ pl['well'],
+ layer="all",
+ colors=["red"],
+ label="captured by well",
+ )
+ mvins.plot_pathline(
+ pl['drain'],
+ layer="all",
+ colors=["green"],
+ label="captured by drain",
+ )
+ mvins.plot_pathline(
+ pl['river'],
+ layer="all",
+ colors=["blue"],
+ label="captured by river",
+ )
+ axins.set_xlim([1000, 2000])
+ axins.set_ylim([8500, 9500])
+ ax.indicate_inset_zoom(axins)
+
+def plot_pathlines_xc(ax, gwf, hd, pl, line, title):
+ # ax.set_aspect("equal")
+ ax.set_title(title)
+ mv = flopy.plot.PlotCrossSection(model=gwf, ax=ax, line=line)
+ mv.plot_grid(lw=0.5)
+ mv.plot_bc("WEL")
+ mv.plot_bc("DRN")
+ mv.plot_bc("RIV")
+ hd = mv.plot_array(hd, alpha=0.2)
+ cb = plt.colorbar(hd, shrink=0.5, ax=ax)
+ cb.set_label("Head")
+ mv.plot_pathline(
+ pl['well'],
+ colors=["red"],
+ label="captured by well",
+ )
+ mv.plot_pathline(
+ pl['drain'],
+ colors=["green"],
+ label="captured by drain",
+ )
+ mv.plot_pathline(
+ pl['river'],
+ colors=["blue"],
+ label="captured by river",
+ )
+ mv.ax.legend()
+
+def plot_results_2d(gwf, heads, ep1, ep2, pl1, pl2, title1="MODFLOW 6 PRT", title2="MODPATH 7"):
+ fig, axes = plt.subplots(ncols=2, nrows=5, figsize=(18, 18))
+ plt.suptitle(
+ t="Example 3: Structured transient forward-tracking model with multiple release times",
+ fontsize=14,
+ y=0.92,
+ )
+
+ # map view pathlines
+ plot_pathlines(axes[0][0], gwf, heads, pl1, f"{title1} map view")
+ plot_pathlines(axes[0][1], gwf, heads, pl2, f"{title2} map view")
+
+ # map view end points
+ plot_endpoints(axes[1][0], gwf, ep1, f"{title1} end points")
+ plot_endpoints(axes[1][1], gwf, ep2, f"{title2} end points")
+
+ # cross sections
+ xc_row = 11
+ xc_line = {"row": xc_row}
+ plot_pathlines_xc(axes[2][0], gwf, heads, pl1, xc_line, f"{title1} cross section, row {xc_row}")
+ plot_pathlines_xc(axes[2][1], gwf, heads, pl2, xc_line, f"{title2} cross section, row {xc_row}")
+ xc_row = 12
+ xc_line = {"row": xc_row}
+ plot_pathlines_xc(axes[3][0], gwf, heads, pl1, xc_line, f"{title1} cross section, row {xc_row}")
+ plot_pathlines_xc(axes[3][1], gwf, heads, pl2, xc_line, f"{title2} cross section, row {xc_row}")
+ xc_row = 14
+ xc_line = {"row": xc_row}
+ plot_pathlines_xc(axes[4][0], gwf, heads, pl1, xc_line, f"{title1} cross section, row {xc_row}")
+ plot_pathlines_xc(axes[4][1], gwf, heads, pl2, xc_line, f"{title2} cross section, row {xc_row}")
+
+ plt.show()
+
+def plot_results_3d(gwf, heads, path1, path2, title1="MODFLOW 6 PRT", title2="MODPATH 7"):
+ # read the VTK files into PyVista meshes
+ grid = pv.read(path1.parent / f"{path1.stem}_000000.vtk")
+ mf6_pls = pv.read(path1.parent / f"{path1.stem}_pathline.vtk")
+ mp7_pls = pv.read(path2.parent / f"{path2.stem}_pathline.vtk")
+
+ # rotate and scale the plot
+ axes = pv.Axes(show_actor=True, actor_scale=2.0, line_width=5)
+ grid.rotate_z(160, point=axes.origin, inplace=True)
+ mf6_pls.rotate_z(160, point=axes.origin, inplace=True)
+ mp7_pls.rotate_z(160, point=axes.origin, inplace=True)
+
+ # check grid properties
+ assert grid.n_cells == gwf.modelgrid.nnodes
+ print("Grid has", grid.n_cells, "cells")
+ print("Grid has", grid.n_arrays, "arrays")
+
+ # set PyVista theme and backend
+ pv.set_plot_theme("document")
+ pv.set_jupyter_backend('trame')
+ # pv.set_jupyter_backend('static')
+
+ # create the plot and add the grid and pathline meshes
+ p = pv.Plotter(shape=(1, 2))
+
+ # todo: when PRT updated to output capture zone,
+ # add scalars and use them to color pathline pts
+
+ # add subplot for MODFLOW 6 PRT
+ p.subplot(0, 0)
+ p.add_text(title1, font_size=20)
+ p.add_mesh(grid, opacity=0.05)
+ p.add_mesh(mf6_pls)
+
+ # add subplot for MODPATH 7
+ p.subplot(0, 1)
+ p.add_text(title2, font_size=20)
+ p.add_mesh(grid, opacity=0.05)
+ p.add_mesh(mp7_pls)
+
+ # zoom in and show the plot
+ p.camera.zoom(2.4)
+ p.show()
+
+
+# -
+
+# Define a function to wrap the entire scenario.
+
+# +
+import flopy.utils.binaryfile as bf
+
+def scenario():
+ # build models
+ mf6sim, mp7sim = build_models()
+ gwf = mf6sim.get_model(nm_mf6)
+ grid = gwf.modelgrid
+
+ # run models
+ run_models(mf6sim, mp7sim)
+
+ # load results
+ hds = flopy.utils.HeadFile(sim_ws / headfile).get_data()
+ mf6_pl = load_mf6_pathlines(sim_ws / trackcsvfile_prt)
+ mf6_ep = load_mf6_endpoints(sim_ws / trackcsvfile_prt)
+ mp7_pl = load_mp7_pathlines(sim_ws / pathlinefile_mp7, grid)
+ mp7_ep = load_mp7_endpoints(sim_ws / endpointfile_mp7, grid)
+
+ # export pathline results to VTK files
+ vtk = Vtk(model=gwf, binary=True, vertical_exageration=True, smooth=False)
+ vtk.add_model(gwf)
+ vtk.add_pathline_points(mf6_pl["well"] + mf6_pl["river"])
+ mf6_vtk_path = mf6sim.sim_path / f"{sim_name}_mf6.vtk"
+ vtk.write(mf6_vtk_path)
+ vtk = Vtk(model=gwf, binary=True, vertical_exageration=True, smooth=False)
+ vtk.add_model(gwf)
+ vtk.add_pathline_points([mp7_pl["well"], mp7_pl["river"]])
+ mp7_vtk_path = mf6sim.sim_path / f"{sim_name}_mp7.vtk"
+ vtk.write(mp7_vtk_path)
+
+ # plot results
+ plot_results_2d(gwf, hds, mf6_ep, mp7_ep, mf6_pl, mp7_pl)
+ plot_results_3d(gwf, hds, mf6_vtk_path, mp7_vtk_path)
+
+
+# -
+
+# Run the MODPATH 7 example problem 3 scenario.
+
+scenario()
+
+# +
+# wdir = Path("/Users/wes/dev/modpath-v7/examples/ex03_mf6")
+#
+# mf6sim, mp7sim = build_models()
+# gwf = mf6sim.get_model(nm_mf6)
+# grid = gwf.modelgrid
+#
+# # load results
+# hds = flopy.utils.HeadFile(sim_ws / headfile).get_data()
+# cbb = bf.CellBudgetFile(sim_ws / budgetfile_prt)
+#
+# mf6_ep = load_mf6_endpoints(cbb, gwf)
+# mp7_ep = load_mp7_endpoints(sim_ws / endpointfile_mp7, grid)
+# mf6_pl = load_mf6_pathlines(cbb, gwf)
+# mp7_pl = load_mp7_pathlines(sim_ws / pathlinefile_mp7, grid)
+#
+# ref_ep = load_mp7_endpoints(wdir / "ex03a_mf6.endpoint7", grid)
+# ref_pl = load_mp7_pathlines(wdir / "ex03a_mf6.pathline7", grid)
+#
+# plot_results_2d(gwf, hds, ref_ep, mp7_ep, ref_pl, mp7_pl, "MODPATH 7 (reference)", "MODPATH 7")