Skip to content

Commit 476dd4c

Browse files
aleafaleaf
committed
feat(10a_prt_particle_tracking-demo.ipynb): add iflowface
1 parent 83cca9a commit 476dd4c

1 file changed

Lines changed: 219 additions & 11 deletions

File tree

‎notebooks/part1_flopy/10a_prt_particle_tracking-demo.ipynb‎

Lines changed: 219 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -133,7 +133,8 @@
133133
"cell_type": "markdown",
134134
"metadata": {},
135135
"source": [
136-
"### Change the workspace"
136+
"### Change the workspace\n",
137+
"(so that we don't modify the original flow model)"
137138
]
138139
},
139140
{
@@ -180,6 +181,172 @@
180181
"sim.write_simulation()"
181182
]
182183
},
184+
{
185+
"cell_type": "markdown",
186+
"metadata": {},
187+
"source": [
188+
"### Add `IFLOWFACE=6` information to the SFR Package\n",
189+
"This may not be necessary here because of the quadtree grid (that is refined in the stream cells), but it is often best practice.\n",
190+
"\n",
191+
"Particle tracking with numerical models is prone to the problem of weak sinks, where a sink (flux of water leaving a cell) is not strong enough to capture all of the flow through that cell. In other words, one or more cell faces still has outward flow, allowing water to pass through the cell. By default, PRT and MODPATH assume that sinks are distributed uniformly throughout a cell (`IFLOWFACE=0` or `IFACE=0`; see below). This setting is appropriate for wells, but for surface water features (especially small streams), this often prevents particles from discharging, leading to unrealistic pathlines and travel times. Some possible solutions:\n",
192+
"* reduce the grid size until all cells of interest are strong sinks\n",
193+
"* set `STOP_AT_WEAK_SINK=True` in the PRT Particle Release (PRP) input (or the equivalent setting in MODPATH), causing all particles to terminate at weak sinks. This option is often not appropriate for small streams, which in many cases should act as weak sinks (only capturing the shallowest particles).\\\n",
194+
"* set `IFLOWFACE=-1` or `IFACE=6` for cells with surface water boundary conditions, to instruct PRT or MODPATH to apply the sink discharge to the top face. In cells with small (weak) streams, assuming the stream bottom corresponds to the cell top, shallow particles that are near the top of the cell will discharge, while deeper particles will be allowed to pass through realistically. Abrams and others (2013) discuss this concept in more detail.\n",
195+
"\n",
196+
"**From the MODFLOW 6 I/O guide:**\n",
197+
"By default, flows associated with stress packages are assumed to be distributed uniformly over the volume of a cell.\n",
198+
"Distributed external inflows and outflows are reflected in the cell-cell flows calculated by the GWF Model, but they are\n",
199+
"not directly involved in determining the normal face velocities used to track a particle through the cell. The user can\n",
200+
"optionally assign a flow associated with a stress package to any face of the cell. Assignment of external flows is done by\n",
201+
"setting the value of an auxiliary package variable called IFLOWFACE to a value that corresponds to one of the cell faces.\n",
202+
"To accommodate non-rectangular cells, the IFLOWFACE numbering scheme in the PRT Model is based on clockwise\n",
203+
"ordering of the lateral cell faces. For a DIS-grid cell, IFLOWFACE = 1 corresponds to the “western” face, i.e., the face parallel to the y axis that has the lesser x coordinate. For a DISV-grid cell, IFLOWFACE = 1 corresponds to the face\n",
204+
"that extends in the clockwise direction from the first vertex listed for that cell in the CELL2D input block. IFLOWFACE\n",
205+
"numbering of lateral cell faces proceeds in clockwise direction. IFLOWFACE = -2 corresponds to the cell bottom, and\n",
206+
"IFLOWFACE = -1 corresponds to the cell top.\n",
207+
"\n",
208+
"**In MODPATH 7** (and earlier)\n",
209+
"`IFACE` is equivalent to `IFLOWFACE` in MODFLOW 6 PRT. Default `IFACE` values can be specified by package, or `IFACE` can be specified by cell (see the docs). A key difference is that `IFACE=6` corresponds to the top of the cell, and `IFACE=5` corresponds to the bottom."
210+
]
211+
},
212+
{
213+
"cell_type": "markdown",
214+
"metadata": {},
215+
"source": [
216+
"#### Two ways to add `IFLOWFACE`\n",
217+
"\n",
218+
"**1) Using Flopy**\n",
219+
"\n",
220+
"Add an auxiliary column named `'iflowface'` to the Options block"
221+
]
222+
},
223+
{
224+
"cell_type": "code",
225+
"execution_count": null,
226+
"metadata": {},
227+
"outputs": [],
228+
"source": [
229+
"from copy import deepcopy\n",
230+
"\n",
231+
"gwf2 = deepcopy(gwf)"
232+
]
233+
},
234+
{
235+
"cell_type": "code",
236+
"execution_count": null,
237+
"metadata": {},
238+
"outputs": [],
239+
"source": [
240+
"gwf2.sfr.auxiliary = ['iflowface']\n",
241+
"gwf2.sfr.auxiliary"
242+
]
243+
},
244+
{
245+
"cell_type": "markdown",
246+
"metadata": {},
247+
"source": [
248+
"Add `iflowface` to the package data"
249+
]
250+
},
251+
{
252+
"cell_type": "code",
253+
"execution_count": null,
254+
"metadata": {},
255+
"outputs": [],
256+
"source": [
257+
"packagedata = pd.DataFrame(gwf2.sfr.packagedata.array)\n",
258+
"# insert the auxiliary column before the existing boundname column\n",
259+
"# (so that the boundname column is still last)\n",
260+
"packagedata.insert(len(packagedata.columns)-1, 'iflowface', -1)\n",
261+
"# make sure the iflowface column is an integer!\n",
262+
"packagedata['iflowface'] = packagedata['iflowface'].astype(int)\n",
263+
"# cast packagedata back to a recarray (required by Flopy)\n",
264+
"# Note: index=False is needed here, \n",
265+
"# otherwise the extra index column inserted by pandas edges out the boundnames column\n",
266+
"# resulting in no boundnames being passed to Flopy\n",
267+
"gwf2.sfr.packagedata.set_data(packagedata.to_records(index=False))"
268+
]
269+
},
270+
{
271+
"cell_type": "code",
272+
"execution_count": null,
273+
"metadata": {},
274+
"outputs": [],
275+
"source": [
276+
"gwf2.sfr.packagedata.array[:5]"
277+
]
278+
},
279+
{
280+
"cell_type": "markdown",
281+
"metadata": {},
282+
"source": [
283+
"**Alternatively, we can write the model files, and then edit the SFR Package file directly.**"
284+
]
285+
},
286+
{
287+
"cell_type": "code",
288+
"execution_count": null,
289+
"metadata": {},
290+
"outputs": [],
291+
"source": [
292+
"write_iflowface = True\n",
293+
"set_iflowface = -1\n",
294+
"#original_sfr_filename = load_ws / gwf.sfr.filename\n",
295+
"sfr_filename = gwf_ws / gwf.sfr.filename\n",
296+
"original_sfr_filename = sfr_filename\n",
297+
"new_sfr_lines = list()\n",
298+
"package_already_has_iflowface = False\n",
299+
"aux_entry_added = False\n",
300+
"boundnames = False\n",
301+
"with open(original_sfr_filename) as src:\n",
302+
" for line in src:\n",
303+
" normalized_line = ' '.join(line.lower().split()).lower()\n",
304+
" if 'begin options' in normalized_line:\n",
305+
" new_sfr_lines.append(line)\n",
306+
" for line in src:\n",
307+
" normalized_line = ' '.join(line.lower().split()).lower()\n",
308+
" if 'boundnames' in normalized_line:\n",
309+
" boundnames = True\n",
310+
" if 'auxiliary' in normalized_line:\n",
311+
" if 'iflowface' not in normalized_line:\n",
312+
" items = line.strip().split()\n",
313+
" items.append('iflowface')\n",
314+
" new_sfr_lines.append(f\" {' '.join(items)}\\n\")\n",
315+
" else:\n",
316+
" new_sfr_lines.append(line)\n",
317+
" package_already_has_iflowface = True\n",
318+
" aux_entry_added = True\n",
319+
" elif not aux_entry_added and 'end options' in normalized_line:\n",
320+
" new_sfr_lines.append(' auxiliary iflowface\\n')\n",
321+
" new_sfr_lines.append(line)\n",
322+
" aux_entry_added = True\n",
323+
" break\n",
324+
" elif 'end options' in normalized_line:\n",
325+
" break\n",
326+
" else:\n",
327+
" new_sfr_lines.append(line)\n",
328+
" elif 'begin packagedata' in normalized_line:\n",
329+
" new_sfr_lines.append(line)\n",
330+
" for line in src:\n",
331+
" if 'end packagedata' in line.lower():\n",
332+
" new_sfr_lines.append(line)\n",
333+
" break\n",
334+
" items = line.strip().split()\n",
335+
" # if there's a boundname column, it needs to be last\n",
336+
" if boundnames:\n",
337+
" items.insert(-1, str(set_iflowface))\n",
338+
" else:\n",
339+
" items.append(set_iflowface)\n",
340+
" new_sfr_lines.append(f\" {' '.join(items)}\\n\")\n",
341+
" else:\n",
342+
" new_sfr_lines.append(line)\n",
343+
"\n",
344+
"if write_iflowface and not package_already_has_iflowface:\n",
345+
" with open(sfr_filename, 'w') as dest:\n",
346+
" for line in new_sfr_lines:\n",
347+
" dest.write(line)"
348+
]
349+
},
183350
{
184351
"cell_type": "markdown",
185352
"metadata": {},
@@ -431,9 +598,9 @@
431598
"outputs": [],
432599
"source": [
433600
"sd = flopy.modpath.CellDataType()\n",
434-
"mp7_particle_data = flopy.modpath.NodeParticleData(subdivisiondata=[sd],\n",
435-
" nodes=list(particle_start_nodes))\n",
436-
"prt_particle_data = list(mp7_particle_data.to_prp(prt.modelgrid))\n",
601+
"mp7_particle_data = flopy.modpath.NodeParticleData(subdivisiondata=sd,\n",
602+
" nodes=[cell2d for cell2d in particle_start_nodes])\n",
603+
"prt_particle_data = list(mp7_particle_data.to_prp(prt.modelgrid, localz=True))\n",
437604
"prt_particle_data"
438605
]
439606
},
@@ -465,6 +632,7 @@
465632
"flopy.mf6.ModflowPrtprp(\n",
466633
" prt,\n",
467634
" nreleasepts=len(prt_start_data),\n",
635+
" #packagedata=prt_particle_data,\n",
468636
" packagedata=prt_start_data.to_records(index=False).tolist(),\n",
469637
" local_z=True,\n",
470638
" boundnames=True,\n",
@@ -483,12 +651,13 @@
483651
"\n",
484652
"**Notes:**\n",
485653
"PRT outputs a record each time a:\n",
486-
" 0: particle was released\n",
487-
" 1: particle exited a cell\n",
488-
" 2: time step ended\n",
489-
" 3: particle terminated\n",
490-
" 4: particle entered a weak sink cell\n",
491-
" 5: user-specified tracking time\n",
654+
"\n",
655+
" 0) particle was released\n",
656+
" 1) particle exited a cell\n",
657+
" 2) time step ended\n",
658+
" 3) particle terminated\n",
659+
" 4) particle entered a weak sink cell\n",
660+
" 5) user-specified tracking time\n",
492661
" \n",
493662
"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."
494663
]
@@ -811,7 +980,9 @@
811980
"for package, izone in izones.items():\n",
812981
" layer_name = f\"Pathlines going to {package.upper()} cells\"\n",
813982
" izone_df = lines.loc[lines['izone'] == izone]\n",
814-
" izone_df.to_file(output_geopackage, index=False, layer=layer_name)"
983+
" izone_df.to_file(output_geopackage, index=False, layer=layer_name)\n",
984+
"izone_df.loc[izone_df['izone'] == 0].to_file(\n",
985+
" output_geopackage, index=False, layer='All other pathlines')"
815986
]
816987
},
817988
{
@@ -838,16 +1009,53 @@
8381009
" izone_df = lines.loc[lines['izone'] == izone]\n",
8391010
" izone_df.plot(ax=ax, ec=line_colors[package],\n",
8401011
" label=f'Pathlines going to {package.upper()} cells')\n",
1012+
"lines.loc[lines['izone'] == 0].plot(\n",
1013+
" ax=ax, ec='g',\n",
1014+
" label=f'All other pathlines')\n",
8411015
"pmv.plot_ibound()\n",
8421016
"pmv.ax.legend()"
8431017
]
1018+
},
1019+
{
1020+
"cell_type": "code",
1021+
"execution_count": null,
1022+
"metadata": {},
1023+
"outputs": [],
1024+
"source": [
1025+
"lines"
1026+
]
1027+
},
1028+
{
1029+
"cell_type": "markdown",
1030+
"metadata": {},
1031+
"source": [
1032+
"### References\n",
1033+
"Abrams, D., Haitjema, H. and Kauffman, L. (2013), On Modeling Weak Sinks in MODPATH. Groundwater, 51: 597-602. https://doi.org/10.1111/j.1745-6584.2012.00995.x"
1034+
]
1035+
},
1036+
{
1037+
"cell_type": "markdown",
1038+
"metadata": {},
1039+
"source": []
8441040
}
8451041
],
8461042
"metadata": {
8471043
"kernelspec": {
8481044
"display_name": "pyclass",
8491045
"language": "python",
8501046
"name": "python3"
1047+
},
1048+
"language_info": {
1049+
"codemirror_mode": {
1050+
"name": "ipython",
1051+
"version": 3
1052+
},
1053+
"file_extension": ".py",
1054+
"mimetype": "text/x-python",
1055+
"name": "python",
1056+
"nbconvert_exporter": "python",
1057+
"pygments_lexer": "ipython3",
1058+
"version": "3.14.3"
8511059
}
8521060
},
8531061
"nbformat": 4,

0 commit comments

Comments
 (0)