{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# 10a: Particle tracking with MODFLOW 6 PRT\n", "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. \n", "\n", "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.\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from IPython.display import clear_output, display\n", "import os\n", "from pathlib import Path\n", "\n", "import numpy as np\n", "import pandas as pd\n", "import geopandas as gpd\n", "import matplotlib.pyplot as plt\n", "from flopy.utils.gridintersect import GridIntersect\n", "import flopy\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The location of the contamination patch and the nodes that the define bounding cells of the patch are calculated below." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# patch upper left and lower right\n", "xmin, xmax = 250. * 1, 250. * 3\n", "ymin, ymax = (40 - 14) * 250., (40 - 11) * 250. \n", "\n", "csx, csy = [xmin, xmin, xmax, xmax, xmin], [ymin, ymax, ymax, ymin, ymin]\n", "polygon = [list(zip(csx, csy))]\n", "(xmin, ymax), (xmax, ymin)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "--------------------------\n", "\n", "### Define the workspace, model names and key options." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "load_ws = Path('data/quadtree')\n", "gwf_ws = Path(\"temp/ex10a\")\n", "gwf_name = \"project\"\n", "#name_mp = f\"{name}_mp\"\n", "exe_name = 'mf6'\n", "\n", "# flow model output files to use with PRT\n", "headfile = f\"{gwf_name}.hds\"\n", "budgetfile = f\"{gwf_name}.cbc\"\n", "\n", "prt_name = f\"{gwf_name}-prt\"\n", "prt_model_ws = gwf_ws / 'prt'\n", "\n", "# PRT output files\n", "budgetfile_prt = f\"{prt_name}.cbc\"\n", "trackfile_prt = f\"{prt_name}.trk\"\n", "trackcsvfile_prt = f\"{prt_name}.trk.csv\"\n", "\n", "# if using \"local_z\" option\n", "# 1=start particles at top of cell\n", "# 0=start particles at bottom of cell\n", "# 1 resulted in Error: release point (z=29.899951638778955) is above grid top 29.886768340000000\n", "particle_release_zrpt = .99\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Load the MODFLOW 6 Model\n", "\n", "Load a simulation object using `flopy.mf6.MFSimulation().load()`." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%capture\n", "sim = flopy.mf6.MFSimulation.load(sim_name=gwf_name, exe_name=exe_name,\n", " sim_ws=load_ws)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Load the groundwater flow model" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gwf = sim.get_model(gwf_name)\n", "gwf.modelgrid" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Change the workspace" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sim.set_sim_path(gwf_ws)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Add ``save_saturation`` and ``save_flows`` to the NPF Package\n", "(required for PRT)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gwf.npf.save_saturation = True\n", "gwf.npf.save_flows = True" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Write the model files" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%capture\n", "sim.write_simulation()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Run the simulation." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sim.run_simulation()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Create and Run the PRT model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plot the model grid and the location of the contamination patch" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(5, 9))\n", "mm = flopy.plot.PlotMapView(gwf, layer=0, ax=ax)\n", "\n", "mm.plot_bc('SFR', color=\"b\", plotAll=True)\n", "mm.plot_bc('WEL', plotAll=True)\n", "mm.plot_inactive(alpha=0.75)\n", "\n", "mm.plot_grid(lw=0.25, color='grey')\n", "\n", "ax.fill(csx, csy, color='#e534eb');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Get the model cells intersecting the contamination patch\n", "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." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gx = GridIntersect(gwf.modelgrid)\n", "results = gx.intersect(polygon, 'Polygon')\n", "particle_start_nodes = results.cellids.astype(int)\n", "particle_start_nodes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Instantiate the MODFLOW 6 prt model\n", "and discretization packages" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "prt_sim = flopy.mf6.MFSimulation(sim_name=prt_name, sim_ws=prt_model_ws)\n", "\n", "# Instantiate the MODFLOW 6 temporal discretization package\n", "flopy.mf6.modflow.mftdis.ModflowTdis(\n", " prt_sim,\n", " time_units=\"DAYS\",\n", " nper=1,\n", " perioddata=[(1, 1, 1)], # perlen, nstp, tsmult\n", ")\n", "prt = flopy.mf6.ModflowPrt(\n", " prt_sim, modelname=prt_name, model_nam_file=f\"{prt_name}.nam\"\n", ")\n", "\n", "# Instantiate the MODFLOW 6 prt discretization package\n", "nlay, ncells_per_layer = gwf.dis.botm.array.shape\n", "disv = flopy.mf6.ModflowGwfdisv(\n", " prt,\n", " nlay=nlay,\n", " ncpl=gwf.dis.ncpl.array,\n", " nvert=gwf.dis.nvert.array,\n", " length_units=gwf.dis.length_units.array,\n", " top=gwf.dis.top.array,\n", " botm=gwf.dis.botm.array,\n", " vertices=gwf.dis.vertices.array,\n", " cell2d=gwf.dis.cell2d.array,\n", " idomain=gwf.dis.idomain.array,\n", " xorigin=gwf.dis.xorigin.array,\n", " yorigin=gwf.dis.yorigin.array,\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Instantiate the MODFLOW 6 PRT Model Input Package.\n", "\n", "**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.\n", "\n", "In this example, we'll assign different izone numbers to SFR and Well Package cells." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# start with a default zone of 0\n", "izone_array = np.zeros((nlay, ncells_per_layer), dtype=int)\n", "\n", "# get the locations of SFR cells\n", "sfr_k, sfr_cellid = zip(*gwf.sfr.packagedata.array['cellid'])\n", "sfr_k[:10], sfr_cellid[:10]\n", "\n", "izones = {\n", " 'wel': 1,\n", " 'sfr': 2\n", "}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Assign SFR cells to zone 2**" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "izone_array[sfr_k, sfr_cellid] = izones['sfr']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Assign wells to zone 1**\n", "Using the stress period data. Note that in models with more than one stress period, different wells may be represented in different stress periods." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "well_k = list()\n", "well_cellid = list()\n", "for per, recarray in gwf.wel.stress_period_data.data.items():\n", " well_k_per, well_cellid_per = zip(*recarray['cellid'])\n", " well_k.append(well_k_per)\n", " well_cellid.append(well_cellid_per)\n", "\n", "izone_array[well_k, well_cellid] = izones['wel']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Make the package**" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "flopy.mf6.ModflowPrtmip(prt,\n", " porosity=0.1, \n", " izone=izone_array\n", " ) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Create the particle start data\n", "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. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# release particles at each row, column location\n", "prt_start_x, prt_start_y = gwf.modelgrid.get_local_coords(\n", " gwf.modelgrid.xcellcenters[particle_start_nodes], \n", " gwf.modelgrid.ycellcenters[particle_start_nodes])\n", "\n", "prt_start_data = pd.DataFrame({\n", " 'irptno': np.arange(len(prt_start_x.ravel())),\n", " 'k': 0,\n", " 'cell2d': particle_start_nodes,\n", " 'xrpt': prt_start_x.ravel(),\n", " 'yrpt': prt_start_y.ravel(),\n", " # start particles near top of saturation in cell with local_z=True\n", " 'zrpt': particle_release_zrpt, #[1] * len(prt_i) ,\n", " # generic boundname for now;\n", " # could be used to differentiate source areas\n", " 'boundname': [f'prt_{cell2d}' for cell2d in particle_start_nodes]\n", "})\n", "#prt_start_data = prt_particle_data\n", "prt_start_data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Alternatively, the [Modflow 6 Examples](https://github.com/MODFLOW-ORG/modflow6-examples/blob/develop/scripts/ex-prt-mp7-p02.py) 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``\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sd = flopy.modpath.CellDataType()\n", "mp7_particle_data = flopy.modpath.NodeParticleData(subdivisiondata=[sd],\n", " nodes=list(particle_start_nodes))\n", "prt_particle_data = list(mp7_particle_data.to_prp(prt.modelgrid))\n", "prt_particle_data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Instantiate the MODFLOW 6 prt particle release point (PRP) package\n", "\n", "**Notes:**\n", "\n", "``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\n", "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\n", "considered to be the water table elevation (the head in the cell) rather than the top defined by the user. \n", "\n", "``istopzone``—integer value defining the stop zone number. If cells have been assigned IZONE values in the\n", "GRIDDATA block, a particle terminates if it enters a cell whose IZONE value matches ISTOPZONE. An\n", "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. \n", "\n", "See the MODFLOW 6 Description of Input and Output for more options\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "flopy.mf6.ModflowPrtprp(\n", " prt,\n", " nreleasepts=len(prt_start_data),\n", " packagedata=prt_start_data.to_records(index=False).tolist(),\n", " local_z=True,\n", " boundnames=True,\n", " perioddata={0: [\"FIRST\"]},\n", " exit_solve_tolerance=1e-5,\n", " extend_tracking=True,\n", " #istopzone=istopzone\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Instantiate the MODFLOW 6 PRT output control package\n", "\n", "**Notes:**\n", "PRT outputs a record each time a:\n", " 0: particle was released\n", " 1: particle exited a cell\n", " 2: time step ended\n", " 3: particle terminated\n", " 4: particle entered a weak sink cell\n", " 5: user-specified tracking time\n", " \n", "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." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "budget_record = [budgetfile_prt]\n", "track_record = [trackfile_prt]\n", "trackcsv_record = [trackcsvfile_prt]\n", "# track positions every year for 100 years\n", "track_nyears = 100\n", "tracktimes = np.linspace(0, track_nyears*365.25, track_nyears+1)\n", "flopy.mf6.ModflowPrtoc(\n", " prt,\n", " budget_filerecord=budget_record,\n", " track_filerecord=track_record,\n", " trackcsv_filerecord=trackcsv_record,\n", " ntracktimes=0,#len(tracktimes),\n", " tracktimes=None,#[(t,) for t in tracktimes],\n", " saverecord=[(\"BUDGET\", \"ALL\")],\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Instantiate the PRT Flow Model Interface and Explicit Model Solution Packages\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fmi_pd = [\n", " (\"GWFHEAD\", os.path.relpath(gwf_ws / headfile, prt_model_ws)), #Path(f\"../{gwf_ws.name}/{headfile}\")),\n", " (\"GWFBUDGET\", os.path.relpath(gwf_ws / budgetfile, prt_model_ws)) #Path(f\"../{gwf_ws.name}/{budgetfile}\")),\n", "]\n", "flopy.mf6.ModflowPrtfmi(prt, packagedata=fmi_pd)\n", "\n", "# Create an explicit model solution (EMS) for the MODFLOW 6 prt model\n", "ems = flopy.mf6.ModflowEms(\n", " prt_sim,\n", " filename=f\"{prt_name}.ems\",\n", ")\n", "sim.register_solution_package(ems, [prt.name])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Write the PRT input files and run the model" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "prt_sim.write_simulation()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "prt_sim.run_simulation()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Post-Process the MODFLOW and PRT Results\n", "\n", "\n", "### Load MODFLOW and PRT results from the heads and pathline files\n", "\n", "Load the MODFLOW heads" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "hobj = gwf.output.head()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "hds = hobj.get_data()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Load the tracking CSV file" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "particle_data = pd.read_csv(prt_model_ws / trackcsvfile_prt)\n", "particle_data.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plot the heads and pathlines" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(5, 9))\n", "mm = flopy.plot.PlotMapView(model=gwf, layer=0, ax=ax)\n", "mm.plot_array(hds, masked_values=[1e30])\n", "\n", "mm.plot_bc('SFR', color='b', plotAll=True)\n", "mm.plot_bc('WEL', plotAll=True)\n", "mm.plot_ibound()\n", "mm.plot_pathline(particle_data, layer='all', color='blue', lw=1)\n", "mm.plot_grid(lw=0.2, color=\"0.5\")\n", "\n", "ax = plt.gca()\n", "ax.fill(csx, csy, color='#e534eb', zorder=100, alpha=.75);" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(5, 9))\n", "\n", "mm = flopy.plot.PlotMapView(model=gwf, ax=ax, layer=0)\n", "mm.plot_array(hds, masked_values=[1e30])\n", "mm.plot_bc('SFR', color='b', plotAll=True)\n", "mm.plot_bc('WEL', plotAll=True)\n", "mm.plot_ibound()\n", "mm.plot_grid(lw=0.2, color=\"0.5\")\n", "\n", "df = particle_data\n", "vmin, vmax = df['t'].min(), df['t'].max()\n", "\n", "times = list(range(0, 74001, 1000))\n", "for ix in range(1, len(times)):\n", " tmp = df[(df['t'] >= times[ix - 1]) & (df['t'] < times[ix])]\n", " s = ax.scatter(tmp['x'].values, tmp['y'].values, s=5, c=tmp['t'].values, \n", " vmin=vmin, vmax=vmax, cmap=\"magma\")\n", " ax.set_title(f\"{times[ix - 1]} - {times[ix]} days\")\n", " if ix == 1:\n", " cbar = fig.colorbar(s, shrink=0.7)\n", " cbar.set_label('particle travel time from source', rotation=270, labelpad=14)\n", " display(fig)\n", " clear_output(wait=True)\n", " plt.pause(0.1) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Export the pathlines to a GeoPackage\n", "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)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "output_geopackage = Path(prt_model_ws / 'Pathlines.gpkg')\n", "output_geopackage.unlink(missing_ok=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### First the PRT results DataFrame needs to be converted to a ``GeoDataFrame``\n", "* In the real world, the model coordinates also need to be converted to real world coordinates. \n", "* We may also want to convert the travel times from the typical model time units of days to years.\n", "* We may want to decompose the ``cellid`` into model layer (`k`) and ``cell2d`` (location in each layer, in this case), or layer, row column if the model has a regular structured grid.\n", "* 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)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gwf.modelgrid.crs # no CRS attached" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gwf.modelgrid.crs = 26915" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "df['t_years'] = df['t'] * 325.25\n", "df['k'] = df['icell'].values[-1] // ncells_per_layer\n", "df['cell2d'] = df['icell'].values[-1] % ncells_per_layer\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x_crs, y_crs = gwf.modelgrid.get_coords(df['x'], df['y'])\n", "\n", "gdf = gpd.GeoDataFrame(df, \n", " geometry=gpd.points_from_xy(x_crs, y_crs),\n", " crs=gwf.modelgrid.crs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### write the GeoPackage\n", "Write particles terminating in Well and SFR cells to separate layers" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "for package, izone in izones.items():\n", " layer_name = f\"Particles going to {package.upper()} cells\"\n", " izone_df = gdf.loc[gdf['izone'] == izone]\n", " izone_df.to_file(output_geopackage, index=False, layer=layer_name)\n", "\n", "gdf.loc[gdf['izone'] == 0].to_file(\n", " output_geopackage, index=False, layer='All other particles')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Add a layer with pathlines\n", "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." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from shapely.geometry import LineString\n", "\n", "# create set of pathlines with start and end information\n", "# first group by particle\n", "by_particle = gdf.sort_values(by=['irpt', 't']).groupby('irpt')\n", "lines = by_particle.last()\n", "lines.crs = gdf.crs # CRS isn't retained in groupby\n", "lines.rename(columns={'k': 'end_k', 'cell2d': 'end_cell2d'}, inplace=True)\n", "line_starts = by_particle.first()\n", "lines['start_k'] = line_starts['k']\n", "lines['start_cell2d'] = line_starts['cell2d']\n", "\n", "# above we sorted the results by particle and then time\n", "# use the .agg method on the grouby object to create a LineString\n", "# from the (sorted) sequence of points defining each particle path\n", "linestring_geoms = by_particle['geometry'].agg(lambda x: LineString(x))\n", "lines['geometry'] = linestring_geoms" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "for package, izone in izones.items():\n", " layer_name = f\"Pathlines going to {package.upper()} cells\"\n", " izone_df = lines.loc[lines['izone'] == izone]\n", " izone_df.to_file(output_geopackage, index=False, layer=layer_name)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plot pathlines going to SFR and WEL cells" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(6, 10))\n", "pmv = flopy.plot.PlotMapView(gwf, ax=ax)\n", "pmv.plot_grid(lw=0.5)\n", "line_colors = {\n", " 'sfr': 'b',\n", " 'wel': 'r'\n", "}\n", "for package, izone in izones.items():\n", " izone_df = lines.loc[lines['izone'] == izone]\n", " izone_df.plot(ax=ax, ec=line_colors[package],\n", " label=f'Pathlines going to {package.upper()} cells')\n", "pmv.plot_ibound()\n", "pmv.ax.legend()" ] } ], "metadata": { "kernelspec": { "display_name": "pyclass", "language": "python", "name": "python3" } }, "nbformat": 4, "nbformat_minor": 4 }