User's Guide for mod-PATH3DU

A groundwater path and travel-time simulator

  • Christopher Muffels
  • Leland Scantlebury
  • Xiaomin Wang
  • Matthew Tonkin
  • Christopher Neville
  • Muhammad Ramadhan
  • James R. Craig
  • Bethesda, Maryland
  • Waterloo, Ontario
  • Last Update: November, 2018

Disclaimer

S.S. Papadopulos & Associates Inc. (SSP&A) and the University of Waterloo (UW) make no representations or warranties with respect to the contents hereof and specifically disclaim any implied warranties of fitness for any particular purpose. Further, SSP&A and UW reserve the right to revise this publication and software, and to make changes from time to time in the content hereof without obligation of SSP&A or UW to notify any person of such revisions or changes.


Support

S.S. Papadopulos & Associates Inc. in collaboration with the University of Waterloo developed mod-PATH3DU to support specific project applications, but is distributing it free-of-charge as a service to our industry. As such, only limited technical support can be provided by S.S. Papadopulos & Associates Inc. through the corresponding author:

Christopher Muffels
7944 Wisconsin Ave.
Bethesda, MD 20895
301-718-8900
cmuffels@sspa.com
www.sspa.com


Acknowledgements

The authors would like to thank the following for their support in developing mod-PATH3DU. A little goes a long way and it is greatly appreciated.

  • Michael Kulbersh and the New England District Office of the U.S. Army Corps of Engineers (Contract No. W912WJ-16-P-0040, and W912WJ-16-D-0005).
  • The team at Montgomery & Associates, especially Derek Blazer, Colin Kikuchi and Timothy Bayley.

The authors acknowledge the following interns for their contributions to the development of mod-PATH3DU.

  • Graham Stonebridge (University of Waterloo) developed the proof-of-concept application of the SSP&A tracking method to MODFLOW,
  • Karl Mihm (University of Maryland) adapted Graham's code to MODPATH v.6 which became mod-PATH3DU v.1.0.0,
  • Robert Gordon (Drexel University) optimized a data tree structure to efficiently determine the model cell corresponding to a particle's location which substantially improved mod-PATH3DU's run time, and
  • Tim Wilson (Drexel University) in researching, developing, and testing a random walk methodology that will be included in a future release.

The first users of mod-PATH3DU were the developers and product leads for the top MODFLOW GUIs - Groundwater Vistas, GMS, and Visual MODFLOW Flex. Their support and encouragement provided a lot of the motivation to advance mod-PATH3DU, especially in the development of this second version. We thank the following from the different teams.

Finally, we acknowledge all our colleagues at S.S. Papadopulos & Assoc. Inc., especially Chunmiao Zheng, Charlie Andrews, Steve Larson, Marinko Karanovic, Jim Cousins and Doug Hayes. This software is a company wide effort drawing upon years, decades even, of research and testing in various forms and applications in support of our clients.

Quick Start Guide

  1. Download mod-PATH3DU and its support programs from www.sspa.com/software
  2. Read this manual in its entirety :)
  3. But seriously, read the IFACE section of this manual. IFACE controls particle capture. If IFACE is not set correctly mod-PATH3DU cannot terminate particles.
  4. If using MODFLOW, ensure all boundary packages write to a compact composite budget file. The "COMPACT BUDGET" option is specified in the MODFLOW output control (OC) file. Check the MODFLOW manual for more details. mod-PATH3DU will terminate in error if the compact budget option is not specified.
  5. Use writeP3DGSF.exe, a support program provided with mod-PATH3DU, to create/convert the grid-specification (GSF) file if necessary.
  6. Create the required mod-PATH3DU input files (Input Instructions):
    • Primary file - download an example from the Input Instructions and modify as needed following the provided instructions.
    • Prepare the Per-cell property file as needed.
    • Prepare starting locations point shapefile. Ensure it has attributes for starting location cell id, release time, and local-z elevation.
  7. Ensure all the model files and mod-PATH3DU input files and executable are in the same folder or directory.
  8. At a command prompt, change to the folder with the files and enter the name of the mod-PATH3DU executable, followed by the name of the primary input file. For example:
    mp3du.exe NAME.json colorcode
  9. Review mod-PATH3DU.log and MP3DU_log.lst to ensure the program executed as expected.
  10. Convert the binary output file into a table or shapefile using writeP3DOutput.exe, a support program provided with mod-PATH3DU.

Release Notes & Compatibility

Version 2.0.0

  • MODFLOW Compatibility
    • MODFLOW-2000, MODFLOW-2005, and MODFLOW-USG (v.1.4 and later)
    • Single or double precision
    • Multiple cell-by-cell flow files
    • MNW1, MNW2, and GLO package support
    • New PATH file type in MODFLOW-USG v.1.4
  • Parallelized
    • MODFLOW binary output file processing - this really speeds up processing of transient models
    • Particle paths
  • New, easy to read input file format (not backward compatible with mod-PATH3DU v.1.x.x)
  • Linear sorption via a retarded velocity
  • Dispersion via Random Walk, for enhanced visualization of transport processes

Overview

mod-PATH3DU is a particle tracking code for calculating the three-dimensional flow pathlines and travel time of solute particles. It supports both MODFLOW-USG (Panday et al., 2013) and MEUK (Tonkin et al., 2016). MODFLOW-USG is an unstructured-grid version of MODFLOW, particles are tracked within the flow solution calculated by MODFLOW models using the Waterloo method (Ramadhan, 2015). The Waterloo method is a semi-analytical approach, akin to the Pollock method (Pollock, 1988), that is based on an analytical Taylor series reconstruction of the intra-cell velocity. It is applicable to cells with arbitrary geometry, including quadpatch, quadtree, nested, and Voronoi grids. MEUK (multi-event universal kriging) is a method of interpolation rooted in universal kriging with hydrogeology-specific drift terms to map hydraulic-head. Particles are tracked on MEUK-calculated head surfaces using the SSP&A method (Tonkin and Larson, 2002). The SSP&A method evaluates the velocity components (x and y-directions) of a particle using Darcy's law, where the gradient in the vicinity of the particle is determined from the MEUK interpolated heads. mod-PATH3DU uses a higher order numerical Runge-Kutta scheme (Zheng, 1992) to move a particle.

This document describes the mod-PATH3DU particle tracking program and provides instructions for using the program. Some of the required input files are created outside of MODFLOW-USG and it is important to read the instructions to make sure that all the required input is specified. mod-PATH3DU does not include a graphical user interface (GUI). It is recommended that it be run from a command prompt. Users will need to devote care to preparing the input and analyzing the output. All the files for the test cases detailed in the Examples section are available for download with the program. These files illustrate the structure of the input files and the procedures required to execute mod-PATH3DU.

Compatibility

MODFLOW

  • MODFLOW-USG v.1.4 or later
  • MODFLOW-2005
  • MODFLOW-2000

mod-PATH3DU uses information in both the cell-by-cell flow (CBB) and head-save (HDS) files written by MODFLOW to calculate velocity. These files list the components of flow and hydraulic head for each model cell. In structured grid mode, MODFLOW-USG writes these two files in the same format as MODFLOW-2000 and MODFLOW-2005; hence, mod-PATH3DU supports these versions of MODFLOW as well.

MODFLOW-USG does not require information about cell shapes or how cells are positioned in space (Panday et al. 2013). Spatial information about cells, however, is critical to a particle tracking model. A grid specification file (GSF), supported by the popular MODFLOW GUIs, was decided upon that provides x, y and z coordinates for the vertices of a cell. mod-PATH3DU requires a modified version of this file for both structured and unstructured grid models. A pre-processing program, writeP3DGSF.exe, is included with mod-PATH3DU to convert a GSF file in the original format into the mod-PATH3DU specific format. Additionally, it will write the required GSF file for a structured MODFLOW model defined by DELR and DELC. The format for the modified GSF file is provided in the Input Instructions section. This modified GSF format makes it easier for mod-PATH3DU to identify the shared face between two cells that are connected and is included as a pre-processing step to improve runtime as this identification need only be done once for a model grid.

MEUK

The hydraulic-head maps written by MEUK are in the ArcMap/ArcInfo ASCII grid format (lower-left corner designation for the origin) - mod-PATH3DU uses these maps, as well as information about the grid, hydraulic conductivity, porosity, and events and drift terms used in MEUK, to track particles. It supports the multi-event feature of MEUK by tracking on each event's head surface for a specified duration, which is similar to the treatment of transient MODFLOW models as a series of steady-state flow periods. Particles are stopped or considered captured if they travel within a user-specified distance of an extraction or injection well when forward or backward tracking, respectively, and when they cross a line-sink. The MEUK grid must be written to a GSF file prior to executing mod-PATH3DU; the pre-processing program writeP3DGSF.exe can be used to do this. Hydraulic conductivity and porosity can be specified as constant across the entire MEUK domain, or be input from a file with non-uniform values. MEUK relies on a series of keywords to indicate the event and group a drift or observation set belong to, as well as for coordinates and other information required to calculate its maps; mod-PATH3DU uses the same keywords. The Input Instructions detail how to relate MEUK keywords to mod-PATH3DU.

Background

Advection

A chapter in Zheng and Bennett (2002) provides an excellent discussion on advective transport, including particle tracking solution methods. We summarize, in part, this chapter here. It is recommended that users interested in additional detail or information review the relevant chapter in this book.

The pathlines of purely advective solute particles are governed by

\begin{equation} \frac{d \mathbf X_p} {dt} = \mathbf V(\mathbf X_p,t) \label 1 \end{equation}

where \(\mathbf X_p(t)\) is the position of the particle at time \(t\), and \(\mathbf V(\mathbf X_p,t)\) is the seepage velocity of the particle at position \(\mathbf X_p\) and time \(t\). The solution to () for a particle at time \(t\) is

\begin{equation} \mathbf X_p(t) = \mathbf X_p(t_0) + \int_{t_0}^t \mathbf V(\mathbf X_p,t)dt \end{equation}

where \(\mathbf X_p(t_0)\) is the position at time \(t_0\). Equation () can be integrated directly, but is typically solved numerically given the complexity of velocity fields in most groundwater systems. The Taylor series expansion of () is given by:

\begin{equation} \mathbf X_p(t + \Delta t) = \mathbf X_p(t) + \frac{d \mathbf X_p} {dt} (\mathbf X_p,t) \Delta t + \frac{d^2 \mathbf X_p} {dt^2} (\mathbf X_p,t) \Delta t^2 + \ldots \end{equation}

where \(\mathbf X_p(t+\Delta t)\) is the position of the particle at the start of the next tracking step, \(t+\Delta t\). Euler's method,

\begin{equation} \mathbf X_p(t + \Delta t) = \mathbf X_p(t) + \mathbf V(\mathbf X_p,t) \Delta t \end{equation}

is a simple first-order solution to (). It relies only on the velocity at the starting point of each tracking step, which requires small tracking steps to be accurate. Accuracy can be increased by using higher order schemes like Runge-Kutta methods.

Runge-Kutta Numerical Scheme

The fourth-order Runge-Kutta method.

The fourth-order Runge-Kutta method (RK4), , solves the particle tracking equation by taking several Euler-like trial steps, \( \mathbf P_2 \), \( \mathbf P_3 \), and \( \mathbf P_4 \), where

\begin{equation} \begin{aligned} \mathbf P_2 & = \mathbf X_p(t) + \mathbf V \biggr (\mathbf X_p,t \biggr ) \frac { \Delta t } { 2 } \\ \mathbf P_3 & = \mathbf X_p(t) + \mathbf V \biggr (\mathbf P_2,t + \frac { \Delta t } { 2 } \biggr ) \frac { \Delta t } { 2 } \\ \mathbf P_4 & = \mathbf X_p(t) + \mathbf V \biggr (\mathbf P_3,t + \frac { \Delta t} {2} \biggr )\Delta t \end{aligned} \end{equation}

The velocity at each trial point, \( \mathbf V(\mathbf X_p,t) \), \( \mathbf V(\mathbf P_2,t+\Delta t/2) \), \( \mathbf V(\mathbf P_3,t+\Delta t/2) \), and \( \mathbf V(\mathbf P_4,t+\Delta t) \), respectively, are combined, as a weighted sum,

\begin{equation} \mathbf V_{RK4}(\mathbf X_p,t) = \frac {1}{6} \Biggr [ \mathbf V(\mathbf X_p,t) + 2 \mathbf V \biggr (\mathbf P_2,t + \frac { \Delta t } { 2 } \biggr ) + 2 \mathbf V \biggr (\mathbf P_3,t + \frac { \Delta t } { 2 } \biggr ) + \mathbf V \biggr (\mathbf P_4,t + \Delta t \biggr ) \Biggr ] \end{equation}

to yield an effective velocity, \( \mathbf V_{RK4}(\mathbf X_p,t) \), that results in a fourth-order accurate tracking step:

\begin{equation} \mathbf X_p(t + \Delta t) = \mathbf X_p(t) + \mathbf V_{RK4}(\mathbf X_p,t) \Delta t \end{equation}

The velocity components are calculated using either the Waterloo or SSP&A methods.

The accuracy of the RK4 method depends on the tracking stepsize: if it is too large, the calculated particle path may deviate from the actual flow path, and if it is too small, unnecessary steps are taken, which reduces performance. A step doubling technique (Press et al., XXXX) is implemented in mod-PATH3DU. It continuously evaluates the error in the calculated paths and updates the stepsize to maintain a user-specified level of accuracy.

Adaptive Stepsize Control

The adaptive stepsize control procedure.

With "step doubling", , the tracking step is taken as a full step, denoted as \( \mathbf X_1(t+\Delta t) \), and then again as two half-steps, the final endpoint denoted as \( \mathbf X_2(t+\Delta t)\), where each step is calculated with the Runge-Kutta method. Because the Runge-Kutta method employed is fourth order, the exact solution, denoted by \( \overline {\mathbf X_p}(t+\Delta t) \), is related to the two approximate solutions by:

\begin{align} \overline {\mathbf X_p}(t+\Delta t) &= \mathbf X_1(t+\Delta t) + {\Delta t}^5 \phi + O \bigr( \Delta t^6 \bigr ) + \ldots \\ \overline { \mathbf X_p}(t+\Delta t) &= \mathbf X_2(t+\Delta t) + 2 \biggr ( \frac {\Delta t} {2} \biggr )^5 \phi + O \bigr( \Delta t^6 \bigr ) + \ldots \end{align}

where \( \Delta t^5 \phi \) and \( 2 \bigr ( \frac {\Delta t} {2} \bigr )^5 \phi \) are the fifth-order error terms that are incurred taking the full step and two half-steps, respectively, and \( \phi \), to order \( {\Delta t}^5 \), is constant over the step. The distance between the two endpoints,

\begin{equation} \mathbf \Delta_s = \mathbf X_2(t+\Delta t)-\mathbf X_1(t+\Delta t) \end{equation}

is an estimate of the error in taking a step of size \( \Delta t \). It is a vector, with elements in each of the principal directions:

\begin{equation} \mathbf \Delta_s= \begin{bmatrix} \Delta X_s \\ \Delta Y_s \\ \Delta Z_s \end{bmatrix} \end{equation}

where \(\Delta X_s \), \(\Delta Y_s \), and \(\Delta Z_s \) are the distance between the full step and two half-steps in the x, y, and z directions, respectively. Further, by rearranging () and () and substituting into (), it is recognized that:

\begin{equation} \mathbf \Delta_s = 15 \cdot {\Biggr ( \frac {\Delta t} {2} \Biggr )}^5 \phi \end{equation}

The allowable error is

\begin{equation} \mathbf \Delta_\varepsilon= \begin{bmatrix} \Delta X_\varepsilon \\ \Delta Y_\varepsilon \\ \Delta Z_\varepsilon \end{bmatrix} \end{equation}

where

\begin{equation} \begin{aligned} \Delta X_\varepsilon = \varepsilon \cdot X_{max} \\ \Delta Y_\varepsilon = \varepsilon \cdot Y_{max} \\ \Delta Z_\varepsilon = \varepsilon \cdot Z_{max} \end{aligned} \end{equation}

The user specifies \(\varepsilon\). It is a single tolerance level that is applied to relative measures in the x, y, and z directions. These measures are the total length, \(X_{max}\), \(Y_{max}\), and \(Z_{max}\), respectively, of the flow domain in each direction. With this approach, the relative scale of each dimension is accounted for. Below is a portion of the Primary Input file indicating how a user provides the error criterion.

{
	.
	.
	.
	
"SIMULATIONS"
: [ {
"PATHLINE"
: { . . .
"ADAPTIVE_STEP_ERROR"
: 1.0e-6, . . . } ] }

The fifth-order error associated with moving a particle a distance of \(\mathbf \Delta_\varepsilon\) is:

\begin{equation} \mathbf {\Delta}_{\varepsilon}= {\biggr( {\Delta t_{\varepsilon}} \biggr)}^5 \phi \end{equation}

Taking the ratio of () to (),

\begin{equation} \frac {\mathbf {\Delta}_{\varepsilon}} {\mathbf \Delta_S} \approx \Biggr ( \frac {\Delta t_{\varepsilon}}{\Delta t} \Biggr )^5 \end{equation}

yields an estimate of the maximum allowable stepsize that satisfies the user-specified stability criterion:

\begin{equation} \Delta t_{\varepsilon} = f_s {\Delta t} \Biggr ( \frac {\mathbf {\Delta}_{\varepsilon}} {\mathbf \Delta_S} \Biggr )^{\frac{1}{5}} \end{equation}

() is evaluated every tracking step to estimate the required stepsize adjustment. It is evaluated in each of the three primary directions, and set to the shortest stepsize. If the tracking error exceeds the allowable tolerance, the stepsize is reduced to \( \Delta t_{\varepsilon} \), and the tracking step is repeated. If the tracking error is less than the specified tolerance, a step of \( \Delta t \) is taken, and the stepsize for the next tracking step is increased to the maximum of \( 4 \cdot \Delta t \) and \( \Delta t_{\varepsilon} \). A safety factor, \(f_s\), less than unity, is incorporated into the timestep adjustment because the error in () is approximate. The safety factor is fixed at 0.9 in mod-PATH3DU. The series of calculations is repeated to move a particle, step-by-step, until it reaches a boundary.

With "step doubling", the updated particle position is taken as the endpoint of the two half-steps. This position is corrected, however, using the full-step by recognizing, from (), that the fifth-order error of the two half-steps is \( \mathbf \Delta_S / 15 \); thus

\begin{equation} \mathbf X_p(t+\Delta t) = \mathbf X_2(t+\Delta t) + \frac {\mathbf \Delta_S} {15} \end{equation}

This estimate of the particle position is fifth-order accurate.

Advection and Dispersion

The dispersion capabilities of mod-PATH3DU are still under development. It has been tested and verified for simple 2D homogeneous flow models, however, dispersion in the z direction has not been verified. The primary motivation for developing a dispersion option in mod-PATH3DU is as a screening level tool and enhanced visualization of transport processes without the use of a traditional groundwater flow or transport model. It is for this reason that the focus has been on 2D dispersion and its application to interpolated hydraulic head surfaces calculated with MEUK.

In one dimension, the advection-dispersion equation is (e.g. Zheng and Bennett, 2002)

\begin{equation} \frac {\partial C}{\partial t}=\frac {\partial}{\partial x} \left(D \frac {\partial C}{\partial x}\right)- \frac {\partial}{\partial x}\left(vC\right) \end{equation}

where \(C\) is concentration, \(D\) is dispersion, and \(v\) is velocity. The analytical solution to (), following an instantaneous injection of mass at \(x=0\) and assuming uniform flow, is:

\begin{equation} C \left( x,t \right) = \frac {C_0}{\sqrt{4 \pi Dt}} e ^ {-(x-vt)^2/4Dt} \end{equation}

As Prickett et al. (1981) and others note, the spreading of solutes by dispersion can be described by a Gaussian or normal distribution, whose mean, \(\mu\), is the advective movement of a particle and whose standard deviation, \(\sigma\), is the dispersive movement; the density function of which is:

\begin{equation} f(x)=\frac {1}{2 \pi \sigma ^2}e^{-(x- \mu )^2 / 2 \sigma ^2} \end{equation}

Equating () and () is the foundation of the "Random Walk" solute transport model described by Prickett et al. (1981). However, flow in most natural groundwater systems is not uniform, as such, the Prickett et al. (1981) model is only approximate.

Kinzelbach (1988) and others demonstrate that for non-uniform flow fields, specifically when dispersion is not constant, the appropriate density function that satisfies (), for many particles, satisfies the Fokker-Planck equation. The Fokker-Planck equation in one-dimensional form is:

\begin{equation} \frac {\partial f}{\partial t}=\frac {\partial}{\partial x} \left( D \frac {\partial f}{\partial x} \right) - \frac {\partial}{\partial x} \left[ \left( v- \frac {\partial D}{\partial x} \right)f \right] \end{equation}

Equations () and () are not identical if the dispersion coefficients are space dependent. Adding a correction term

\begin{equation} v \to v+\frac{\partial D}{\partial x}=v^\prime \end{equation}

to the velocity in (), however, does make these two equations equivalent. Thus, to correctly solve the advection-dispersion equation with the random walk method, the velocity term is the seepage velocity plus the spatial derivative of the dispersion coefficient (Zheng and Bennett, 2002).

Generally, solute particles are advected at the average groundwater velocity and dispersed randomly, per (e.g. Salamon, 2006)

\begin{equation} \mathbf X_p(t+ \Delta t)=\mathbf X_p(t)+\mathbf A(\mathbf X_p,t) \Delta t + \mathbf B(\mathbf X_p,t)\mathbf Z \sqrt{\Delta t} \end{equation}

where \(\mathbf {BZ}\) is the (random) displacement due to dispersion, \(\mathbf B\) is the displacement matrix and is related to dispersion, \(\mathbf D\), by

\begin{equation} 2 \mathbf D=\mathbf {BB}^T \end{equation}

Dispersion, in three-dimensions, is defined following Burnett and Frind (1987) as

\begin{equation} \mathbf D = \frac {1}{\left\vert v \right\vert} \begin{bmatrix} \left( \alpha_l {v_x}^2 + \alpha_t^h {v_y}^2+ \alpha_t^v {v_z}^2 \right) & v_x v_y (\alpha_l-\alpha_t^h ) & v_x v_z (\alpha_l-\alpha_t^v ) \\\ v_x v_y (\alpha_l-\alpha_t^h ) & (\alpha_t^h {v_x}^2+\alpha_l {v_y}^2+\alpha_t^v {v_z}^2 ) & v_y v_z (\alpha_l-\alpha_t^v ) \\\ v_x v_z (\alpha_l-\alpha_t^v ) & v_y v_z (\alpha_l-\alpha_t^v ) & \alpha_t^v ({v_x}^2+{v_y}^2 )+\alpha_l {v_z}^2 \end{bmatrix} \end{equation}

Thus,

\begin{equation} \mathbf B = \begin{bmatrix} \frac {v_x \sqrt { 2 \left (\alpha_l |v| \right) }}{\left\vert v \right\vert} & - \frac {v_x v_z \sqrt {2 \left (\alpha_t^v \right) \left \vert v \right \vert}}{\left \vert v \right \vert \sqrt {{v_x}^2 + {v_y}^2}} & {- \frac {v_y}{\sqrt {{v_x}^2 + {v_y}^2}}} {\sqrt {2 \left( \frac {\alpha_t^h \left({v_x}^2+{v_y}^2\right)+\alpha_t^v{v_z}^2}{\left\vert v \right\vert} \right) } } \\\ \frac {v_y \sqrt { 2 \left (\alpha_l |v| \right) }}{\left\vert v \right\vert} & - \frac {v_y v_z \sqrt {2 \left (\alpha_t^v \right) \left \vert v \right \vert}}{\left \vert v \right \vert \sqrt {{v_x}^2 + {v_y}^2}} & {- \frac {v_x}{\sqrt {{v_x}^2 + {v_y}^2}}} {\sqrt {2 \left( \frac {\alpha_t^h \left({v_x}^2+{v_y}^2\right)+\alpha_t^v{v_z}^2}{\left\vert v \right\vert} \right) } } \\\ \frac {v_z \sqrt {2 \left(\alpha_l \left\vert v \right\vert \right)}}{\left\vert v \right\vert} & \sqrt { \frac {2 \alpha_v^t \left\vert v \right\vert \left( {v_x}^2+{v_y}^2 \right) } { {\left\vert v \right\vert}^2 } } & 0 \end{bmatrix} \end{equation}

And, the velocity correction, (), termed the "drift" vector, \(\mathbf A \left( \mathbf X_p,t \right)\), is

\begin{equation} \mathbf A \left( \mathbf X_p,t \right)=\mathbf V \left( \mathbf X_p,t \right) + \nabla \cdot \mathbf D \left( \mathbf X_p,t \right) \end{equation}

mod-PATH3DU uses the solution to () given by Labolle et al. (2000). Their generalized stochastic differential solution (GSDE) assumes constant porosity and isotropic dispersion and is given by the equation

\begin{equation} \mathbf X_p \left( \mathbf t + \Delta t \right) = \mathbf X_p \left( t \right) + \mathbf V \left( \mathbf X_p,t \right) \Delta t + \mathbf B \left( \mathbf X_p + \Delta \mathbf Y,t \right) \mathbf Z \sqrt {\Delta t} \end{equation}

where

\begin{equation} \Delta \mathbf Y = \mathbf B \left( \mathbf X_p,t \right) \mathbf Z \sqrt {\Delta t} \end{equation}

Implementation of this approach is a two-step process. First, () is evaluated to determine the particle displacement due to dispersion, and second, this point is used to determine the dispersion displacement in (). The velocity is calculated using either the SSP&A method or the Waterloo method. The velocity calculated by the Waterloo method is discontinuous at the interface between two cells which is problematic for the random walk method when evaluating the drift vector (e.g. Salamon et al. 2006). This issue is illustrated in (Example 4b). Future work will evaluate replacing the Waterloo method in the drift calculation with the SSP&A method, which is approximate, but continuous across cell interfaces.

For MODFLOW models, the dispersivity parameters are specified in the Per-cell Property File, while for MEUK models they are specified in the Primary File. Dispersion is not simulated until the "DISPERSION" option is specified in the Primary File, a portion of this file is given below to indicate where this option is specified.

{
	.
	.
	.
	
"SIMULATIONS"
: [ {
"PATHLINE"
: { . . .
"OPTIONS"
: [
"DISPERSION"
], } . . . } ] }

Solute Retardation

Adsorption according to a linear isotherm is incorporated into mod-PATH3DU by replacing the pore velocity with a retarded velocity,

$$ \frac {\mathbf V(\mathbf X_p,t)}{R} $$

everywhere, i.e. (), (), and ().

For MODFLOW models, the retardation factor, \(R\), is specified in the Per-cell Property File, while for MEUK models it is specified in the Primary File. Even if a retardation factor is specified, adsorption is not simulated until the "RETARDATION" option is specified in the Primary File, a portion of this file is given below to indicate where this option is specified.

{
	.
	.
	.
	
"SIMULATIONS"
: [ {
"PATHLINE"
: { . . .
"OPTIONS"
: [
"RETARDATION"
], } . . . } ] }

The Waterloo Method

The Waterloo method is detailed by (Ramadhan, 2015). It is based upon reconstructing the intra-cell velocity field using local analytical solutions to the Laplace and Poisson equations for each cell. Solutions to Laplace's equation can be formulated in terms of a complex potential function, \(\Omega\) (e.g. Haitjema, 2005):

\begin{equation} \Omega \left( z \right) = \Phi \left( z \right) + i \Psi \left( z \right) \end{equation}

The real and imaginary parts of \(\Omega\) are the discharge potential \(\Phi\), and the stream function \(\Psi\), respectively, and \(z\) is the complex coordinate, where \(z=x+iy\). The Waterloo method makes use of complex discharge potential \(W\), defined as (e.g. Haitjema, 1995; Ramadhan, 2015):

\begin{equation} W = - \frac {d \Omega}{dz} = - \frac {\partial \Phi}{\partial x} - i \frac {\partial \Psi}{\partial x} = Q_x - i Q_y \end{equation}

where \(Q_x\) and \(Q_y\) are the volumetric discharge (\(L^3/t\)) in the x- and y- directions, respectively. A special property of complex coordinates is that any function of \(z\), \(f(z)\), has real and imaginary parts that automatically satisfy the Laplace equation. The Waterloo method, for example, uses two such functions, the complex Taylor series solution of a circular inhomogeneity (Jankovic and Barnes, 1999b) to account for flow through the sides of a cell, and the complex logarithmic function to represent flow to a point sink or source. These functions are respectively:

\begin{equation} f \left( z \right) = \sum_{n=0}^\infty {a_nz^n} \end{equation}

and

\begin{equation} f \left( z \right) = \frac {Q}{2 \pi} \ln \left( z \right) \end{equation}

where \(a_n\) are unknown complex coefficients. These functions form the basis by which the Waterloo method calculates a local exact solution to the groundwater flow equation for any cell. The form of the solution in terms of a complex potential is (Ramadhan, 2015):

\begin{equation} \Omega \left( Z \right) = \Phi \left( Z \right) + i \Psi \left( Z \right) = \sum_{n=0}^N {a_n Z^n - \frac {q_v}{2} \Re {\left( RZ \right)}^2 + \frac {Q_w}{2 \pi} \log { \left( R \left \vert Z - Z_w \right \vert \right) } } \end{equation}

where \(Z=\frac{z-z_c}{R}\) is the local complex coordinate of the particle, \(z\), relative to the complex location of the cell center \(z_c\), \(R\) is the maximum radius of a circle centered at the cell center that fully encompasses the cell, \(N\) is termed the order of approximation (Jankovic and Barnes, 1999a), \(q_v\) is the vertical flow to the cell from recharge, creeks or streams, cells above or below, etc., \(Q_w\) is the pumping rate of a well in the cell located at \(Z_w=\frac{z_w-z_c}{R}\). The only unknowns in () are the \(N\) \(a\) parameters. These parameters are solved for using the overspecification principal discussed by Jankovic and Barnes (1999a and 1999b).

The SSP&A Method

The SSP&A method was originally developed by Tonkin and Larson (2002) to track particles and estimate hydraulic capture on a mapping of ground water level data under the influence of pumping. Later it was used by Karanovic et al. (2009) to calculate capture frequency maps. It is a grid independent approach and therefore, sufficiently flexible to be used with either structured or unstructured grids. The SSP&A method evaluates the velocity components of a particle in the x and y directions using the interpolated head based on the distribution solved for by MEUK, for example. The head interpolation is implemented at points in the immediate vicinity of a particle using universal kriging (), from which the velocity components in the x and y directions are calculated according to Darcy's law.

Schematic of SSP&A method velocity calculation

The SSP&A method uses bilinear interpolation, by default, to determine these gradients, except near a pumping well, where kriging is awlays used. The user can choose to use kriging everywhere.

Kriging is often used to interpolate irregularly spaced measurement data to unsampled locations (typically a grid of points suitable for contouring). The values at unknown locations are then estimated by minimizing the error variance of predicted values based on their spatial distribution. In the context of mod-PATH3DU, the heads calculated by MEUK that are used to determine velocity constitute the "measured" data, while the position of a particle is the "unsampled" location. Kriging is an exact interpolator in the absence of measurement error or co-located data. Two popular forms of kriging are simple and ordinary. With simple kriging, the mean of the data is assumed to be constant everywhere and its value known a-priori. With ordinary kriging, the mean is assumed to be unknown, but is estimable using some function of the measured data. In addition, ordinary kriging can support a spatially varying mean that is not only a function of the data, but includes some trend. This form of ordinary kriging is called universal kriging. Universal kriging is a means of incorporating the shape associated with different hydrologic features into the estimate at an unsampled location. For example, near a pumping well, a point sink/source of known strength, derived from the Theis equation, can be used to incorporate the logarithmic drawdown shape (Tonkin and Larson, 2002). The mathematical details on kriging have been well documented by Journel and Huijbregts (1992), Cressie (1993) and Deutsch and Journel (1998)[](). The kriging algorithm used in the SSP&A method is adapted from the public domain software Geostatistical Software Library (GSLIB) (Deutsch and Journel, 1998) and Skrivan and Karlinger (1980).

Using mod-PATH3DU

IFACE and Particle Termination

Like the Pollock method, the Waterloo method computes the velocity vector using the flow components of the cell faces (sides, top and bottom) calculated by MODFLOW. The Waterloo method is predicated on a net flow of zero for each cell. It is important that, in addition to the flux across each cell face, boundary condition flows are assigned to a face as well - which face being dictated by the purpose of the boundary condition. Table 1 lists the MODFLOW boundary packages supported by mod-PATH3DU. The auxiliary variable IFACE (Pollock, 2012) is the mechanism by which boundary flows are assigned to a particular cell face and included in the correct term in Equation . Assigning ALL boundary package flows to the correct term on the right-hand side (RHS) of Equation via IFACE is very important. For example, consider recharge, its flux is vertical and enters a cell through the top face. If the recharge flux is not assigned to the top face, the Waterloo method will not account for it in the vertical flow term of Equation and erroneous particle paths will be calculated or the program will terminate in error because the unknown coefficients in Equation cannot be solved for. Additionally, the velocity in the z-direction calculated by mod-PATH3DU would be incorrect for the upper layer and a particle could only be reliably tracked within layers below.

The following is a list of the MODFLOW boundary package file types that are supported by mod-PATH3DU.

FTYPE (and IFACE key)
Package
WEL Well package
DRN Drain package
RIV River package
EVT Evapotranspiration package
GHB General-head package
RCH Recharge package
FHB Flow and head package
STR Stream package
CHD Constant-head package
SFR Stream-flow routing package
MNW2 Multi-node well package (v.2)
MNW1 Multi-node well package (v.1)
CLN Connected-linear network

In MODPATH, IFACE ranges from 0 to 6, where 1 through 4 represent the different side faces of a cell, 5 and 6 denote the bottom and top faces, respectively, and 0 is used to specify an internal sink/source (i.e. flow is not assigned to a face). mod-PATH3DU supports a limited number of these original IFACE values; notably boundary package flows cannot be assigned to side faces. A list of supported IFACE values is given in the table below. The default IFACE for packages not assigned is 7.

Table 2 - Supported IFACE values and their meaning
IFACE Meaning
0 Internal sink/source; flow for which is included in the \(Q_w\) term in Equation .
2 Implicit side-face; flow is not explicitly assigned to a cell face, rather mod-PATH3DU distributes the flow to any side faces with zero flow for inclusion in the first term on the RHS of Equation .
5 Same as MODPATH. Flow through the bottom face of a cell; which is included in the \(q_v\) term in Equation and in the calculation of velocity in the z-direction.
6 Same as MODPATH. Flow through the top face of a cell; which is included in the \(q_v\) term in Equation and in the calculation of velocity in the z-direction.
7 Internal sink/source. Flow is assigned to the top face (i.e. included in the \(q_v\) term in Equation ), but a particle is terminated as soon as it enters the cell.

IFACE can be specified one of two ways: first via MODFLOW/MODFLOW-USG using the AUX IFACE option, and second using the IFACE setting in the primary mod-PATH3DU input file. Note, if IFACE is specified both ways for a cell, the IFACE specified via MODFLOW is ignored. Review the Input Instructions to learn how to specify IFACE in the primary mod-PATH3DU input file. Generally, there is an IFACE section to this file in which the MODFLOW file type (i.e. CHD, WEL, DRN, RCH or CLN) and the desired IFACE integer are listed. Below is a sample of a mod-PATH3DU primary input file that illustrates how to assign IFACE. To learn more about MODFLOW file types review the NAME file section of the respective manual for each version, some of which are online, for example the Online Guide to MODFLOW-2000 or the Online Guide to MODFLOW-2005.

{
	
"FLOW_MODEL_TYPE"
: {
"MODFLOW"
: { . . .
"IFACE"
: [ {
"WEL"
:
0
}, {
"CHD"
:
2
}, {
"RCH"
:
6
} ], . . . } }, . . . }

In the above example, the MODFLOW model includes the well (WEL), constant head (CHD), and recharge (RCH) packages. The WEL package is assigned an IFACE of 0 to include the local analytic well correction which terminates a particle some user defined distance from the well (see the Input Instructions for details on how to set this radius). The CHD package is assigned an IFACE of 2 so that any unaccounted horizontal flow (the flow in or out of the CHD) is automatically attributed to a side face (for example, implicit no-flow, that is any side faces on the edge of the model domain). The RCH package is assigned an IFACE of 6 to assign any flow to the top face of the cell to ensure vertical velocity is correctly calculated.

An excellent discussion of the impact and importance of IFACE on MODPATH results is given by Abrams et al. (2013). Generally, incorrect assignment (or no specification) of IFACE for a boundary condition will manifest as particle tracks that do not make sense in that they abruptly change direction or oscillate between cells or a cell-corner. When a particle is oscillating it manifests as a long run time, large output file size, or termination due to stagnation. If this is experienced, check the IFACE settings and check the flow balance for the offending cell(s). Also, if AUX IFACE is entered as an option in a boundary package file (through MODFLOW) ensure it is supplied correctly, that is, if AUX IFACE is specified in a MODFLOW file but that variable is not present in the cell listing section of that file, mod-PATH3DU may not perform as expected.

Table 3 - Example IFACE settings
IFACE Key IFACE Value Boundary Purpose Capture (depending on tracking direction)
WEL 0 Extraction or injection well, implicit analytic correction for weak/strong sink. Particle is tracked in boundary cell, captured if within user-specified distance from cell center.
DRN 6 Gaining stream. Particle is tracked in boundary cell, captured when it reaches top face or water table.
7 Horizontal drain. Particle captured when it reaches boundary cell.
0 Level controlled well. Particle is tracked in boundary cell, captured if within user-specified distance from cell center.
RIV 6 Gaining or losing river/stream. Particle is tracked in boundary cell, captured when it reaches top face or water table.
EVT 6 Evaporation. Particle is tracked in boundary cell, captured when it reaches the water table.
GHB 2 Flow laterally into/out of model domain through side face touching implicit or explicit no-flow. Particle is tracked in boundary cell until captured at face.
7 Flow into/out of domain. Particle captured when it reaches boundary cell.
RCH 6 Recharge. Particle is tracked in boundary cell, captured when it reaches water table.
FHB 2 Flow laterally into/out of model domain through side face touching implicit or explicit no-flow. Particle is tracked in boundary cell until captured at face.
7 Flow into/out of model domain. Particle captured when it reaches boundary cell.
STR 6 Gaining or losing river/stream. Particle is tracked in boundary cell, captured when it reaches top face or water table.
CHD 2 Flow laterally into/out of model domain through side face touching implicit or explicit no-flow. Particle is tracked in boundary cell until captured at face.
7 Flow into/out of model domain. Particle captured when it reaches boundary cell.
SFR 6 Gaining or losing river/stream. Particle is tracked in boundary cell, captured when it reaches top face or water table.
MNW2 0 Extraction or injection well, implicit analytic correction for weak/strong sink. Particle is tracked in boundary cell, captured if within user-specified distance from cell center.
MNW1 0 Extraction or injection well, implicit analytic correction for weak/strong sink. Particle is tracked in boundary cell, captured if within user-specified distance from cell center.
CLN - mod-PATH3DU cannot distinguish between different BCs applied to a CLN network. If more than one BC type is applied, use IFACE of 7. If using 7 and a well is being simulated with the CLN, the analytic well correction is not applied, and the well is assumed to be strong and any particles entering the boundary cell are captured.
0 Extraction or injection well, analytic correction. Particle is tracked in boundary cell, captured if within user-specified distance from cell center.
7 Flow into/out of model domain. Particle captured when it reaches boundary cell.

Note - the mod-PATH3DU log file summarizes the boundary packages listed in the MODFLOW cell-by-cell flow file. To just get this list before tracking, remove the SIMULATIONS section from the primary input file (or spell it wrong, i.e. SDIMULATIONS, but not SSIMULATIONS or SIMULATIONSS) and execute. Unless IFACE is being specified via MODFLOW, ensure IFACE is specified in the mod-PATH3DU primary input file for each boundary package listed in the log file summary (do not include STORAGE or the face flows, i.e. FLOW JA/FRONT/SIDE/BOTTOM). Also, if a boundary package is not listed here and an IFACE other than 7 is desired, check your MODFLOW input files and ensure they write to a cell-by-cell flow file.

Weak Sinks and Sources

A limitation of the Pollock method is the treatment of weak sinks and sources. For a weak sink, some of the water flowing to the cell containing it discharges to the sink and some passes through the cell (Pollock, 1994). The Waterloo method does not differentiate between weak or strong sinks - sinks are explicitly represented by the well analytic element term in Equation (third term on the RHS), and accurately tracked to or around.

Backward Tracking and Time Concepts

mod-PATH3DU tracks in the simulation time defined by the flow model. That is, all times input to and output by mod-PATH3DU correspond to the simulation time defined in the flow model. For example, if a MODFLOW or MEUK model has a single stress period or event 1000 days long and a particle's starting time is specified as day 500, mod-PATH3DU will track forward from day 500 to day 1000, and track backwards from day 500 to day 0. It is noted here this is different than MODPATH - which tracks according to a "tracking time" which is 0 at a specified "reference time" in the MODFLOW "simulation time". Technically, mod-PATH3DU does something similar, but this accounting is handled internally and hidden from the user.

To track backwards, mod-PATH3DU uses negative tracking time steps. The user supplies the initial tracking step and mod-PATH3DU internally multiplies this number by -1 when tracking backwards. This approach, again, differs from MODPATH, which reverses the sign of all the velocity components.

Distorted Vertical Discretization & the Water Table

Distorted vertical discretization is handled in mod-PATH3DU in a manner similar to MODPATH. In each cell, elevation, \(z\), is transformed to a local coordinate system, \(z_L\), that ranges from 0 at the bottom of a cell to 1 at the top or water table elevation. By scaling the velocity in the vertical direction by the saturated thickness of the cell, a particle can be tracked in this transformed space provided \(z_L\) is updated when a particle changes layers. The "Non-Rectangular Vertical Discretization" and "Water Table Layers" sections of the MODPATH manual (Pollock, 1994) provide excellent discussions of tracking in this transformed space, as does Zheng (1994).

Quasi Three-dimensional Representation of Confining Layers

Implicit confining layers are not supported now. The program is setup to handle these layers, but at the time of this writing the option was not tested. If there is interest in this option, please contact the corresponding author.

Connected Linear Networks (CLN)

Generally, mod-PATH3DU best supports the CLN package if it is used to simulate a single type of feature, well or drain, for example. Distinguishing between different types of CLN from the GWF cell-by-cell file alone is not possible. So, if the CLN represents a single type, use the appropriate IFACE, if it represents more than one type, it is recommended that an IFACE of 7 is used - note, that if one of the types is a well, an IFACE of 7 assumes this well is a strong sink, and particles are terminated once they enter a cell containing a CLN.

Multithreading (Parallelization) and Reproducibility

With the release of v.2.0.0, mod-PATH3DU is multi-threaded. It is parallelized in two places: when processing the MODFLOW composite budget file and when tracking particles. The number of threads to use for each can be supplied in the primary input file.

When using more than one tracking thread and the dispersion option, reproducibility is not guaranteed because it cannot be guaranteed which thread will be assigned a particle for a time step and each thread is drawing from its own random number generator. If reproducibility is desired in the dispersion calculations, use a single tracking thread.

Limitations

The Waterloo method, because it is semi-analytical, has similar limitations to the Pollock Method. For a detailed description of the Waterloo Method, including limitations, see Ramadhan (2015).

The applicability of the SSP&A method to tracking problems is limited by:

  1. the assumptions in the underlying tracking scheme,
  2. the interpolation method it employs, and
  3. limitations in the groundwater flow model, including discretization and boundary effects.

The accuracy of the universal kriging interpolation is dictated by (a) grid discretization, (b) severe heterogeneities, and (c) proximity to certain boundaries. The refinement obtained by adding more cells spaced closer together will provide the approach with better information to calculate velocity. Nevertheless, the use of universal kriging can overcome many discretization and boundary effects provided an appropriate trend term is used. However, because the coefficients of these drift terms are determined through a regression, they can be made zero or their sign reversed (although this is unlikely). For example, when pumping in the presence of a strong regional gradient, the pumping well signal (its contribution to the hydraulic-head in a cell) may be overwhelmed by the signal from the regional gradient and thus made zero in the regression. Trend terms to overcome hydraulic conductivity dichotomies between two cells and to account for additional boundaries, including no-flow and streams/rivers, are possible, but are not yet available in this version of mod-PATH3DU.

The accuracy of the Runge Kutta scheme depends on the tracking stepsize \(\Delta t\). If \(\Delta t\) is too large, the calculation of the velocity components may be inaccurate and the particle tracking pathline may divert from the actual flow path. However, mod-PATH3DU mitigates this error by implementing the adaptive stepsize control procedure used in PATH3D called "step doubling". For more information on the Runge-Kutta scheme please see the RK4 section in this manual, Zheng (1989;1992;1994) and Zheng and Bennett (2002).

Because the Pollock method is used to calculate velocity in the z-direction, vertical sub-discretization is not supported. Further, the cell structure in each layer must be identical, as in . We have an approach to overcome this, but have not implemented it as we have not required it in our own practice. If you require it and are interested in funding this development, please contact the corresponding author.

Vertical Grid Discretization Limitations

Input Instructions

mod-PATH3DU v.2.0.0 is a rewrite. It is not backward compatible with v.1.x.x and utilities are not provided to convert v.1.x.x inputs into the new format. There are four mod-PATH3DU specific input files:

  1. Primary file,
  2. Particle starting location shapefile,
  3. Per-cell property file, and
  4. Grid specification file (GSF)

The following sections detail these files.

JSON

The format of the primary input file for mod-PATH3DU and its support programs follow the structure and rules of the lightweight data-interchange JSON format. This format has few and simple rules and lends itself to being a flexible and dynamic human readable file.

Example JSON file. Each color represents a unique key:value pair, or object. Purple indicates the primary object (for which a key is not required). Pay particular attention to the use of commas and brackets in this example.

JSON is an organized sequence of key:value pairs. A key must be a string (enclosed in quotation marks), for example "OBJECT_01", and a value can be one of three types,

  1. an integer, real, or string (denoted by double-quotes) value,
  2. a vector of values, or
  3. an object (listing of additional key:value pairs).

Vectors are denoted by square-braces, [ ], objects by curly-braces, { }, and single values by the absence of braces. Multiple objects can be listed for a value, each being separated by a comma.

In the following sections, JSON formatted input instructions are provided as complete examples, context is provided when hovering over a keyword (parameter) or option. Most keywords are optional - mod-PATH3DU will use a default value if a keyword is not specified or terminate in error if it is required. Because the JSON formatted input files are human readable, it is hoped these complete examples make clear the file requirements and flexibility in using the JSON format. In the examples, pay attention to the brackets and commas - this is where most mistakes with the file are made. JSON is supported by a number of open-source libraries that facilitate the reading and writing of files - mod-PATH3DU incorporates JSONCPP for this purpose. Additionally, there are tools available that allow a user to debug the JSON formatting of a file, some of which are online, e.g. JSONEditorOnline.

When checking the JSON object library prepared by JSONCPP for keywords, mod-PATH3DU looks for the word or part of the word, it does not compare two strings. Thus, if you wish to "comment" out a particular setting or option, spell it wrong by putting a dummy character inside the word. For example, mod-PATH3DU will recognize SsIMULATION, dSIMULATION, SIMULATIOnN, and SIMULATIONd as SIMULATION, however it will not recognize SdIMULATION or SIMULATIOdN.

Primary File

The primary input file has two sections - the first indicates the type of model, MODFLOW or MEUK, and the second lists one or more simulations, PATHLINE or ENDPOINT, to execute. Each of these two sections, FLOW_MODEL_TYPE and SIMULATIONS, have required and optional keywords - in the example input files below (choose which to display using the dropdown below), pay particular attention to the first example - in the ENDPOINT simulation section, all of the settings (keywords) and options are listed. These keywords/options are also applicable to the PATHLINE simulation type, but are not listed there to highlight just the required keywords. The second example, is for a MODFLOW model again, but includes only the required keywords. The third is a complete MEUK example. To get started, copy the example to a text editor or use the link below to download the file. Edit as required in your text editor of choice.

Select an example input file to display:
{
	
"FLOW_MODEL_TYPE"
: {
"MODFLOW"
: {
"NAME_FILE"
: "MODFLOW.nam",
"GSF_FILE"
: {
"TYPE"
: "GSF_V.1.1.0",
"FILE_NAME"
: "MP3DU.gsf" },
"OUTPUT_PRECISION"
: "DOUBLE",
"IFACE"
: [ {
"WEL"
:
0
}, {
"CHD"
:
2
}, {
"RCH"
:
6
} ],
"THREAD_COUNT"
: 10 } },
"SIMULATIONS"
: [ {
"ENDPOINT"
: {
"NAME"
: "NAME_OF_EPT_SIMULATION",
"DIRECTION"
: "FORWARD",
"THREAD_COUNT"
: 4,
"INITIAL_STEPSIZE"
: 0.1,
"MAX_STEPSIZE"
: 1.0e6,
"STAGNATION_DT"
: 1.0e-15,
"EULER_DT"
: 1.0e-4,
"ADAPTIVE_STEP_ERROR"
: 1.0e-6,
"SIMULATION_END_TIME"
: 3650.0,
"CAPTURE_RADIUS"
: 10.0,
"OPTIONS"
: [
"DISPERSION"
,
"RETARDATION"
,
"TRACK_TO_TERMINATION"
, ],
"PARTICLE_START_LOCATIONS"
: {
"REPEAT"
: 1,
"REPEAT_DT"
: 1.0e30,
"SHAPEFILE"
: {
"FILE_NAME"
: "EPT_PartStart.shp",
"CELLID_ATTR"
: "P3D_CellID",
"TIME_ATTR"
: "TREL",
"ZLOC_ATTR"
: "ZLOC",
"ADDTL_ATTR"
: [ "ATTR_01", "ATTR_02" ] } } } }, {
"PATHLINE"
: {
"NAME"
: "NAME_OF_PTL_SIMULATION",
"PARTICLE_START_LOCATIONS"
: {
"SHAPEFILE"
: {
"FILE_NAME"
: "PTL_PartStart.shp",
"CELLID_ATTR"
: "P3D_CellID",
"TIME_ATTR"
: "TREL",
"ZLOC_ATTR"
: "ZLOC" } } } } ] }
Click here to download the Complete MODFLOW example file.

Particle Starting Locations Shapefile

Particle starting locations are provided in a shapefile. The shapefile must be a POINT shapefile and contain attributes indicating:

  1. cell id - the id number of the cell containing the particle,
  2. local z-elevation - the relative elevation of the particle, and
  3. the release time.

Each mod-PATH3DU "SIMULATION" must include a "PARTICLE_START_LOCATIONS" section. This section requires a "SHAPEFILE" keyword through which the user supplies the filename and the attribute headers corresponding to the required "CELLID_ATTR", "TIME_ATTR", and "ZLOC_ATTR" keywords. See examples in the Primary File section.

Additional attributes from this shapefile can be passed through to the mod-PATH3DU binary output file. For example, if particles are released from a particular well or area of a site, an identifier attribute - string (upto 50 characters), integer or real number - can be associated with each particle. The writeP3DOutput.exe program, supplied with mod-PATH3DU, which writes a shapefile of pathlines, for example, will then write these pass-through attributes to the final shapefile output. This feature can be useful to colorcode pathlines for different wells, for example. An example of this feature is provided in the example Primary File. There is an optional keyword "ADDTL_ATTR" in the "SHAPEFILE" subsection of the "PARTICLE STARTING LOCATIONS" section. Any pass-through attribute headers are listed here.

Per-cell Property File (MODFLOW only)

The per-cell property file lists values for properties required by the particle tracking that may vary in each model cell, such as porosity and dispersivity parameters. This file is only required for MODFLOW and must be included in the MODFLOW NAME file (blue text in example below) with the PATH type. MODFLOW-USG v.1.4 and later ignore the PATH file type if it is listed in the NAME file; older versions, and MODFLOW-2000 and MODFLOW-2005, terminate in error when unknown types are listed. As such, a user may have to comment (#) the PATH type in the NAME file when running the flow model.

Example MODFLOW NAME file with PATH type:


LIST  7 dispdemoMNW.lst
BAS6  1 dispdemo.bas
DIS  29 dispdemoWEL.dis
LPF  11 dispdemo.lpf
CHD  26 dispdemoWEL.chd
OC   22 dispdemo.oc
MST  21 dispdemoWEL.mst
ORT  23 dispdemoWEL.ort
WEL  25 dispdemoWEL.wel
PATH
24
mp3du.p3d
DATA 27 retardation.dat DATA(BINARY) 50 dispdemoWEL.cbb DATA(BINARY) 30 dispdemoWEL.hds DATA(BINARY) 31 dispdemoWEL.ddn LMT6 33 dispdemoWEL.lmt

The following is a complete example Per-cell Property file for a model with 2 layers and 5 cells per layer. This file is structured as a MODFLOW input file, which allows properties to be specified on a per-cell basis. It uses MODFLOW routines, UXDINT and UXDREL, to read in the property arrays. Context for the different parts of the file is provided when hovering over a particular parameter or file part below. MODFLOW allows for per-cell arrays to be provided several ways, including an option to specify a constant value for each model cell. The example input file below attempts to illustrate each of these options. However, in preparing this file for mod-PATH3DU, any one or combination of these headers can be used as needed.


# PATH3D input file
CONSTANT
3
(FREE)
-1
VELOCITY METHOD LAYER 1
CONSTANT
3
(FREE)
-1
VELOCITY METHOD LAYER 2
INTERNAL
1.00
(FREE)
-1
POROSITY L1
0.10 0.20 0.30 0.40 0.50
INTERNAL
1.00
(FREE)
-1
POROSITY L2
0.60 0.70 0.80 0.90 1.00
EXTERNAL
27
1.0
(FREE)
-1
RETARDATION L1
EXTERNAL
27
1.0
(FREE)
-1
RETARDATION L2
OPEN/CLOSE
FNAME
1.0
(FREE)
-1
DISPH L1
OPEN/CLOSE
FNAME
1.0
(FREE)
-1
DISPH L2
CONSTANT
1.0
(FREE)
-1
DISPT L1
CONSTANT
1.0
(FREE)
-1
DISPT L2
CONSTANT
0.0
(FREE)
-1
DISPV L1
CONSTANT
0.0
(FREE)
-1
DISPV L2

Grid Specification File (GSF)

The grid specification file provides x and y-coordinates for the vertices of each cell in a MODFLOW-USG or MEUK grid. The format of this file required by mod-PATH3DU is a modified form of the original format written by most GUIs and required for PEST, for example. The mod-PATH3DU-specific GSF format is less general than the original format because the requirements for mod-PATH3DU are more restrictive - the restrictions are (cell and vertex numbers correspond to ):

  1. Duplicate vertices are not listed. For example, between two identical square cells, cells 1 and 2, there are two shared vertices, 2 and 5. These shared vertices must be listed only once - and referenced by both cells.
  2. The z-coordinate is not needed when listing vertex coordinates. mod-PATH3DU gets elevations from the discretization and head-save files.
  3. For each cell, the vertices must be listed clockwise. It does not matter which vertex is listed first for a cell, but the subsequent vertices must be listed clockwise from that one. For example, for cell 2 in the figure below, its vertices can be listed as 4,5,2,3 or 2,3,4,5, etc.
  4. For quad-based grids - if a cell has more than one connection on a side, all of the shared vertices must be listed for that cell. For example, for cell 1 in the figure below, its vertices are 1,2,5,7,10 - vertex 7 must be included.
  5. For each vertex pair making a cell "face", the connected cell that shares that face must be listed. For example, for cell 1, if the vertices are listed as 5,7,10,1,2, the cell that shares the face given by vertices 5 and 7 is 4, and for the face given by 7 and 10, the shared cell is 3. For the face given by 10 and 1 a -999 is required to indicate there is not a connected cell; the same for the face defined by 1 and 2. The face between vertices 2 and 5 share cell 2.
Cells and vertices referenced in GSF specifications listed above.

A utility is provided with mod-PATH3DU that will convert a GSF file in the original format into the mod-PATH3DU specific one. This program, writeP3DGSF.exe, is described in detail below, but in summary, will:

  1. Check for and exclude duplicate vertices - its tolerance is 1.e-4; any vertices within 1.e-4 units of each other are considered identical.
  2. Remove any reference to z-coordinates.
  3. List vertices for each cell clockwise.
  4. Ensure all vertices are listed for any cells with a side connected to more than one other cell. For example, in above, if the GSF in the original format only lists vertices 1,2,5,10 for cell 1, writeP3DGSF.exe will automatically add vertex 7 to those listed for this cell.
  5. Determine the connected cell that shares each face of a cell.
  6. Write a new GSF file in the required format.

The following is a complete example GSF file corresponding to the grid shown in . This file has two main sections, first a listing of the unique vertices, followed by a listing of the vertices and neighboring cells for each model cell. Context for the different parts of the file is provided when hovering over a particular parameter or file part below.


mod-PATH3DUv110
10
1
0.0
7.5
2
5.0
7.5
3
10.0
7.5
4
10.0
2.5
5
5.0
2.5
6
5.0
0.0
7
2.5
2.5
8
2.5
0.0
9
0.0
0.0
10
0.0
2.5
1
2.5
5.0
5
5 7 10 1 2
4 3 -999 -999 2
2
7.5
5.0
4
4 5 2 3
-999 1 -999 -999
3
1.25
1.25
4
8 9 10 7
-999 -999 1 4
4
3.75
1.25
4
6 8 7 5
-999 3 1 -999

WriteP3DGSF.exe

Regardless of whether a model is structured or unstructured, mod-PATH3DU requires a GSF file in a program specific format. writeP3DGSF.exe is provided with mod-PATH3DU to translate grids defined in other formats into the required GSF format. Primarily, writeP3DGSF.exe is used to convert a GSF written by one of the primary GUIs in the original format, herein refered to as v.1.0.0, into this new format. Users of mod-PATH3DU v.1.1.x will be familiar with this program. Its role has not changed - but its functionality has. In addition to converting a GSF from one format to another, writeP3DGSF.exe will also write a GSF file for a structured grid defined by a MODFLOW discretization file with an option to include an offset and rotation to write the GSF in real-world coordinates. Recall, mod-PATH3DU tracks on the grid as it is defined in the GSF. There is no need to track in "model" coordinates.

The program is executed from the command line as follows:

 
writeP3DGSF.exe
NAME.json
colorcode

The program has a single input file - using the same formatting as the primary input file for mod-PATH3DU. Below, complete example input files are provided. Each file indicates all of the required and optional keywords. Context for each keyword is provided when hovering over it. The input file requires a FLOW_MODEL_TYPE, which can be MODFLOW or MEUK. If the type is MODFLOW, the type of conversion is specified through the GSF_FILE section.

Select an example input file to display:
{
  
"FLOW_MODEL_TYPE"
: {
"MODFLOW"
: {
"NAME_FILE"
: "MODFLOW.nam",
"GSF_FILE"
: {
"TYPE"
: "GSF_V.1.0.0",
"FILE_NAME"
: "StandardFormat.gsf" } } },
"OUTPUT_FILENAME"
: "mp3du.gsf" }
Click here to download this example file.

In the above example, a GSF in the standard format, indicated in the GSF_FILE section as GSF_V.1.0.0, will be converted to the mod-PATH3DU specific format. The MODFLOW name file is required because writeP3DGSF.exe uses the cell connections in the JA array to verify shared vertices.

Output

Binary Format

mod-PATH3DU writes a single binary file regardless of simulation type or options. Writing to the output file is the sole bottlekneck when tracking particles in parallel - each thread must wait its turn for access to the output file. The longer a thread takes to write its particle position to a file, the longer other threads sit idle waiting for their turn. A simple binary file has little overhead, taking less time to write. The mod-PATH3DU suite includes a program for writing other file formats, especially shapefiles, from this file. This program, writeP3DOutput.exe, is described in the next section.

The following table details the format of the binary file. The format is likely to change over time - if processing this file directly - be sure to note the version and order of attributes. The first part of the file is a header section that gives the version followed by an ordered listing of the attributes for each particle position in the file. This listing includes the number of attributes, followed by the name, size in bytes, and an integer flag indicating the type of each attribute. The second part of the file lists these attributes for each calculated particle position. Any pass-through attributes from the particle starting locations file are included and listed here as well.


Variable       Type      Size(Bytes)   Description
------------   -------   -----------   ---------------------------------------------
VERSION CHAR 16 File version
NATTR INTEGER 4 Number of attributes for each particle point
|- ATTR_NAME CHAR 10 Name of attribute |- ATTR_SIZE INTEGER 4 Size (bytes) of the attribute |- ATTR_TYPE INTEGER 4 Attribute type; 1:Integer, 2:Double, 3:Char
ATTRIBUTE INTEGER 4 | DOUBLE 8 |- Each attribute is one of these types. CHAR 51 |

In the case of an endpoint simulation, in which the number of pathline positions is 1 per particle, the total number of attributes listed in the output file is NATTR x NPARTICLE x 1 POSITION / PARTICLE. Table 4 lists the attributes currently written to the binary file.

Table 4 - Output Attributes
Name Type Size (bytes)> Description
PID INTEGER 4 Particle id
CELLID INTEGER 4 Cell containing particle
PTIME DOUBLE 8 Time
PX DOUBLE 8 X-coordinate
PY DOUBLE 8 Y-coordinate
PZGLO DOUBLE 8 Z-coordinate (global)
PZLOC DOUBLE 8 Z-coordinate (local)
PVX_ADV DOUBLE 8 Advective velocity in the x-direction
PVY_ADV DOUBLE 8 Advective velocity in the y-direction
PVZ_ADV DOUBLE 8 Advective velocity in the z-direction
PLAYER INTEGER 4 Layer of cell containing particle
PTSTART DOUBLE 8 Release time of particle
PTERM CHAR 51 Termination reason

WriteP3DOutput.exe

WriteP3DOutput.exe is a program to write more commonly used output file types from the binary file produced by mod-PATH3DU. It can write the following

  • ASCII Table - A listing of each of the calculated particle positions, sorted by particle id (PID), in a common text file format.
  • DBF Table - A listing of each of the calculated particle positions, unsorted, in the DBF file format. This file can be posted in GIS, or viewed as a spreadsheet.
  • Shapefile - A shapefile, in one of the following formats:
    • Pathline whole - a single POLYLINE for each particle.
    • Pathline parts - each part of a particle's path is a POLYLINE.
    • Points in time - POINTS along a pathline, interpolated in time.
    • Endpoint - The vertex of every calculated particle position as a POINT.
  • Summary - Nothing is written, except a summary of the mod-PATH3DU binary output file to the console.

The program is executed from the command line as follows:

 
writeP3DOutput.exe
NAME.json
colorcode

The program has a single input file - using the same formatting as the primary input file for mod-PATH3DU. Below, is a complete example input file; it indicates all of the required and optional keywords. Context for each keyword is provided when hovering over it. It lists all of the different output options. Users can select one or more output type to write as needed.

{
  
"MP3DU_BIN"
: "NAME_OF_PTL_SIMULATION_PATHLINE.bin",
"OUTPUTS"
: [ {
"SUMMARY"
: { } }, {
"ASCII_TABLE"
: {
"FILE_NAME"
: "NAME_OF_OUTPUT.dat" } }, {
"DBF_TABLE"
: {
"FILE_NAME"
: "NAME_OF_OUTPUT.dbf" } }, {
"PATHLINE_WHOLE"
: {
"FILE_NAME"
: "NAME_01_OF_OUTPUT.shp",
"TIME_CONVERSION"
: 2.7397E-03,
"MAX_TIME"
: 1.0,
"MIN_TIME"
: 3.0 } }, {
"PATHLINE_PARTS"
: {
"FILE_NAME"
: "NAME_02_OF_OUTPUT.shp",
"TIME_CONVERSION"
: 1.0,
"MAX_TIME"
: 365.0,
"MIN_TIME"
: 1095.0 } }, {
"POINTS_IN_TIME"
: {
"FILE_NAME"
: "NAME_03_OF_OUTPUT.shp",
"TIME_CONVERSION"
: 2.7397E-03,
"TIME_DT"
: 1.0 } }, {
"ENDPOINT"
: {
"FILE_NAME"
: "NAME_04_OF_OUTPUT.shp" } } ] }
Click here to download this example file.

Examples

Select an example:

Example 1a

Tracking Near a "Weak" Pumping Well (Structured Grid)

The purpose of this example is to highlight the strength of the Waterloo and SSP&A methods in tracking near a pumping well - the analytical correction term implicitly handles weak and strong sinks. The groundwater system is confined (50 feet thick), steady-state and two-dimensional, in plan-view. The ambient groundwater flow field is uniform throughout the whole domain and flow is from left to right. The fully penetrating well is pumped at a constant rate of 50 feet3/day. Hydraulic conductivity is 10 feet/day and porosity is 0.3. The numerical domain consists of 21 rows, each 5 feet wide, and 87 columns with varying widths ranging from 5 to 20 feet, for a total model width of 500 feet. The conceptual model is shown in , below.

Conceptual model - structured grid

Particle paths, near the pumping well, calculated by MODPATH6 and both the Waterloo and SSP&A methods of mod-PATH3DU are shown in . The Waterloo method tracks the same as the Pollock method, but also calculates pathlines in the cell containing the weak well – unlike the Pollock method. For both the Waterloo and SSPA methods, pathlines within the analytically calculated capture zone (using Jacob, 1949; p.344.), are captured by the well, while those outside, flow through the cell and continue tracking to their proper terminus. The head contours presented are an expression of what mod-PATH3DU tracks on in the vicinity of a well. They were calculated by interpolating from the MODFLOW calculated heads to a very fine grid using the universal kriging algorithm available with mod-PATH3DU. The shape visible around the well is determined by the linear-log drift term applied in this case.

Particle paths calculated using MODPATH6 and mod-PATH3DU near a "weak" sink

Because both the Waterloo and SSP&A methods include an analytic solution when tracking near a well, they can mitigate some coarse grid discretization issues. Further, they both can simulate off-center pumping wells. An example of how mod-PATH3DU can be used to calculate capture zones for off-center pumping wells (with one or more in a cell) using the Area Based Redistribution (ABRD) approach of Pinales et al. (2003;2005) is presented in Muffels et al. 2011.

Click here to download this example.

References

Abrams, D., H. Haitjema, and L. Kauffman, 2012, On Modeling Weak Sinks in MODPATH: Ground Water 51, no. 4: 597-602, doi 10.1111/j.1745-6584.2012.00995.x.
Burnett, R.D., and E.O. Frind, 1987, Simulation of contaminant transport in three dimensions, 2, Dimensionality effects. Water Resources Research, Vol. 23, No. 4, 695.
Cressie, N., 1993, Statistics for Spatial Data, revised ed., 900pp., Wiley, New York.
Deutsch, C. V. and Journel, A. G., 1998, GSLIB: Geostatistical Software Library and User's Guide, 2nd ed. Oxford University Press, Oxford.
Haitjema, H.M., 1995, Analytic Element Modeling of Groundwater Flow, Academic Press Inc., San Diego, California.
Harbaugh, A.W., 2005, MODFLOW-2005, the U.S. Geological Survey modular ground-water model—The Ground-Water Flow Process: U.S. Geological Survey Techniques and Methods 6–A16, variously paginated.
Jacob, C.E., 1949, Flow of ground water, in Engineering Hydraulics, Proceedings of the 4th Hydraulics Conference, Iowa Institute of Hydraulic Research, H. Rouse (ed.), June 12-15, 1949, John Wiley & Sons, Inc., New York.
Jankovic, I., and Barnes R., 1999a, High-order line elements in modeling two-dimensional groundwater flow, Journal of Hydrology 226, pages 224-233.
Jankovic, I., and Barnes R., 1999b, Three-dimensional flow through large numbers of spheroidal inhomogeneities, Journal of Hydrology 226, pages 224-233.
Journel, A.G. and Huijbregts, C.J., 1992, Mining Geostatistics. Academic Press, New York.
Karanovic, M., M. Tonkin, and D. Wilson, 2009, KT3D_H20: A Program for Kriging Water-Level Data Using Hydrologic Drift Terms: Ground Water 45, no. 4: 580-586, doi 10.1111/j.1745-6584.2009.00565.x.
Kinzelbach, W., 1990, Simulation of pollutant transport in groundwater with the random walk method. Groundwater Monitoring and Management (Proceedings of the Dresden Symposium, March 1987). IAHS Publ. no. 173, 1990.
LaBolle, E.M., Quastel, J., Fogg, G.E., and Gravner, J., 2000, Diffusion processes in composite porous media and their numerical integration by random walks: Generalized stochastic differential equations with discontinuous coefficients. Water Resources Research, Vol. 36, No. 3, Pages 651-662.
Muffels, C., Stonebridge, G., Tonkin, M.J., and Karanovic, M., 2011, An Unstructured Version of PATH3D, PATH3DU. Presentation at MODFLOW and More 2011, Integrated Hydrologic Modeling, International Ground Water Modeling Center, Colorado School of Mines, Golden, CO, June 6-8, 2011, v. 1, pp. 272-275.
Neville, C.J., 2006, Analytical Solutions for Solute Transport with One-Dimensional Flow: Brick Sources in an Infinite Aquifer, S.S. Papadopulos & Associates, Inc., Waterloo, Ontario.
Panday, S., Langevin, C.D., Niswonger, R.G., Ibaraki, M., and Hughes, J.D., 2013, MODFLOW-USG version 1: An unstructured gird version of MODFLOW for simulating groundwater flow and tightly coupled processes using a control volume finite difference formulation: U.S. Geological Survey Techniques and Methods 6 A45.
Pinales, A., A. Keer, F. Espinosa, L. Manzanares, G. Llerar, and A. Chavez, 2003, An Alternative Approach for Assigning Well Rates in a Block-Centered Finite-Difference Grid. In MODFLOW and More 2003: Understanding through Modeling - Conference Proceedings.
Pinales, A. Chavez, G. Llerar, L. Manzanares, and A. Keer, 2005, An Improved Approach for Assigning Pumping Rates to Heterogeneous Aquifer Models. Ground Water, vol. 43:274-279.
Pollock, D.W., 1988, Semi-analytical computation of pathlines for finite difference models, Ground Water vol. 26, no. 6: 743-750.
Pollock, D.W., 1989, Documentation of a computer program to compute and display pathlines using results from the U.S. Geological Survey modular three-dimensional finite-difference ground-water flow model: U.S. Geological Survey Open-File Report 89-381.
Pollock, D.W., 1994, User's guide for MODPATH/MODPATH-PLOT, version 3: A particle-tracking post-processing package for MODFLOW, the U.S. Geological Survey finite-difference ground-water flow model: U.S. Geological Survey Open-File Report 94-464.
Pollock, D.W., 2012, User's guide for MODPATH version 6 A particle-tracking model for MODFLOW: U.S. Geological Survey Techniques and Methods 6-A41, 58p.
Pollock, D.W., 2015, Extending the MODPATH algorithm to Rectangular Unstructured Grids, Ground Water vol. 54, no. 1: 121–125. doi: 10.1111/gwat.12328.
Prickett, T.A., Naymik, T.G., and Lonnquist, C.G., 1981, A “Random-Walk” Solute Transport Model for Selected Groundwater Quality Evaluations. Illinois State Water Survey, Champaign, Bulletin 65, 1981.
Rhamadhan, M., 2015, A Semi-Analytical Particle Tracking Algorithm for Arbitrary Unstructured Grids. Unpublished MASc. Thesis. University of Waterloo, Waterloo Ontario.
Salamon, P., 2006, On modeling contaminant transport in complex porous media using random walk particle tracking. Phd Thesis. Universidad Politecnica de Valencia.
Skrivan, J. A., and Karlinger, M. R., 1980, Semi-variogram estimation and universal kriging program; U. S. Geological Survey Computer Contribution, 98 p. Tacoma, Washington. (Computer Program K603).
Tonkin, M.J., and Larson, S. P., 2002, Kriging Water Levels with a Regional-linear and Point-logarithmic Drift. Ground Water, vol. 40, no. 2: 185-193.
Zheng, C., 1989, PATH3D, A ground-water path and travel-time simulator, version 3.0 user’s manual, S.S. Papadopulos & Associates, Inc., Bethesda, MD.
Zheng, C., 1992, PATH3D, A ground-water path and travel-time simulator, version 3.2 user’s manual, S.S. Papadopulos & Associates, Inc., Bethesda, MD.
Zheng, C., 1994, Analysis of particle tracking errors associated with spatial discretization, Ground Water, vol. 32 no. 5: 821-828.
Zheng, C., and Bennett, G., 2002. Applied Contaminant Transport Modeling, 2nd ed., John Wiley & Sons, New York.