10a: Particle tracking with MODFLOW 6 PRT
In this exercise, we will use the MODFLOW 6 Particle Tracking (PRT) Model to simulate advective transport with a quadtree version of the Freyberg flow model.
The PRT Model calculates three-dimensional, advective particle trajectories in flowing groundwater. The PRT Model is designed to work with the MODFLOW 6 Groundwater Flow (GWF) Model and uses the same spatial discretization, which may be represented using either a structured (DIS) or an unstructured (DISV) grid. The PRT Model replicates much of the functionality of MODPATH 7 and offers support for a much broader class of unstructured grids. The PRT Model can be run in the same simulation as the associated GWF Model or in a separate simulation that reads previously calculated flows from a binary budget file. Currently, the PRT Model documented does not support grids of DISU type, tracking of particles through advanced stress package features such as lakes or streams reaches, or exchange of particles between PRT models.
This exercise demonstrates setting up and running the PRT Model in a separate simulation contained in a subfolder of the groundwater flow model workspace. We will also do some basic post-processing- making plots and exporting results to a GeoPackage for visualization in a GIS environment.
[1]:
from IPython.display import clear_output, display
import os
from pathlib import Path
import numpy as np
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
from flopy.utils.gridintersect import GridIntersect
import flopy
The location of the contamination patch and the nodes that the define bounding cells of the patch are calculated below.
[2]:
# patch upper left and lower right
xmin, xmax = 250. * 1, 250. * 3
ymin, ymax = (40 - 14) * 250., (40 - 11) * 250.
csx, csy = [xmin, xmin, xmax, xmax, xmin], [ymin, ymax, ymax, ymin, ymin]
polygon = [list(zip(csx, csy))]
(xmin, ymax), (xmax, ymin)
[2]:
((250.0, 7250.0), (750.0, 6500.0))
Define the workspace, model names and key options.
[3]:
load_ws = Path('data/quadtree')
gwf_ws = Path("temp/ex10a")
gwf_name = "project"
#name_mp = f"{name}_mp"
exe_name = 'mf6'
# flow model output files to use with PRT
headfile = f"{gwf_name}.hds"
budgetfile = f"{gwf_name}.cbc"
prt_name = f"{gwf_name}-prt"
prt_model_ws = gwf_ws / 'prt'
# PRT output files
budgetfile_prt = f"{prt_name}.cbc"
trackfile_prt = f"{prt_name}.trk"
trackcsvfile_prt = f"{prt_name}.trk.csv"
# if using "local_z" option
# 1=start particles at top of cell
# 0=start particles at bottom of cell
# 1 resulted in Error: release point (z=29.899951638778955) is above grid top 29.886768340000000
particle_release_zrpt = .99
Load the MODFLOW 6 Model
Load a simulation object using flopy.mf6.MFSimulation().load().
[4]:
%%capture
sim = flopy.mf6.MFSimulation.load(sim_name=gwf_name, exe_name=exe_name,
sim_ws=load_ws)
Load the groundwater flow model
[5]:
gwf = sim.get_model(gwf_name)
gwf.modelgrid
[5]:
xll:0.0; yll:0.0; rotation:0.0; units:meters; lenuni:2
Change the workspace
[6]:
sim.set_sim_path(gwf_ws)
Add save_saturation and save_flows to the NPF Package
(required for PRT)
[7]:
gwf.npf.save_saturation = True
gwf.npf.save_flows = True
Write the model files
[8]:
%%capture
sim.write_simulation()
Run the simulation.
[9]:
sim.run_simulation()
FloPy is using the following executable to run the model: ../../../../../../../../../.local/bin/mf6
MODFLOW 6
U.S. GEOLOGICAL SURVEY MODULAR HYDROLOGIC MODEL
VERSION 6.6.2 05/12/2025
MODFLOW 6 compiled May 24 2025 11:41:46 with Intel(R) Fortran Intel(R) 64
Compiler Classic for applications running on Intel(R) 64, Version 2021.7.0
Build 20220726_000000
This software has been approved for release by the U.S. Geological
Survey (USGS). Although the software has been subjected to rigorous
review, the USGS reserves the right to update the software as needed
pursuant to further analysis and review. No warranty, expressed or
implied, is made by the USGS or the U.S. Government as to the
functionality of the software and related material nor shall the
fact of release constitute any such warranty. Furthermore, the
software is released on condition that neither the USGS nor the U.S.
Government shall be held liable for any damages resulting from its
authorized or unauthorized use. Also refer to the USGS Water
Resources Software User Rights Notice for complete use, copyright,
and distribution information.
MODFLOW runs in SEQUENTIAL mode
Run start date and time (yyyy/mm/dd hh:mm:ss): 2025/09/26 14:26:11
Writing simulation list file: mfsim.lst
Using Simulation name file: mfsim.nam
Solving: Stress period: 1 Time step: 1
Run end date and time (yyyy/mm/dd hh:mm:ss): 2025/09/26 14:26:11
Elapsed run time: 0.396 Seconds
WARNING REPORT:
1. OPTIONS BLOCK VARIABLE 'UNIT_CONVERSION' IN FILE 'project.sfr' WAS
DEPRECATED IN VERSION 6.4.2. SETTING UNIT_CONVERSION DIRECTLY.
Normal termination of simulation.
[9]:
(True, [])
Create and Run the PRT model
Plot the model grid and the location of the contamination patch
[10]:
fig, ax = plt.subplots(figsize=(5, 9))
mm = flopy.plot.PlotMapView(gwf, layer=0, ax=ax)
mm.plot_bc('SFR', color="b", plotAll=True)
mm.plot_bc('WEL', plotAll=True)
mm.plot_inactive(alpha=0.75)
mm.plot_grid(lw=0.25, color='grey')
ax.fill(csx, csy, color='#e534eb');
Get the model cells intersecting the contamination patch
The GridIntersect utility in Flopy has an intersect method that can find the model cells intersecting a set of points, lines, or polygons. Since this is a DISV grid, these will be cell2d values (that are the same across the model layers). Similarly, for a structured grid, GridIntersect would return row, column locations.
[11]:
gx = GridIntersect(gwf.modelgrid)
results = gx.intersect(polygon, 'Polygon')
particle_start_nodes = results.cellids.astype(int)
particle_start_nodes
[11]:
array([1352, 1353, 1459, 1460, 1581, 1582])
Instantiate the MODFLOW 6 prt model
and discretization packages
[12]:
prt_sim = flopy.mf6.MFSimulation(sim_name=prt_name, sim_ws=prt_model_ws)
# Instantiate the MODFLOW 6 temporal discretization package
flopy.mf6.modflow.mftdis.ModflowTdis(
prt_sim,
time_units="DAYS",
nper=1,
perioddata=[(1, 1, 1)], # perlen, nstp, tsmult
)
prt = flopy.mf6.ModflowPrt(
prt_sim, modelname=prt_name, model_nam_file=f"{prt_name}.nam"
)
# Instantiate the MODFLOW 6 prt discretization package
nlay, ncells_per_layer = gwf.dis.botm.array.shape
disv = flopy.mf6.ModflowGwfdisv(
prt,
nlay=nlay,
ncpl=gwf.dis.ncpl.array,
nvert=gwf.dis.nvert.array,
length_units=gwf.dis.length_units.array,
top=gwf.dis.top.array,
botm=gwf.dis.botm.array,
vertices=gwf.dis.vertices.array,
cell2d=gwf.dis.cell2d.array,
idomain=gwf.dis.idomain.array,
xorigin=gwf.dis.xorigin.array,
yorigin=gwf.dis.yorigin.array,
)
Instantiate the MODFLOW 6 PRT Model Input Package.
First make an ``izone`` array with a zone number for each cell. The “izone” for the model cell will be reported in the PRT output, allowing us to more easily track where particles go, and where they ultimately discharge. For example, if we start particles at the water table, we can determine the contributing area for a stream or other boundary condition, by assigning an izone to those cells. One of the izone numbers can be designated as an istopzone; particles entering this zone will be
terminated regardless of whether the zone is a strong sink or not.
In this example, we’ll assign different izone numbers to SFR and Well Package cells.
[13]:
# start with a default zone of 0
izone_array = np.zeros((nlay, ncells_per_layer), dtype=int)
# get the locations of SFR cells
sfr_k, sfr_cellid = zip(*gwf.sfr.packagedata.array['cellid'])
sfr_k[:10], sfr_cellid[:10]
izones = {
'wel': 1,
'sfr': 2
}
Assign SFR cells to zone 2
[14]:
izone_array[sfr_k, sfr_cellid] = izones['sfr']
Assign wells to zone 1 Using the stress period data. Note that in models with more than one stress period, different wells may be represented in different stress periods.
[15]:
well_k = list()
well_cellid = list()
for per, recarray in gwf.wel.stress_period_data.data.items():
well_k_per, well_cellid_per = zip(*recarray['cellid'])
well_k.append(well_k_per)
well_cellid.append(well_cellid_per)
izone_array[well_k, well_cellid] = izones['wel']
Make the package
[16]:
flopy.mf6.ModflowPrtmip(prt,
porosity=0.1,
izone=izone_array
)
[16]:
package_name = mip
filename = project-prt.mip
package_type = mip
model_or_simulation_package = model
model_name = project-prt
Block griddata
--------------------
porosity
{constant 0.1}
izone
Layer_1{internal}
(array([0, 0, 0, ..., 0, 0, 0], shape=(5108,), dtype=int32))
Layer_2{internal}
(array([0, 0, 0, ..., 0, 0, 0], shape=(5108,), dtype=int32))
Layer_3{internal}
(array([0, 0, 0, ..., 0, 0, 0], shape=(5108,), dtype=int32))
Create the particle start data
From the nodes intersecting the polygon we made above, develop a DataFrame of particle start locations. For simplicity, we are just starting a single particle in each cell, but since the x and y are in local (model) coordinates, it would be fairly straightforward to start multiple particles in each cell on a finer spacing.
[17]:
# release particles at each row, column location
prt_start_x, prt_start_y = gwf.modelgrid.get_local_coords(
gwf.modelgrid.xcellcenters[particle_start_nodes],
gwf.modelgrid.ycellcenters[particle_start_nodes])
prt_start_data = pd.DataFrame({
'irptno': np.arange(len(prt_start_x.ravel())),
'k': 0,
'cell2d': particle_start_nodes,
'xrpt': prt_start_x.ravel(),
'yrpt': prt_start_y.ravel(),
# start particles near top of saturation in cell with local_z=True
'zrpt': particle_release_zrpt, #[1] * len(prt_i) ,
# generic boundname for now;
# could be used to differentiate source areas
'boundname': [f'prt_{cell2d}' for cell2d in particle_start_nodes]
})
#prt_start_data = prt_particle_data
prt_start_data
[17]:
| irptno | k | cell2d | xrpt | yrpt | zrpt | boundname | |
|---|---|---|---|---|---|---|---|
| 0 | 0 | 0 | 1352 | 375.0 | 7125.0 | 0.99 | prt_1352 |
| 1 | 1 | 0 | 1353 | 625.0 | 7125.0 | 0.99 | prt_1353 |
| 2 | 2 | 0 | 1459 | 375.0 | 6875.0 | 0.99 | prt_1459 |
| 3 | 3 | 0 | 1460 | 625.0 | 6875.0 | 0.99 | prt_1460 |
| 4 | 4 | 0 | 1581 | 375.0 | 6625.0 | 0.99 | prt_1581 |
| 5 | 5 | 0 | 1582 | 625.0 | 6625.0 | 0.99 | prt_1582 |
Alternatively, the Modflow 6 Examples illustrate functionality in Flopy to generate evenly spaced particles along cell faces for MODPATH models, and then convert the results to PRT format (if you can get it to work; search the MODFLOW 6 Examples repository for “ModflowPrtprp”)
[18]:
sd = flopy.modpath.CellDataType()
mp7_particle_data = flopy.modpath.NodeParticleData(subdivisiondata=[sd],
nodes=list(particle_start_nodes))
prt_particle_data = list(mp7_particle_data.to_prp(prt.modelgrid))
prt_particle_data
[18]:
[(0,
(0, 1352),
np.float64(291.6666666666667),
np.float64(7041.666666666667),
np.float64(25.157748113333334)),
(1,
(0, 1352),
np.float64(375.0),
np.float64(7041.666666666667),
np.float64(25.157748113333334)),
(2,
(0, 1352),
np.float64(458.33333333333337),
np.float64(7041.666666666667),
np.float64(25.157748113333334)),
(3,
(0, 1352),
np.float64(291.6666666666667),
np.float64(7125.0),
np.float64(25.157748113333334)),
(4,
(0, 1352),
np.float64(375.0),
np.float64(7125.0),
np.float64(25.157748113333334)),
(5,
(0, 1352),
np.float64(458.33333333333337),
np.float64(7125.0),
np.float64(25.157748113333334)),
(6,
(0, 1352),
np.float64(291.6666666666667),
np.float64(7208.333333333334),
np.float64(25.157748113333334)),
(7,
(0, 1352),
np.float64(375.0),
np.float64(7208.333333333334),
np.float64(25.157748113333334)),
(8,
(0, 1352),
np.float64(458.33333333333337),
np.float64(7208.333333333334),
np.float64(25.157748113333334)),
(9,
(0, 1352),
np.float64(291.6666666666667),
np.float64(7041.666666666667),
np.float64(26.95494016)),
(10,
(0, 1352),
np.float64(375.0),
np.float64(7041.666666666667),
np.float64(26.95494016)),
(11,
(0, 1352),
np.float64(458.33333333333337),
np.float64(7041.666666666667),
np.float64(26.95494016)),
(12,
(0, 1352),
np.float64(291.6666666666667),
np.float64(7125.0),
np.float64(26.95494016)),
(13,
(0, 1352),
np.float64(375.0),
np.float64(7125.0),
np.float64(26.95494016)),
(14,
(0, 1352),
np.float64(458.33333333333337),
np.float64(7125.0),
np.float64(26.95494016)),
(15,
(0, 1352),
np.float64(291.6666666666667),
np.float64(7208.333333333334),
np.float64(26.95494016)),
(16,
(0, 1352),
np.float64(375.0),
np.float64(7208.333333333334),
np.float64(26.95494016)),
(17,
(0, 1352),
np.float64(458.33333333333337),
np.float64(7208.333333333334),
np.float64(26.95494016)),
(18,
(0, 1352),
np.float64(291.6666666666667),
np.float64(7041.666666666667),
np.float64(28.752132206666666)),
(19,
(0, 1352),
np.float64(375.0),
np.float64(7041.666666666667),
np.float64(28.752132206666666)),
(20,
(0, 1352),
np.float64(458.33333333333337),
np.float64(7041.666666666667),
np.float64(28.752132206666666)),
(21,
(0, 1352),
np.float64(291.6666666666667),
np.float64(7125.0),
np.float64(28.752132206666666)),
(22,
(0, 1352),
np.float64(375.0),
np.float64(7125.0),
np.float64(28.752132206666666)),
(23,
(0, 1352),
np.float64(458.33333333333337),
np.float64(7125.0),
np.float64(28.752132206666666)),
(24,
(0, 1352),
np.float64(291.6666666666667),
np.float64(7208.333333333334),
np.float64(28.752132206666666)),
(25,
(0, 1352),
np.float64(375.0),
np.float64(7208.333333333334),
np.float64(28.752132206666666)),
(26,
(0, 1352),
np.float64(458.33333333333337),
np.float64(7208.333333333334),
np.float64(28.752132206666666))]
Instantiate the MODFLOW 6 prt particle release point (PRP) package
Notes:
local_z—indicates that “zrpt” [in the particle start data] defines the local z coordinate of the release point within the cell, with value of 0 at the bottom and 1 at the top of the cell. If the cell is partially saturated at release time, the top of the cell is considered to be the water table elevation (the head in the cell) rather than the top defined by the user.
istopzone—integer value defining the stop zone number. If cells have been assigned IZONE values in the GRIDDATA block, a particle terminates if it enters a cell whose IZONE value matches ISTOPZONE. An ISTOPZONE value of zero indicates that there is no stop zone. The default value is zero. This way we can allow weak sinks but have exceptions where the particles will stop anyways.
See the MODFLOW 6 Description of Input and Output for more options
[19]:
flopy.mf6.ModflowPrtprp(
prt,
nreleasepts=len(prt_start_data),
packagedata=prt_start_data.to_records(index=False).tolist(),
local_z=True,
boundnames=True,
perioddata={0: ["FIRST"]},
exit_solve_tolerance=1e-5,
extend_tracking=True,
#istopzone=istopzone
)
[19]:
package_name = prp_0
filename = project-prt.prp
package_type = prp
model_or_simulation_package = model
model_name = project-prt
Block options
--------------------
boundnames
{internal}
(True)
exit_solve_tolerance
{internal}
(1e-05)
local_z
{internal}
(True)
extend_tracking
{internal}
(True)
Block dimensions
--------------------
nreleasepts
{internal}
(6)
Block packagedata
--------------------
packagedata
{internal}
(rec.array([(0, (0, 1352), 375., 7125., 0.99, 'prt_1352'),
(1, (0, 1353), 625., 7125., 0.99, 'prt_1353'),
(2, (0, 1459), 375., 6875., 0.99, 'prt_1459'),
(3, (0, 1460), 625., 6875., 0.99, 'prt_1460'),
(4, (0, 1581), 375., 6625., 0.99, 'prt_1581'),
(5, (0, 1582), 625., 6625., 0.99, 'prt_1582')],
dtype=[('irptno', '<i8'), ('cellid', 'O'), ('xrpt', '<f8'), ('yrpt', '<f8'), ('zrpt', '<f8'), ('boundname', 'O')]))
Block period
--------------------
perioddata
{internal}
(rec.array([('FIRST', None)],
dtype=[('releasesetting', 'O'), ('releasesetting_data', 'O')]))
Instantiate the MODFLOW 6 PRT output control package
Notes: PRT outputs a record each time a: 0: particle was released 1: particle exited a cell 2: time step ended 3: particle terminated 4: particle entered a weak sink cell 5: user-specified tracking time
The code below shows how to input user-specified tracking times to Flopy. Depending on the problem, this may not be necessary in a real-world context, as there may already be many records from particles exiting cells. Therefore, tracktimes=None (no user-specified times) is ultimately input to Flopy.
[20]:
budget_record = [budgetfile_prt]
track_record = [trackfile_prt]
trackcsv_record = [trackcsvfile_prt]
# track positions every year for 100 years
track_nyears = 100
tracktimes = np.linspace(0, track_nyears*365.25, track_nyears+1)
flopy.mf6.ModflowPrtoc(
prt,
budget_filerecord=budget_record,
track_filerecord=track_record,
trackcsv_filerecord=trackcsv_record,
ntracktimes=0,#len(tracktimes),
tracktimes=None,#[(t,) for t in tracktimes],
saverecord=[("BUDGET", "ALL")],
)
[20]:
package_name = oc
filename = project-prt.oc
package_type = oc
model_or_simulation_package = model
model_name = project-prt
Block options
--------------------
budget_filerecord
{internal}
(rec.array([('project-prt.cbc',)],
dtype=[('budgetfile', 'O')]))
track_filerecord
{internal}
(rec.array([('project-prt.trk',)],
dtype=[('trackfile', 'O')]))
trackcsv_filerecord
{internal}
(rec.array([('project-prt.trk.csv',)],
dtype=[('trackcsvfile', 'O')]))
Block dimensions
--------------------
ntracktimes
{internal}
(0)
Block period
--------------------
saverecord
{internal}
(rec.array([('BUDGET', 'ALL', None)],
dtype=[('rtype', 'O'), ('ocsetting', 'O'), ('ocsetting_data', 'O')]))
printrecord
None
Instantiate the PRT Flow Model Interface and Explicit Model Solution Packages
[21]:
fmi_pd = [
("GWFHEAD", os.path.relpath(gwf_ws / headfile, prt_model_ws)), #Path(f"../{gwf_ws.name}/{headfile}")),
("GWFBUDGET", os.path.relpath(gwf_ws / budgetfile, prt_model_ws)) #Path(f"../{gwf_ws.name}/{budgetfile}")),
]
flopy.mf6.ModflowPrtfmi(prt, packagedata=fmi_pd)
# Create an explicit model solution (EMS) for the MODFLOW 6 prt model
ems = flopy.mf6.ModflowEms(
prt_sim,
filename=f"{prt_name}.ems",
)
sim.register_solution_package(ems, [prt.name])
Write the PRT input files and run the model
[22]:
prt_sim.write_simulation()
writing simulation...
writing simulation name file...
writing simulation tdis package...
writing solution package ems_-1...
writing model project-prt...
writing model name file...
writing package disv...
writing package mip...
writing package prp_0...
writing package oc...
writing package fmi...
[23]:
prt_sim.run_simulation()
FloPy is using the following executable to run the model: ../../../../../../../../../../.local/bin/mf6
MODFLOW 6
U.S. GEOLOGICAL SURVEY MODULAR HYDROLOGIC MODEL
VERSION 6.6.2 05/12/2025
MODFLOW 6 compiled May 24 2025 11:41:46 with Intel(R) Fortran Intel(R) 64
Compiler Classic for applications running on Intel(R) 64, Version 2021.7.0
Build 20220726_000000
This software has been approved for release by the U.S. Geological
Survey (USGS). Although the software has been subjected to rigorous
review, the USGS reserves the right to update the software as needed
pursuant to further analysis and review. No warranty, expressed or
implied, is made by the USGS or the U.S. Government as to the
functionality of the software and related material nor shall the
fact of release constitute any such warranty. Furthermore, the
software is released on condition that neither the USGS nor the U.S.
Government shall be held liable for any damages resulting from its
authorized or unauthorized use. Also refer to the USGS Water
Resources Software User Rights Notice for complete use, copyright,
and distribution information.
MODFLOW runs in SEQUENTIAL mode
Run start date and time (yyyy/mm/dd hh:mm:ss): 2025/09/26 14:26:12
Writing simulation list file: mfsim.lst
Using Simulation name file: mfsim.nam
Solving: Stress period: 1 Time step: 1
Run end date and time (yyyy/mm/dd hh:mm:ss): 2025/09/26 14:26:12
Elapsed run time: 0.101 Seconds
Normal termination of simulation.
[23]:
(True, [])
Post-Process the MODFLOW and PRT Results
Load MODFLOW and PRT results from the heads and pathline files
Load the MODFLOW heads
[24]:
hobj = gwf.output.head()
[25]:
hds = hobj.get_data()
Load the tracking CSV file
[26]:
particle_data = pd.read_csv(prt_model_ws / trackcsvfile_prt)
particle_data.head()
[26]:
| kper | kstp | imdl | iprp | irpt | ilay | icell | izone | istatus | ireason | trelease | t | x | y | z | name | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1353 | 0 | 1 | 0 | 0.0 | 0.000000 | 375.000000 | 7125.0 | 29.843871 | PRT_1352 |
| 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1353 | 0 | 1 | 1 | 0.0 | 0.000000 | 375.000000 | 7125.0 | 29.650728 | PRT_1352 |
| 2 | 1 | 1 | 1 | 1 | 1 | 1 | 1353 | 0 | 1 | 1 | 0.0 | 5600.493437 | 378.953010 | 7250.0 | 29.650728 | PRT_1352 |
| 3 | 1 | 1 | 1 | 1 | 1 | 1 | 1198 | 0 | 1 | 1 | 0.0 | 11433.142331 | 402.113867 | 7500.0 | 29.417818 | PRT_1352 |
| 4 | 1 | 1 | 1 | 1 | 1 | 1 | 1043 | 0 | 1 | 1 | 0.0 | 15357.870522 | 425.233970 | 7750.0 | 29.079618 | PRT_1352 |
Plot the heads and pathlines
[27]:
fig, ax = plt.subplots(figsize=(5, 9))
mm = flopy.plot.PlotMapView(model=gwf, layer=0, ax=ax)
mm.plot_array(hds, masked_values=[1e30])
mm.plot_bc('SFR', color='b', plotAll=True)
mm.plot_bc('WEL', plotAll=True)
mm.plot_ibound()
mm.plot_pathline(particle_data, layer='all', color='blue', lw=1)
mm.plot_grid(lw=0.2, color="0.5")
ax = plt.gca()
ax.fill(csx, csy, color='#e534eb', zorder=100, alpha=.75);
[28]:
fig, ax = plt.subplots(figsize=(5, 9))
mm = flopy.plot.PlotMapView(model=gwf, ax=ax, layer=0)
mm.plot_array(hds, masked_values=[1e30])
mm.plot_bc('SFR', color='b', plotAll=True)
mm.plot_bc('WEL', plotAll=True)
mm.plot_ibound()
mm.plot_grid(lw=0.2, color="0.5")
df = particle_data
vmin, vmax = df['t'].min(), df['t'].max()
times = list(range(0, 74001, 1000))
for ix in range(1, len(times)):
tmp = df[(df['t'] >= times[ix - 1]) & (df['t'] < times[ix])]
s = ax.scatter(tmp['x'].values, tmp['y'].values, s=5, c=tmp['t'].values,
vmin=vmin, vmax=vmax, cmap="magma")
ax.set_title(f"{times[ix - 1]} - {times[ix]} days")
if ix == 1:
cbar = fig.colorbar(s, shrink=0.7)
cbar.set_label('particle travel time from source', rotation=270, labelpad=14)
display(fig)
clear_output(wait=True)
plt.pause(0.1)
Export the pathlines to a GeoPackage
GeoPackages allow for subsequent adding of layers; delete any pre-existing version of this GeoPackage to clear any pre-existing layers (often good practice in real projects).
[29]:
output_geopackage = Path(prt_model_ws / 'Pathlines.gpkg')
output_geopackage.unlink(missing_ok=True)
In the real world, the model coordinates also need to be converted to real world coordinates.
We may also want to convert the travel times from the typical model time units of days to years.
We may want to decompose the
cellidinto model layer (k) andcell2d(location in each layer, in this case), or layer, row column if the model has a regular structured grid.Finally, the GeoDataFrame will need a coordinate reference (to write to the GeoPackage) defining the coordinate reference system for the real-world coordinates (for example, UTM zone 15 north, represented by EPSG code 29615).
[30]:
gwf.modelgrid.crs # no CRS attached
[31]:
gwf.modelgrid.crs = 26915
[32]:
df['t_years'] = df['t'] * 325.25
df['k'] = df['icell'].values[-1] // ncells_per_layer
df['cell2d'] = df['icell'].values[-1] % ncells_per_layer
[33]:
x_crs, y_crs = gwf.modelgrid.get_coords(df['x'], df['y'])
gdf = gpd.GeoDataFrame(df,
geometry=gpd.points_from_xy(x_crs, y_crs),
crs=gwf.modelgrid.crs)
Write particles terminating in Well and SFR cells to separate layers
[34]:
for package, izone in izones.items():
layer_name = f"Particles going to {package.upper()} cells"
izone_df = gdf.loc[gdf['izone'] == izone]
izone_df.to_file(output_geopackage, index=False, layer=layer_name)
gdf.loc[gdf['izone'] == 0].to_file(
output_geopackage, index=False, layer='All other particles')
/home/runner/micromamba/envs/pyclass-docs/lib/python3.11/site-packages/pyogrio/__init__.py:7: DeprecationWarning: The 'shapely.geos' module is deprecated, and will be removed in a future version. All attributes of 'shapely.geos' are available directly from the top-level 'shapely' namespace (since shapely 2.0.0).
import shapely.geos # noqa: F401
The previous layers that we added to the GeoPackage were points representing particle locations at discrete points in time. We can also combine these in to pathlines (linestrings) to better visualize the paths taken by the particles. In real-world projects, for example where are particle is started in every cell of a large model, we may also have too many points to work with easily; showing a single pathline for each particle may be advantageous.
[35]:
from shapely.geometry import LineString
# create set of pathlines with start and end information
# first group by particle
by_particle = gdf.sort_values(by=['irpt', 't']).groupby('irpt')
lines = by_particle.last()
lines.crs = gdf.crs # CRS isn't retained in groupby
lines.rename(columns={'k': 'end_k', 'cell2d': 'end_cell2d'}, inplace=True)
line_starts = by_particle.first()
lines['start_k'] = line_starts['k']
lines['start_cell2d'] = line_starts['cell2d']
# above we sorted the results by particle and then time
# use the .agg method on the grouby object to create a LineString
# from the (sorted) sequence of points defining each particle path
linestring_geoms = by_particle['geometry'].agg(lambda x: LineString(x))
lines['geometry'] = linestring_geoms
[36]:
for package, izone in izones.items():
layer_name = f"Pathlines going to {package.upper()} cells"
izone_df = lines.loc[lines['izone'] == izone]
izone_df.to_file(output_geopackage, index=False, layer=layer_name)
Plot pathlines going to SFR and WEL cells
[37]:
fig, ax = plt.subplots(figsize=(6, 10))
pmv = flopy.plot.PlotMapView(gwf, ax=ax)
pmv.plot_grid(lw=0.5)
line_colors = {
'sfr': 'b',
'wel': 'r'
}
for package, izone in izones.items():
izone_df = lines.loc[lines['izone'] == izone]
izone_df.plot(ax=ax, ec=line_colors[package],
label=f'Pathlines going to {package.upper()} cells')
pmv.plot_ibound()
pmv.ax.legend()
[37]:
<matplotlib.legend.Legend at 0x7f6dfb37a7d0>
