10: Particle tracking with MODPATH¶
In this exercise, we will use MODPATH to simulate advective transport with the Freyberg flow model.
MODPATH is a particle tracking software that simulates the advective transport of particles using the output from a MODFLOW groundwater flow model. The MODPATH user guide, documentation, and executable can be found here. The MODPATH executable is downloaded with the stack of modflow related programs when get_modflow() was run in setup for this course. Students should already have this program on their machines.
For this exercise, we will use a quadtree version of the Freyberg model.
[1]:
from IPython.display import clear_output, display
import pathlib as pl
import numpy as np
import matplotlib.pyplot as plt
import geopandas as gpd
from flopy.utils.gridintersect import GridIntersect
import flopy
import pandas as pd
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 and model names.¶
[3]:
load_ws = pl.Path('data/quadtree')
ws = pl.Path("temp/ex10a")
name = "project"
name_mp = f"{name}_mp"
exe_name = 'mf6'
Load the MODFLOW 6 Model¶
Load a simulation object using flopy.mf6.MFSimulation().load().
[4]:
%%capture
sim = flopy.mf6.MFSimulation.load(sim_name=name, exe_name=exe_name,
sim_ws=load_ws)
Load the groundwater flow model¶
[5]:
gwf = sim.get_model(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(ws)
Write the model files¶
[7]:
%%capture
sim.write_simulation()
Run the simulation.
[8]:
sim.run_simulation()
FloPy is using the following executable to run the model: ../../../../../../../.local/share/flopy/bin/mf6
MODFLOW 6
U.S. GEOLOGICAL SURVEY MODULAR HYDROLOGIC MODEL
VERSION 6.4.2 06/28/2023
MODFLOW 6 compiled Jul 05 2023 20:29:14 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.
Run start date and time (yyyy/mm/dd hh:mm:ss): 2024/03/20 14:00:05
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): 2024/03/20 14:00:06
Elapsed run time: 0.603 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.
[8]:
(True, [])
Create and Run the MODPATH model¶
Lets plot the model grid and the location of the contamination patch.¶
[9]:
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');
Find the node numbers of the contamination patch using the groundwater model grid object and the GridIntersect utility. The GridIntersect utility has an intersect method that can be used to identified cellids from points, lines, and polygons.
[10]:
gx = GridIntersect(gwf.modelgrid)
results = gx.intersect(polygon, 'Polygon')
nodes = results.cellids
nodes
[10]:
array([1352, 1353, 1459, 1460, 1581, 1582], dtype=object)
Create the MODPATH particle data¶
[11]:
# create
sd = flopy.modpath.CellDataType()
p = flopy.modpath.NodeParticleData(subdivisiondata=[sd],
nodes=list(nodes))
# create forward particle group
fpth = name_mp + '.sloc'
pg = flopy.modpath.ParticleGroupNodeTemplate(particlegroupname='contaminant',
particledata=p,
filename=fpth)
Create the MODPATH 7 files¶
[12]:
# create modpath files
mp = flopy.modpath.Modpath7(modelname=name_mp, flowmodel=gwf,
exe_name='mp7', model_ws=ws)
mpbas = flopy.modpath.Modpath7Bas(mp, porosity=0.1)
mpsim = flopy.modpath.Modpath7Sim(mp, simulationtype='pathline',
trackingdirection='forward',
weaksinkoption='pass_through',
weaksourceoption='pass_through',
referencetime=0.,
stoptimeoption='extend',
particlegroups=pg)
help(flopy.modpath.Modpath7Sim)
Help on class Modpath7Sim in module flopy.modpath.mp7sim:
class Modpath7Sim(flopy.pakbase.Package)
| Modpath7Sim(model, mpnamefilename=None, listingfilename=None, endpointfilename=None, pathlinefilename=None, timeseriesfilename=None, tracefilename=None, simulationtype='pathline', trackingdirection='forward', weaksinkoption='stop_at', weaksourceoption='stop_at', budgetoutputoption='no', traceparticledata=None, budgetcellnumbers=None, referencetime=None, stoptimeoption='extend', stoptime=None, timepointdata=None, zonedataoption='off', stopzone=None, zones=0, retardationfactoroption='off', retardation=1.0, particlegroups=None, extension='mpsim')
|
| MODPATH Simulation File Package Class.
|
| Parameters
| ----------
| model : model object
| The model object (of type :class:`flopy.modpath.Modpath7`) to
| which this package will be added.
| mpnamefilename : str
| Filename of the MODPATH 7 name file. If mpnamefilename is not
| defined it will be generated from the model name
| (default is None).
| listingfilename : str
| Filename of the MODPATH 7 listing file. If listingfilename is not
| defined it will be generated from the model name
| (default is None).
| endpointfilename : str
| Filename of the MODPATH 7 endpoint file. If endpointfilename is
| not defined it will be generated from the model name
| (default is None).
| pathlinefilename : str
| Filename of the MODPATH 7 pathline file. If pathlinefilename is
| not defined it will be generated from the model name
| (default is None).
| timeseriesfilename : str
| Filename of the MODPATH 7 timeseries file. If timeseriesfilename
| is not defined it will be generated from the model name
| (default is None).
| tracefilename : str
| Filename of the MODPATH 7 tracefile file. If tracefilename is not
| defined it will be generated from the model name
| (default is None).
| simulationtype : str
| MODPATH 7 simulation type. Valid simulation types are 'endpoint',
| 'pathline', 'timeseries', or 'combined' (default is 'pathline').
| trackingdirection : str
| MODPATH 7 tracking direction. Valid tracking directions are
| 'forward' or 'backward' (default os 'forward').
| weaksinkoption : str
| MODPATH 7 weak sink option. Valid weak sink options are
| 'pass_through' or 'stop_at' (default value is 'stop_at').
| weaksourceoption : str
| MODPATH 7 weak source option. Valid weak source options are
| 'pass_through' or 'stop_at' (default value is 'stop_at').
| budgetoutputoption : str
| MODPATH 7 budget output option. Valid budget output options are
| 'no' - individual cell water balance errors are not computed
| and budget record headers are not printed, 'summary' - a summary
| of individual cell water balance errors for each time step is
| printed in the listing file without record headers, or
| 'record_summary' - a summary of individual cell water balance
| errors for each time step is printed in the listing file with
| record headers (default is 'summary').
| traceparticledata : list or tuple
| List or tuple with two ints that define the particle group and
| particle id (zero-based) of the specified particle that is
| followed in detail. If traceparticledata is None, trace mode is
| off (default is None).
| budgetcellnumbers : int, list of ints, tuple of ints, or np.ndarray
| Cell numbers (zero-based) for which detailed water budgets are
| computed. If budgetcellnumbers is None, detailed water budgets are
| not calculated (default is None).
| referencetime : float, list, or tuple
| Specified reference time if a float or a list/tuple with a single
| float value is provided (reference time option 1). Otherwise a
| list or tuple with a zero-based stress period (int) and time
| step (int) and a float defining the relative time position in the
| time step is provided (reference time option 2). If referencetime
| is None, reference time is set to 0 (default is None).
| stoptimeoption : str
| String indicating how a particle tracking simulation is
| terminated based on time. If stop time option is 'total', particles
| will be stopped at the end of the final time step if 'forward'
| tracking is simulated or at the beginning of the first time step
| if backward tracking. If stop time option is 'extend', initial or
| final steady-state time steps will be extended and all particles
| will be tracked until they reach a termination location. If stop
| time option is 'specified', particles will be tracked until they
| reach a termination location or the specified stop time is reached
| (default is 'extend').
| stoptime : float
| User-specified value of tracking time at which to stop a particle
| tracking simulation. Stop time is only used if the stop time option
| is 'specified'. If stoptime is None and the stop time option is
| 'specified' particles will be terminated at the end of the last
| time step if 'forward' tracking or the beginning of the first time
| step if 'backward' tracking (default is None).
| timepointdata : list or tuple
| List or tuple with 2 items that is only used if simulationtype is
| 'timeseries' or 'combined'. If the second item is a float then the
| timepoint data corresponds to time point option 1 and the first
| entry is the number of time points (timepointcount) and the second
| entry is the time point interval. If the second item is a list,
| tuple, or np.ndarray then the timepoint data corresponds to time
| point option 2 and the number of time points entries
| (timepointcount) in the second item and the second item is an
| list, tuple, or array of user-defined time points. If Timepointdata
| is None, time point option 1 is specified and the total simulation
| time is split into 100 intervals (default is None).
| zonedataoption : str
| If zonedataoption is 'off', zone array data are not read and a zone
| value of 1 is applied to all cells. If zonedataoption is 'on',
| zone array data are read (default is 'off').
| stopzone : int
| A zero-based specified integer zone value that indicates an
| automatic stopping location for particles and is only used if
| zonedataoption is 'on'. A value of -1 indicates no automatic stop
| zone is used. Stopzone values less than -1 are not allowed. If
| stopzone is None, stopzone is set to -1 (default is None).
| zones : float or array of floats (nlay, nrow, ncol)
| Array of zero-based positive integer zones that are only used if
| zonedataoption is 'on' (default is 0).
| retardationfactoroption : str
| If retardationfactoroption is 'off', retardation array data are not
| read and a retardation factor of 1 is applied to all cells. If
| retardationfactoroption is 'on', retardation factor array data are
| read (default is 'off').
| retardation : float or array of floats (nlay, nrow, ncol)
| Array of retardation factors that are only used if
| retardationfactoroption is 'on' (default is 1).
| particlegroups : ParticleGroup or list of ParticleGroups
| ParticleGroup or list of ParticlesGroups that contain data for
| individual particle groups. If None is specified, a
| particle in the center of node 0 will be created (default is None).
| extension : string
| Filename extension (default is 'mpsim')
|
| Examples
| --------
|
| >>> import flopy
| >>> m = flopy.modflow.Modflow.load('mf2005.nam')
| >>> mp = flopy.modpath.Modpath7('mf2005_mp', flowmodel=m)
| >>> mpsim = flopy.modpath.Modpath7Sim(mp)
|
| Method resolution order:
| Modpath7Sim
| flopy.pakbase.Package
| flopy.pakbase.PackageInterface
| builtins.object
|
| Methods defined here:
|
| __init__(self, model, mpnamefilename=None, listingfilename=None, endpointfilename=None, pathlinefilename=None, timeseriesfilename=None, tracefilename=None, simulationtype='pathline', trackingdirection='forward', weaksinkoption='stop_at', weaksourceoption='stop_at', budgetoutputoption='no', traceparticledata=None, budgetcellnumbers=None, referencetime=None, stoptimeoption='extend', stoptime=None, timepointdata=None, zonedataoption='off', stopzone=None, zones=0, retardationfactoroption='off', retardation=1.0, particlegroups=None, extension='mpsim')
| Initialize self. See help(type(self)) for accurate signature.
|
| write_file(self, check=False)
| Write the package file
|
| Parameters
| ----------
| check : boolean
| Check package data for common errors. (default False)
|
| Returns
| -------
| None
|
| ----------------------------------------------------------------------
| Methods inherited from flopy.pakbase.Package:
|
| __getitem__(self, item)
|
| __repr__(self)
| Return repr(self).
|
| __setattr__(self, key, value)
| Implement setattr(self, name, value).
|
| __setitem__(self, key, value)
|
| export(self, f, **kwargs)
| Method to export a package to netcdf or shapefile based on the
| extension of the file name (.shp for shapefile, .nc for netcdf)
|
| Parameters
| ----------
| f : str
| filename
| kwargs : keyword arguments
| modelgrid : flopy.discretization.Grid instance
| user supplied modelgrid which can be used for exporting
| in lieu of the modelgrid associated with the model object
|
| Returns
| -------
| None or Netcdf object
|
| level1_arraylist(self, idx, v, name, txt)
|
| plot(self, **kwargs)
| Plot 2-D, 3-D, transient 2-D, and stress period list (MfList)
| package input data
|
| Parameters
| ----------
| **kwargs : dict
| filename_base : str
| Base file name that will be used to automatically generate file
| names for output image files. Plots will be exported as image
| files if file_name_base is not None. (default is None)
| file_extension : str
| Valid matplotlib.pyplot file extension for savefig(). Only used
| if filename_base is not None. (default is 'png')
| mflay : int
| MODFLOW zero-based layer number to return. If None, then all
| all layers will be included. (default is None)
| kper : int
| MODFLOW zero-based stress period number to return. (default is
| zero)
| key : str
| MfList dictionary key. (default is None)
|
| Returns
| ----------
| axes : list
| Empty list is returned if filename_base is not None. Otherwise
| a list of matplotlib.pyplot.axis are returned.
|
| See Also
| --------
|
| Notes
| -----
|
| Examples
| --------
| >>> import flopy
| >>> ml = flopy.modflow.Modflow.load('test.nam')
| >>> ml.dis.plot()
|
| set_cbc_output_file(self, ipakcb, model, fname)
|
| to_shapefile(self, filename, **kwargs)
| Export 2-D, 3-D, and transient 2-D model data to shapefile (polygons).
| Adds an attribute for each layer in each data array
|
| Parameters
| ----------
| filename : str
| Shapefile name to write
|
| Returns
| ----------
| None
|
| See Also
| --------
|
| Notes
| -----
|
| Examples
| --------
| >>> import flopy
| >>> ml = flopy.modflow.Modflow.load('test.nam')
| >>> ml.lpf.to_shapefile('test_hk.shp')
|
| webdoc(self)
| Open the web documentation.
|
| ----------------------------------------------------------------------
| Static methods inherited from flopy.pakbase.Package:
|
| add_to_dtype(dtype, field_names, field_types)
| Add one or more fields to a structured array data type
|
| Parameters
| ----------
| dtype : numpy.dtype
| Input structured array datatype to add to.
| field_names : str or list
| One or more field names.
| field_types : numpy.dtype or list
| One or more data types. If one data type is supplied, it is
| repeated for each field name.
|
| load(f: Union[str, bytes, os.PathLike], model, pak_type, ext_unit_dict=None, **kwargs)
| Default load method for standard boundary packages.
|
| ----------------------------------------------------------------------
| Readonly properties inherited from flopy.pakbase.Package:
|
| data_list
|
| package_type
|
| plottable
|
| ----------------------------------------------------------------------
| Data descriptors inherited from flopy.pakbase.Package:
|
| name
|
| parent
|
| ----------------------------------------------------------------------
| Methods inherited from flopy.pakbase.PackageInterface:
|
| check(self, f=None, verbose=True, level=1, checktype=None)
| Check package data for common errors.
|
| Parameters
| ----------
| f : str or file handle
| String defining file name or file handle for summary file
| of check method output. If a sting is passed a file handle
| is created. If f is None, check method does not write
| results to a summary file. (default is None)
| verbose : bool
| Boolean flag used to determine if check method results are
| written to the screen
| level : int
| Check method analysis level. If level=0, summary checks are
| performed. If level=1, full checks are performed.
| checktype : check
| Checker type to be used. By default class check is used from
| check.py.
|
| Returns
| -------
| None
|
| Examples
| --------
|
| >>> import flopy
| >>> m = flopy.modflow.Modflow.load('model.nam')
| >>> m.dis.check()
|
| ----------------------------------------------------------------------
| Readonly properties inherited from flopy.pakbase.PackageInterface:
|
| has_stress_period_data
|
| ----------------------------------------------------------------------
| Data descriptors inherited from flopy.pakbase.PackageInterface:
|
| __dict__
| dictionary for instance variables
|
| __weakref__
| list of weak references to the object
Write MODPATH 7 files and run the model¶
[13]:
# write modpath datasets
mp.write_input()
# run modpath
mp.run_model()
FloPy is using the following executable to run the model: ../../../../../../../.local/share/flopy/bin/mp7
MODPATH Version 7.2.001
Program compiled Jul 05 2023 20:35:22 with IFORT compiler (ver. 20.21.7)
Run particle tracking simulation ...
Processing Time Step 1 Period 1. Time = 1.00000E+00 Steady-state flow
Particle Summary:
0 particles are pending release.
0 particles remain active.
0 particles terminated at boundary faces.
0 particles terminated at weak sink cells.
0 particles terminated at weak source cells.
162 particles terminated at strong source/sink cells.
0 particles terminated in cells with a specified zone number.
0 particles were stranded in inactive or dry cells.
0 particles were unreleased.
0 particles have an unknown status.
Normal termination.
[13]:
(True, [])
Post-Process the MODFLOW and MODPATH Results¶
Load MODFLOW and MODPATH results from the heads and pathline files¶
Load the MODFLOW heads
[14]:
hobj = gwf.output.head()
[15]:
hds = hobj.get_data()
Load the pathline file
[16]:
ppth =ws / f"{name_mp}.mppth"
p = flopy.utils.PathlineFile(ppth)
p0 = p.get_alldata()
[17]:
p0[0:2]
[17]:
[array([( 291.66666, 7041.6665, 25.157747 , 0. , 0, 0),
( 292.26523, 7053.434 , 24.259151 , 1066.6232, 0, 0),
( 296.87158, 7187.315 , 18.867577 , 9176.828 , 1, 0),
( 298.29572, 7250. , 17.83786 , 11656.656 , 2, 0),
( 327.38013, 7500. , 17.041315 , 17482.574 , 2, 0),
( 351.95312, 7750. , 14.886883 , 21403.549 , 2, 0),
( 384.25983, 8000. , 14.738331 , 24241.266 , 2, 0),
( 437.98865, 8250. , 15.478686 , 26416.283 , 2, 0),
( 500. , 8325.292 , 14.165776 , 27109.031 , 2, 0),
( 692.9299 , 8500. , 15.878951 , 28298.334 , 2, 0),
( 750. , 8540.724 , 17.176243 , 28605.803 , 2, 0),
(1000. , 8630.868 , 14.297673 , 30072.422 , 2, 0),
(1250. , 8669.552 , 14.05823 , 31359.709 , 2, 0),
(1500. , 8683.525 , 15.129906 , 32383.266 , 2, 0),
(1750. , 8648.858 , 15.427488 , 33216.24 , 2, 0),
(2000. , 8599.475 , 13.062349 , 34121.06 , 2, 0),
(2250. , 8531.675 , 14.423839 , 34772.07 , 2, 0),
(2322.0684 , 8500. , 10.878183 , 35016.754 , 2, 0),
(2500. , 8421.822 , 9.152523 , 35671.594 , 2, 0),
(2588.7498 , 8375. , 3.7609835, 36132.93 , 2, 0),
(2625. , 8356.87 , 3.756292 , 36310.203 , 2, 0),
(2750. , 8295.934 , 3.7417758, 36877.59 , 2, 0),
(2849.241 , 8250. , 3.6520696, 37296.168 , 2, 0),
(2875. , 8235.725 , 3.6121962, 37420.03 , 2, 0),
(2937.5 , 8205.341 , 3.7405884, 37695.703 , 2, 0),
(2971.6401 , 8187.5 , 3.7372005, 37843.93 , 2, 0),
(3000. , 8175.1055, 3.7348478, 37948.184 , 2, 0),
(3062.5 , 8156.6216, 3.9858716, 38161.88 , 2, 0),
(3063.7656 , 8156.25 , 4.1625876, 38166.094 , 2, 0),
(3093.75 , 8147.8945, 4.1604686, 38263.008 , 2, 0),
(3125. , 8139.558 , 4.1584296, 38357.234 , 2, 0),
(3140.625 , 8136.6265, 4.202556 , 38401.426 , 2, 0),
(3156.25 , 8133.5854, 4.201312 , 38444.7 , 2, 0),
(3171.875 , 8131.6494, 4.052273 , 38486.99 , 2, 0),
(3187.5 , 8127.9507, 3.803092 , 38593.96 , 2, 0),
(3191.0186 , 8125. , 3.708662 , 38683.953 , 2, 0),
(3203.125 , 8118.2197, 3.697419 , 38890.438 , 2, 0),
(3218.75 , 8110.778 , 3.696752 , 39123.6 , 2, 0),
(3250. , 8098.4146, 3.542844 , 39516.027 , 2, 0),
(3259.4897 , 8093.75 , 3.3440852, 39625.527 , 2, 0),
(3281.25 , 8082.3413, 3.34398 , 39905.668 , 2, 0),
(3312.5 , 8070.0215, 3.3438506, 40262.293 , 2, 0),
(3333.566 , 8062.5 , 3.0362823, 40485.49 , 2, 0),
(3375. , 8047.49 , 3.1637154, 40932.04 , 2, 0),
(3437.5 , 8013.0425, 2.9553792, 41570.11 , 2, 0),
(3464.91 , 8000. , 2.9553003, 41827.246 , 2, 0),
(3500. , 7990.7407, 3.2104254, 42010.31 , 2, 0),
(3625. , 7967.9907, 3.2902322, 42530.023 , 2, 0),
(3680.1426 , 7937.5 , 3.2615664, 42698.07 , 2, 0),
(3687.5 , 7934.937 , 3.2615619, 42714.484 , 2, 0),
(3750. , 7915.0894, 3.2615318, 42823.734 , 2, 0),
(3773.753 , 7906.25 , 3.2376952, 42854.773 , 2, 0),
(3781.25 , 7904.3774, 3.23769 , 42862.258 , 2, 0),
(3812.5 , 7896.45 , 3.237673 , 42886.637 , 2, 0),
(3828.125 , 7891.537 , 3.2256062, 42895.133 , 2, 0),
(3830.2957 , 7890.625 , 3.225588 , 42896.133 , 2, 0),
(3843.75 , 7887.5903, 3.2255182, 42900.066 , 2, 0),
(3859.375 , 7885.2886, 3.2016592, 42902.145 , 2, 0),
(3859.375 , 7885.2886, 3.2188113, 42902.145 , 2, 0)],
dtype={'names': ['x', 'y', 'z', 'time', 'k', 'particleid'], 'formats': ['<f4', '<f4', '<f4', '<f4', '<i4', '<i4'], 'offsets': [20, 24, 28, 16, 32, 0], 'itemsize': 60}),
array([( 375. , 7041.6665, 25.157747 , 0. , 0, 1),
( 375.75186, 7053.434 , 24.259151 , 1066.6232, 0, 1),
( 381.52545, 7187.315 , 18.867577 , 9176.828 , 1, 1),
( 383.30856, 7250. , 17.83786 , 11656.656 , 2, 1),
( 406.18857, 7500. , 17.041315 , 17482.574 , 2, 1),
( 429.23462, 7750. , 14.886883 , 21403.549 , 2, 1),
( 466.02756, 8000. , 14.738331 , 24241.266 , 2, 1),
( 500. , 8134.05 , 15.576507 , 25436.912 , 2, 1),
( 518.4427 , 8250. , 15.027732 , 26109.572 , 2, 1),
( 750. , 8465.881 , 15.941434 , 27393.213 , 2, 1),
( 894.6647 , 8500. , 12.070309 , 28152.492 , 2, 1),
(1000. , 8535.509 , 14.394118 , 28748.236 , 2, 1),
(1250. , 8562.502 , 14.136991 , 30035.523 , 2, 1),
(1500. , 8564.003 , 15.185995 , 31059.08 , 2, 1),
(1750. , 8521.553 , 15.469522 , 31892.055 , 2, 1),
(1840.7961 , 8500. , 13.12992 , 32226.053 , 2, 1),
(2000. , 8464.708 , 13.11985 , 32755.252 , 2, 1),
(2250. , 8386.366 , 14.19391 , 33350.297 , 2, 1),
(2500. , 8257.273 , 9.194104 , 34267.223 , 2, 1),
(2514.2393 , 8250. , 3.8298988, 34341.61 , 2, 1),
(2625. , 8181.72 , 3.796223 , 35022.78 , 2, 1),
(2709.697 , 8125. , 3.7807877, 35533.742 , 2, 1),
(2750. , 8103.26 , 3.7751324, 35727.203 , 2, 1),
(2875. , 8047.326 , 3.7963305, 36297.71 , 2, 1),
(2937.5 , 8022.2764, 3.7470431, 36574.035 , 2, 1),
(2998.7302 , 8000. , 3.7404735, 36823.684 , 2, 1),
(3000. , 7999.431 , 3.7184064, 36830.277 , 2, 1),
(3062.5 , 7978.786 , 3.7055447, 37123.875 , 2, 1),
(3093.75 , 7972.582 , 3.6036906, 37257.2 , 2, 1),
(3112.059 , 7968.75 , 3.601786 , 37334.098 , 2, 1),
(3125. , 7966.217 , 3.6005259, 37385.312 , 2, 1),
(3140.625 , 7964.6694, 3.5051665, 37446.508 , 2, 1),
(3156.25 , 7963.018 , 3.5032501, 37508.477 , 2, 1),
(3171.875 , 7961.6426, 3.5620072, 37567.363 , 2, 1),
(3187.5 , 7959.37 , 3.3982596, 37676.285 , 2, 1),
(3203.125 , 7953.85 , 3.4009736, 37923.64 , 2, 1),
(3205.5557 , 7953.125 , 3.4006464, 37960.184 , 2, 1),
(3218.75 , 7949.4136, 3.3990479, 38150.58 , 2, 1),
(3250. , 7942.7827, 3.453589 , 38556.527 , 2, 1),
(3281.25 , 7937.636 , 3.60839 , 38914.742 , 2, 1),
(3282.212 , 7937.5 , 3.6083763, 38924.887 , 2, 1),
(3312.5 , 7933.544 , 3.3178153, 39235.176 , 2, 1),
(3375. , 7927.682 , 3.427484 , 39777.83 , 2, 1),
(3437.5 , 7925.3003, 3.2190576, 40175.04 , 2, 1),
(3500. , 7913.6294, 3.2186823, 40493.594 , 2, 1),
(3625. , 7912.4023, 3.2980738, 41013.31 , 2, 1),
(3687.5 , 7899.6323, 3.2693596, 41184.34 , 2, 1),
(3750. , 7891.814 , 3.2692423, 41293.59 , 2, 1),
(3781.25 , 7888.684 , 3.2453935, 41329.97 , 2, 1),
(3812.5 , 7885.8335, 3.2453277, 41354.344 , 2, 1),
(3828.125 , 7884.7036, 3.2328875, 41362.043 , 2, 1),
(3843.75 , 7883.1997, 3.2325552, 41366.78 , 2, 1),
(3859.375 , 7883.1904, 3.2065158, 41368.86 , 2, 1),
(3859.375 , 7883.1904, 3.2236557, 41368.86 , 2, 1)],
dtype={'names': ['x', 'y', 'z', 'time', 'k', 'particleid'], 'formats': ['<f4', '<f4', '<f4', '<f4', '<i4', '<i4'], 'offsets': [20, 24, 28, 16, 32, 0], 'itemsize': 60})]
Plot the heads and pathlines¶
[18]:
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(p0, 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);
[19]:
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")
p = flopy.utils.PathlineFile(ppth)
p0 = p.get_alldata()
df = pd.DataFrame.from_records(p0[0])
for ix, recarray in enumerate(p0):
dft = pd.DataFrame.from_records(recarray)
df = pd.concat((df, dft), ignore_index=True)
vmin, vmax = df.time.min(), df.time.max()
times = list(range(0, 74001, 1000))
for ix in range(1, len(times)):
tmp = df[(df.time >= times[ix - 1]) & (df.time < times[ix])]
s = ax.scatter(tmp.x.values, tmp.y.values, c=tmp.time.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 shapefile¶
[20]:
spth = pl.Path(ws / 'pathline.shp')
p.write_shapefile(p0, mg=gwf.modelgrid, one_per_particle=False, shpname=spth)
{'names': ['x', 'y', 'z', 'time', 'k', 'particleid'], 'formats': ['<f4', '<f4', '<f4', '<f4', '<i4', '<i4'], 'offsets': [20, 24, 28, 16, 32, 0], 'itemsize': 60}
No CRS information for writing a .prj file.
Supply an valid coordinate system reference to the attached modelgrid object or .export() method.
Load the shapefile into geopandas¶
[21]:
rshp = gpd.read_file(spth)
Show the dataframe¶
[22]:
rshp
[22]:
| x | y | z | time | k | particleid | geometry | |
|---|---|---|---|---|---|---|---|
| 0 | 292.265228 | 7053.434082 | 24.259151 | 1066.623169 | 1 | 1 | LINESTRING (291.667 7041.667, 292.265 7053.434) |
| 1 | 296.871582 | 7187.314941 | 18.867577 | 9176.828125 | 2 | 1 | LINESTRING (292.265 7053.434, 296.872 7187.315) |
| 2 | 298.295715 | 7250.000000 | 17.837860 | 11656.656250 | 3 | 1 | LINESTRING (296.872 7187.315, 298.296 7250.000) |
| 3 | 327.380127 | 7500.000000 | 17.041315 | 17482.574219 | 3 | 1 | LINESTRING (298.296 7250.000, 327.380 7500.000) |
| 4 | 351.953125 | 7750.000000 | 14.886883 | 21403.548828 | 3 | 1 | LINESTRING (327.380 7500.000, 351.953 7750.000) |
| ... | ... | ... | ... | ... | ... | ... | ... |
| 7326 | 3312.500000 | 5125.320312 | 1.554207 | 49924.312500 | 3 | 162 | LINESTRING (3304.377 5125.000, 3312.500 5125.320) |
| 7327 | 3328.125000 | 5128.071289 | 1.580646 | 49943.050781 | 3 | 162 | LINESTRING (3312.500 5125.320, 3328.125 5128.071) |
| 7328 | 3343.750000 | 5128.837891 | 1.579881 | 49954.789062 | 3 | 162 | LINESTRING (3328.125 5128.071, 3343.750 5128.838) |
| 7329 | 3359.375000 | 5131.090332 | 1.578388 | 49959.777344 | 3 | 162 | LINESTRING (3343.750 5128.838, 3359.375 5131.090) |
| 7330 | 3359.375000 | 5131.090332 | 1.587395 | 49959.777344 | 3 | 162 | LINESTRING (3359.375 5131.090, 3359.375 5131.090) |
7331 rows × 7 columns
Plot the geopandas dataframe on the modelgrid¶
[23]:
fig, ax = plt.subplots(figsize=(6, 10))
pmv = flopy.plot.PlotMapView(gwf, ax=ax)
pmv.plot_grid(lw=0.5)
pmv.plot_shapes(rshp.geometry.to_list(), edgecolor="grey") #, edgecolor="b")
pmv.plot_ibound();
Export the data to VTK (Visualization ToolKit) format for a 3d representation¶
FloPy’s Vtk module allows us to create Visualization ToolKit (VTK) files that can be opened and explored with pyVISTAs or external software such as ParaView.
[24]:
from flopy.export.vtk import Vtk
vtk = Vtk(model=gwf, binary=False, vertical_exageration=50, smooth=False)
vtk.add_model(gwf)
vtk.add_pathline_points(p0)
Write VTK file for use in ParaView or other VTK software
[25]:
pth = pl.Path("temp")
pth.mkdir(exist_ok=True)
vtk.write(pth / "freyberg")
Alternatively we can visualize the VTK representation using pyvista. To activate these cells, set plotvtk=True¶
[26]:
plotvtk = True
if plotvtk:
import pyvista as pv
Convert the VTK to pyvista meshes and rotate the meshes to match the orientation of our matplotlib plots
[27]:
if plotvtk:
grid, pathlines = vtk.to_pyvista()
[28]:
if plotvtk:
axes = pv.Axes(show_actor=True, actor_scale=2.0, line_width=5)
grid.rotate_z(160, point=axes.origin, inplace=True)
pathlines.rotate_z(160, point=axes.origin, inplace=True)
Select particle release locations and build a dictionary of particle tracks (pathlines). This will be used below for particle labelling, as well as for animation.
Note: while below we construct pathlines manually from data read from the exported VTK files, pathlines may also be read directly from the MODPATH 7 pathline output file (provided the simulation was run in pathline or combined mode, as this one was).
[29]:
if plotvtk:
tracks = {}
particle_ids = set()
release_locs = list()
for i, t in enumerate(pathlines["time"]):
pid = str(round(float(pathlines["particleid"][i])))
loc = pathlines.points[i]
if pid not in tracks:
tracks[pid] = []
particle_ids.add(pid)
release_locs.append(loc)
# store the particle location in the corresponding track
tracks[pid].append((loc, t))
release_locs = np.array(release_locs)
tracks = {k: np.array(v, dtype=object) for k, v in tracks.items()}
max_track_len = max([len(v) for v in tracks.values()])
View the grid and pathlines with PyVista, with particle tracks/locations colored by time. Also add particle ID labels to a few particles’ release locations.
[30]:
if plotvtk:
pv.set_plot_theme("document")
pv.set_jupyter_backend("static")
# create the plot and add the grid and pathline meshes
p = pv.Plotter()
p.add_mesh(grid, opacity=0.05)
p.add_mesh(pathlines, scalars="time")
# add a particle ID label to each 20th particle's starting point
label_coords = []
start_labels = []
for pid, track in tracks.items():
if int(pid) % 20 == 0:
label_coords.append(track[0][0])
start_labels.append(f"Particle {pid}")
p.add_point_labels(
label_coords,
start_labels,
font_size=10,
point_size=15,
point_color="black",
)
# zoom in and show the plot
p.camera.zoom(2.4)
p.show()
Create an animated GIF of the particles traveling along their pathlines, with particles colored by time.
[31]:
if plotvtk:
# create plotter
p = pv.Plotter(notebook=False, off_screen=True)
# open GIF file
gif_path = pth / f"freyberg_tracks.gif"
p.open_gif(str(gif_path))
# create mesh from release locations
spls = pv.PolyData(release_locs)
spls.point_data["time"] = np.zeros(len(spls.points))
# add the underlying grid mesh and particle data, then zoom in
p.add_mesh(grid, opacity=0.05)
p.add_mesh(spls, clim=[0, 1.23e09])
p.camera.zoom(2.4)
# cycle through time steps and update particle location
for i in range(1, max_track_len):
pts = []
times = []
segments = []
for pid in particle_ids:
track = tracks[pid]
npts = len(track)
# use last locn if particle has already terminated
loc, t = track[i] if i < npts else track[npts - 1]
pts.append(loc)
times.append(t)
if i < npts:
segments.append(track[i - 1][0])
segments.append(loc)
p.update_coordinates(np.vstack(pts), render=False)
p.update_scalars(np.array(times), mesh=spls, render=False)
p.add_lines(np.array(segments), width=1, color="black")
p.write_frame() # write frame to file
# close the plotter and the GIF file
p.close()
/Users/mnfienen/miniforge3/envs/pyclass/lib/python3.11/site-packages/pyvista/plotting/plotter.py:4722: PyVistaDeprecationWarning: This method is deprecated and will be removed in a future version of PyVista. Directly modify the points of a mesh in-place instead.
warnings.warn(
/Users/mnfienen/miniforge3/envs/pyclass/lib/python3.11/site-packages/pyvista/plotting/plotter.py:4644: PyVistaDeprecationWarning: This method is deprecated and will be removed in a future version of PyVista. Directly modify the scalars of a mesh in-place instead.
warnings.warn(
[32]:
if plotvtk:
from IPython.core.display import Image
display(Image(data=open(gif_path, "rb").read(), format="gif"))
<IPython.core.display.Image object>
[ ]: