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.
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
1801 Rockville Pike, Suite 220
Rockville, MD 20852
301-718-8900
cmuffels@sspa.com
www.sspa.com
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.
The authors acknowledge the following interns for their contributions to the development of mod-PATH3DU.
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.
mp3du.exe NAME.json colorcode
mod-PATH3DU is a particle tracking code for calculating the three-dimensional flow pathlines and travel time of solute particles. It supports the latest MODFLOW releases: MODFLOW-USG (Panday et al., 2013) and MODFLOW 6 (Langevin et al., 2017). Further, it supports tracking on a map of water levels interpolated from monitoring well data using MEUK (Tonkin et al., 2016). MODFLOW-USG and MODFLOW6 are unstructured-grid versions 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.
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.
The MODFLOW 6 binary output files - head-save and cell-by-cell flow - are largely unchanged from other MODFLOW releases. However, in addition to these files, MODFLOW 6 also writes a Binary Grid File (GRB) that lists the requisite cell vertex coordinates. The writeP3DGSF.exe program was updated to process this file and write the required GSF file for mod-PATH3DU. Further, for other MODFLOW releases, the necessary input and output files are supplied to mod-PATH3DU through the NAME file directly, however, these files must be explicitly listed for MODFLOW 6 as detailed in the Input Instructions section.
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.
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
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
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:
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,
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.
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
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,
to yield an effective velocity, \( \mathbf V_{RK4}(\mathbf X_p,t) \), that results in a fourth-order accurate tracking step:
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.
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:
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,
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:
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:
The allowable error is
where
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:
Taking the ratio of () to (),
yields an estimate of the maximum allowable stepsize that satisfies the user-specified stability criterion:
() 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
This estimate of the particle position is fifth-order accurate.
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)
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:
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:
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:
Equations () and () are not identical if the dispersion coefficients are space dependent. Adding a correction term
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)
where \(\mathbf {BZ}\) is the (random) displacement due to dispersion, \(\mathbf B\) is the displacement matrix and is related to dispersion, \(\mathbf D\), by
Dispersion, in three-dimensions, is defined following Burnett and Frind (1987) as
Thus,
And, the velocity correction, (), termed the "drift" vector, \(\mathbf A \left( \mathbf X_p,t \right)\), is
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
where
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" ],
}
.
.
.
}
]
}
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 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):
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):
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:
and
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):
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 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.
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).
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 |
LAK | Lake package |
SGB | Specified-gradient package |
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.
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.
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. | |
LAK | 7 | Flow into/out of domain through LAK cell (inactive in GWF domain). | Particle captured when it reaches boundary cell. |
SGB | 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. |
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.
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.
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 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).
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.
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.
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.
Version 2.X.X is limited to a maximum of approximately 2 million cells. An expected future release will support much more - email the corresponding author if you are interested in trying a beta version of this release.
MODFLOW 6 is supported through the USGS_HFWK model type in the primary input file. Unlike support for other versions of MODFLOW where the NAME file is used to supply input and output file names to mod-PATH3DU, the required input and output files must be explicitly listed for MODFLOW 6. This requirement is detailed in the Input Instructions. Support is available to convert DIS and DISV type grids to a GSF file. While there is support of DISU types, it is limited because mod-PATH3DU tracks on a per-layer basis and the DISU grid type is not. However, mod-PATH3DU automatically assigns a layer to each model cell by looking at each cell's connections and the top and bottom elevation of each - if a connected cell's top elevation is equal to this cell's bottom elevation, the connected cell is assigned to the next layer below.
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:
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.
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:
The following sections detail these files.
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.
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,
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.
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.
{
"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"
}
}
}
}
]
}
{
"FLOW_MODEL_TYPE" : {
"MODFLOW" : {
"NAME_FILE" : "MODFLOW.nam",
"GSF_FILE" : {
"TYPE" : "GSF_V.1.1.0",
"FILE_NAME" : "MP3DU.gsf"
},
"IFACE" : [ { "WEL" : 0 }, { "CHD" : 2 }, { "RCH" : 6 } ]
}
},
"SIMULATIONS" : [
{
"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"
}
}
}
}
]
}
{
"FLOW_MODEL_TYPE" : {
"USGS_HFWK" : {
"GRB_FILE" : "MODFLOW6.grb",
"TDIS_FILE" : "MODFLOW6.tdis",
"PATH_FILE" : "MP3DU.p3d",
"HDS_FILE" : "MODFLOW6.hds",
"CBB_FILE" : "MODFLOW6.cbb",
"GSF_FILE" : {
"TYPE" : "GSF_V.1.1.0",
"FILE_NAME" : "MP3DU.gsf"
},
"IFACE" : [ { "RCHA" : 6 } ]
}
},
"SIMULATIONS" : [
{
"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"
}
}
}
}
]
}
{
"FLOW_MODEL_TYPE" : {
"MEUK" : {
"GRID" : {
"NCOLS" : 100,
"NROWS" : 100
},
"GSF_FILE" : {
"TYPE" : "GSF_V.1.1.0",
"FILE_NAME" : "MEUK.gsf"
},
"HHK" : {
"FILE_NAME" : "MEUK_HHK.dat"
},
"POROSITY" : {
"CONSTANT" : 0.2
},
"DISPERSION" : {
"LONGITUDINAL" : 1.0,
"TRANSVERSE_HORIZONTAL" : 0.1
},
"OPTIONS" : [ "KRIGE_ONLY" ],
"SEARCH_RADIUS" : 1000.0,
"DRIFTS" : [
{
"WELL" : {
"FILE_NAME" : "MEUK_WELLS.csv",
"XCOORDS" : "XCOORDS",
"YCOORDS" : "YCOORDS",
"VAL" : "VAL",
"TERM" : "TERM",
"NAME" : "NAME",
"EVENT" : "EVENT"
},
"LINESINK" : {
"FILE_NAME" : "MEUK_LINESINKS.csv",
"XCOORDS" : "XCOORDS",
"YCOORDS" : "YCOORDS",
"XCOORDS_2" : "XCOORDS_2",
"YCOORDS_2" : "YCOORDS_2",
"VAL" : "VAL",
"TERM" : "TERM",
"NAME" : "NAME",
"EVENT" : "EVENT"
}
}
],
"EVENTS" : [
{
"Event01" : {
"DURATION" : 365.0,
"FILE_NAME" : "Heads_Event01.csv"
},
"Event02" : {
"DURATION" : 365.0,
"FILE_NAME" : "Heads_Event02.csv"
}
}
]
}
},
"SIMULATIONS" : [
{
"ENDPOINT" : {
"NAME" : "NAME_OF_EPT_SIMULATION",
"DIRECTION" : "FORWARD",
"THREAD_COUNT" : 4,
"INITIAL_STEPSIZE" : 0.1,
"MAX_STEPSIZE" : 1.0e6,
"ADAPTIVE_STEP_ERROR" : 1.0e-6,
"SIMULATION_END_TIME" : 3650.0,
"CAPTURE_RADIUS" : 10.0,
"OPTIONS" : [
"DISPERSION",
"RETARDATION",
"TRACK_TO_TERMINATION",
"DOPRI"
],
"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"
}
}
}
}
]
}
Particle starting locations are provided in a shapefile. The shapefile must be a POINT shapefile and contain attributes indicating:
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.
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
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 ):
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:
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
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.
{
"FLOW_MODEL_TYPE" : {
"MODFLOW" : {
"NAME_FILE" : "MODFLOW.nam",
"GSF_FILE" : {
"TYPE" : "GSF_V.1.0.0",
"FILE_NAME" : "StandardFormat.gsf"
}
}
},
"OUTPUT_FILENAME" : "mp3du.gsf"
}
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.
{
"TRANSFORMATION" : {
"XOFF" : 1234.,
"YOFF" : 5678.,
"ROT" : 9.
},
"FLOW_MODEL_TYPE" : {
"MODFLOW" : {
"NAME_FILE" : "MODFLOW200X.nam",
"GSF_FILE" : {
"TYPE" : "STRUCTURED_DIS"
}
}
},
"OUTPUT_FILENAME" : "mp3du.gsf"
}
In the above example, a structured grid, defined in a MODFLOW DIS file, will be converted to the mod-PATH3DU specific GSF format. The GSF_FILE section indicates a "STRUCTURED_DIS" is to be converted, where the DIS file is identified in the MODFLOW name file.
{
"TRANSFORMATION" : {
"XOFF" : 1234.,
"YOFF" : 5678.,
"ROT" : 9.
},
"FLOW_MODEL_TYPE" : {
"USGS_HFWK" : {
"GRB_FILE" : "MODFLOW6.grb",
"GSF_FILE" : {
"TYPE" : "HFWK_GRB_V.1.0.0"
}
}
},
"OUTPUT_FILENAME" : "mp3du.gsf"
}
In the above example, a structured grid (DIS type) or unstructured grid (DISV or DISU), defined in a MODFLOW6 GRB file, will be converted to the mod-PATH3DU specific GSF format. The GSF_FILE section indicates a "HFWK_GRB_V.1.0.0" is to be converted, where the GRB file is identified in the USGS_HFWK section.
{
"FLOW_MODEL_TYPE" : {
"MEUK" : {
"GRID" : {
"NCOLS" : 250,
"NROWS" : 250,
"XLLCORNER" : 0.,
"YLLCORNER" : 0.,
"CELLSIZE" : 100.
}
}
},
"OUTPUT_FILENAME" : "mp3du.gsf"
}
In the above example, a structured MEUK grid will be converted to the mod-PATH3DU specific GSF format. There is only one conversion/grid definition option for MEUK.
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.
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 is a program to write more commonly used output file types from the binary file produced by mod-PATH3DU. It can write the following
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"
} },
{ "V100_PTL" : {
"FILE_NAME" : "NAME_OF_OUTPUT.ptl"
} },
{ "V100_EPT" : {
"FILE_NAME" : "NAME_OF_OUTPUT.ept"
} }
]
}
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.
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.
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.
The modeling exercise in Example 1a is repeated using an unstructured grid defined by Voronoi cells, .
The mod-PATH3DU calculated particle paths for this model are shown on .
This example is adapted from the documentation of PATH3D ver. 4.6. The example was developed to test the calculation of groundwater pathlines in a heterogeneous cross-section. The groundwater model comprises three strata; the hydraulic conductivity of each stratum is shown in . The model is discretized into 20 layers. The cross-section is bordered by no-flow boundaries along the left, right and bottom edges. The top boundary is assigned specified-heads to represent a linear water table declining from 41 feet at x = 100 feet to 40 feet at x = 0 feet. Porosity is 0.2.
A particle is tracked backwards from x = 1 foot. The particle paths calculated using PATH3D v.4.6 and Waterloo method in mod-PATH3DU for this particle are shown in . The total travel time calculated by each approach are nearly identical; PATH3D calculated 637 days, while mod-PATH3DU calculates 636 days. This example primarily serves to demonstrate the implementation of the Pollock method in the z-direction when used in conjunction with the Waterloo method is correct and that the method can be applied to cross-section models.
Example 3 is taken from Pollock (2015). It is a two-dimensional steady-state model with a quad-based unstructured grid. The system is confined with a constant hydraulic conductivity (100 feet/day) and porosity (0.25). No-flow boundaries surround the model on all four sides. Flow circulation is induced with constant-head cells in the center of the model domain, from 15 feet to 2 feet. There are 352 cells in the unstructured grid. mod-PATH3DU was used to calculate flow paths for the unstructured grid model and these were compared with the results in Pollock (2015) - . The flow paths calculated by mod-PATH3DU are comparable to those calculated by the algorithm described by Pollock (2015).
The purpose of this example is to demonstrate the Random-Walk dispersion module. The conceptual flow model is illustrated in . The grid spacing in both the x and y-directions is 10 m. The grid contains 60 columns extending in the x-direction and 10 rows extending in the y-direction. The system is confined, with two constant head boundaries to produce a uniform groundwater flow field with an advective velocity of 0.25 m/day. Ten-thousand particles are started in the source area, centered about coordinate (50,50), and tracked forwards from their release points for 1000 days. Longitudinal and transverse dispersivities are 5 m and 0.05 m, respectively.
After some time, for the case of an instantaneous release at a point in a uniform flow field, the resulting concentration follows Gaussian (or Normal) distributions in the x and y-directions. The center-of-mass of the plume, or the mean position in x and y, representing the distance traveled based on the average linear groundwater velocity, after 1000 days is 250 m from the source, at x = 300 m, and y = 50 m. The mean x calculated from the MP3DU output is 299.6 m, and the mean y is 50.0 m. The longitudinal standard deviation is \(\sigma=\sqrt{2\alpha_lvt}=\)50 m, and transverse it is \(\sigma=\sqrt{2\alpha_t^hvt}=\)5 m. The corresponding standard deviations calculated from the MP3DU output are 50.4, and 5.0, respectively. The frequency distributions in both x and y are shown in . The distributions approximate closely Gaussian distributions.
Concentration can be calculated by determining the mass at each particle and evaluating the particle frequency on a regular grid. For this calculation a total of 300,000 particles,\(N_{particles}\), were released, and each was assigned a mass
of 0.004 μg, where \(C_0\) is the source concentration, 1000.0 μg/m3, and the pore volume of the source, \( \theta V_{source}\), is \(0.3 \cdot 4=1.2 m^2\) (particles are released within a 2x2 m area). After tracking, the number of particles within each cell of a regular 1x1 m grid is determined, from which the total mass in each cell, \(M_t\), is calculated. The concentration in each cell is calculated as
$$ C_{cell}=\frac{M_t}{\theta V_{source} } $$The results are compared with BRICK (Neville, 2006), a program that implements the analytical solution to three-dimensional transport of an instantaneous brick source within steady, uniform, unidirectional groundwater flow. The mod-PATH3DU and BRICK results compare very well as shown in .
mod-PATH3DU can also simulate the effect of linear-sorption on the apparent velocity of a reactive solute. This option is demonstrated using the model described in the previous example, Example 4a, as BRICK also simulates linear-sorption via a retardation factor. In this example, a retardation factor of 1.5 is specified for both models. From Example 4a, the center-of-mass for an unretarded solute after 1000 days is 250 m from the source in the x-direction; retarded, this position is 166.67 m from the source (216.67 m) - the average x-coordinate in the mod-PATH3DU output for this simulation is 216.97 m. A total of 100,000 particles, \(N_{particles}\) were released, each assigned a mass, \(M_p\),
$$ M_p=\frac{C_0 \theta V_{source} R}{N_{particles}} $$of 0.018 μg. The resulting concentration distribution determined from the frequency of particle endpoints calculated by mod-PATH3DU for a regular 1x1 m grid is compared with the corresponding BRICK solution in the figure below. The mod-PATH3DU and BRICK results compare very well.
Example files will be available soon.
This example illustrates the advection-dispersion functionality of mod-PATH3DU in a non-uniform groundwater flow field. It is a single layer model with 201 rows and 201 columns. Model cells are 100 feet a side. Hydraulic conductivity is uniform at 100 feet/day. Constant-head boundary conditions line the left and right-hand sides; along these boundaries specified-heads vary linearly - on the left between 80 feet at the top and 75 feet at the bottom, and on the right between 65 and 55 feet at the top and bottom, respectively. The model is steady-state with three extraction wells; the first pumping at a rate of 25,000 feet3/day, the second at 100,000 feet3/day, and the third at 25,000 feet3/day.
The transport of a contaminant following an instantaneous injection of 1000 ug/feet3 in 9 model cells is simulated using MT3D. Longitudinal dispersivity is 100 feet, and transverse dispersivity is 10 feet. Porosity is 0.3. mod-PATH3DU is used to simulate advection and dispersion by releasing 144 particles over the source area and repeating each particle 5000 times for a total of 720,000 particles. The correspdoning mass for each particle is calculated; and using the particle count in each model cell at different time intervals, following an endpoint simulation, concentration distributions are calculated and visually compared against the MT3D results in the figures that follow. The first set of figures result from tracking on the MODFLOW solution using the Waterloo method, while the second set of figures result from tracking on the equivalent MEUK model (by exporting the MODFLOW head-solution to the required ArcMap ASCII grid format) using the SSP&A method. For both the MODFLOW and equivalent MEUK model, the correspondence with MT3D is good - and certainly demonstrates the enhanced transport visualization capabilities of the Random Walk implementation.
Example files will be available soon.
This example is based on the same MODFLOW and MT3D models described in Example 5a, with the following changes, 1) no extraction wells were simulated, 2) hydraulic conductivity is varied through zones, 3) transverse dispersivity was tripled from 10 to 30, and 4) fewer particles were tracked - only 2,000 per start location for a total of 288,000 particles.
The following figure illustrates the good correspondence between MT3D and the advection-dispersion tracking of mod-PATH3DU. The hydraulic conductivity zone with a value of 10 feet/day is outlined in red to make clear transport in this zone. However, particles incorrectly stagnate at the interface between the hydraulic conductivity contrast. We are currently working to resolve this issue and improve the robustness of our implementation of the ADE.
Example files will be available soon.
This example is taken from MODPATH Version 7: Example Problems (Pollock, 2017). The model has 3 layers, the second layer is a confining unit. A well is pumping in layer 3. The grid is refined (quadpatch) in the area of the well. The model is detailed in Pollock (2017). The example has two parts - first, backward-tracking pathlines for particles released around the well are calculated, second backward-tracking endpoints are calculated. and are the pathlines and endpoints, respectively, calculated by mod-PATH3DU - these figures are comparable to Figures 18 and 20 in Pollock (2017). Three dimensional backward-tracking pathlines are given in . This figure was produced with Groundwater Desktop (GWD).
This example is taken from MODFLOW 6 - Example Problems (MODFLOW 6 Development Team). The purpose of this example is to demonstrate how a triangular mesh can be used to simulate a circular island. It is a synthetic example and provided here for illustrative purposes. The results are not verified beyond a visual inspection of the pathlines and that they radiate perpindicular to flow as expected and shown on and . These figure were produced using Groundwater Desktop (GWD).
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.
Langevin, C.D., Hughes, J.D., Banta, E.R., Niswonger, R.G., Panday, S., and Provost, A.M., 2017, Documentation for the MODFLOW 6 Groundwater Flow Model. U.S. Geological Survey Techniques and Methods, book 6, chap. A55, 197 p.
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.