docs(PRT): revise mf6io and supptechinfo (#1833)

Incorporate fixes and suggestions from @rbwinst-usgs' review.

Co-authored-by: aprovost-usgs <[email protected]>

doc/MODFLOW6References.bib

@@ -3018,7 +3018,7 @@ @article{zhang2012
 	Year = {2012}}
 
 @article{cordes1992,
-	Author = {Cordes C and Kinzelbach W},
+	Author = {Cordes, C and Kinzelbach, W},
 	Date-Added = {2024-05-15 00:00:00 +0000},
 	Doi = {10.1029/92WR01686},
 	Journal = {Water Resources Research},
@@ -3030,7 +3030,7 @@ @article{cordes1992
 @article{pollock2015,
 	Author = {Pollock, D W},
 	Date-Added = {2024-05-15 00:00:00 +0000},
-	Doi = {10.1111/gwat.123286},
+	Doi = {10.1111/gwat.12328},
 	Journal = {Groundwater},
 	Title = {Extending the MODPATH algorithm to rectangular unstructured grids},
 	Volume = {54},

doc/SuppTechInfo/particle-tracking-model.tex

@@ -22,14 +22,14 @@ \subsection{Tracking Approach} \label{sec:trackingapproach}
 
 \noindent On DISV grids, which can contain a mix of rectangular, quad-refined, and non-rectangular cells, the PRT Model identifies the geometry of each cell and applies the most efficient tracking method possible. Rectagular cells do not need to be aligned with the model coordinate axes to be recognized as rectangular. Because measures of the cell geometry are compared with numerical tolerances, it is recommended that the cell vertex coordinates be written to double precision in the model input for a DISV grid.
 
-By default, flows associated with stress packages are assumed to be distributed uniformly over the volume of a cell, as in MODPATH 7. Distributed external inflows and outflows are reflected in the cell-cell flows calculated by the GWF Model, but they are not directly involved in determining the normal face velocities used to track a particle through the cell. The user can optionally assign a flow associated with a stress package to any face of the cell. In MODPATH 7, this is done by setting the value of an input parameter called IFACE to a value that corresponds to one of the six faces of a rectangular cell (left, right, back, front, bottom, and top). In the PRT Model, assignment of external flows is done by setting the value of an input parameter called IFLOWFACE to a value that corresponds to one of the cell faces. To accommodate non-rectangular cells, the face numbering scheme in the PRT Model is different from that in MODPATH 7 and is based on clockwise ordering of the cell faces, as described in the \mf input instructions.
+By default, flows associated with stress packages are assumed to be distributed uniformly over the volume of a cell, as in MODPATH 7. Distributed external inflows and outflows are reflected in the cell-cell flows calculated by the GWF Model, but they are not directly involved in determining the normal face velocities used to track a particle through the cell. The user can optionally assign a flow associated with a stress package to any face of the cell. In MODPATH 7, this is done by setting the value of an input parameter called IFACE to a value that corresponds to one of the six faces of a rectangular cell (left, right, front, back, bottom, and top). In the PRT Model, assignment of external flows is done by setting the value of an input parameter called IFLOWFACE to a value that corresponds to one of the cell faces. To accommodate non-rectangular cells, the face numbering scheme in the PRT Model is different from that in MODPATH 7 and is based on clockwise ordering of the cell faces, as described in the \mf input instructions.
 
 \subsection{Generalized Pollock's Method For Non-Rectangular Cells} \label{sec:genpollockmethod}
 
 \cite{zhang2012} review particle tracking methods that ``extend the widely used velocity interpolation algorithms, such as Pollock's algorithm, to more complex geometries." They classify methods according to (1) whether the method refines a cell into subcells, (2) the basis functions used to interpolate velocity, (3) whether the interpolated velocity is continuous throughout a cell, (4) and whether the method is locally conservative. The generalization of Pollock's method used in the PRT Model is new as of this writing, but it is related to the methods reviewed by \cite{zhang2012} and fits into their classification as described below:
 
 \begin{itemize}
-\item The new method used in the PRT Model subdivides a polygonal cell (in plan view) radially into triangular subcells. Each subcell has one edge that coincides with a cell edge and two edges that are internal to the cell. The number of subcells equals the number of faces of the polygonal cell, and all subcells share a vertex at the centroid of cell.
+\item The new method used in the PRT Model subdivides a polygonal cell (in plan view) radially into triangular subcells. Each subcell has one edge that coincides with a cell edge and two edges that are internal to the cell. The number of subcells equals the number of faces of the polygonal cell, and all subcells share a vertex at the centroid of cell. The cell centroid calculated by PRT is not necessarily identical to the cell center specified by the user for a DISV grid.
 \item Cell-face normal velocities are used to calculate a velocity vector at each vertex of each subcell, including the shared cell-centroid vertex. The basis functions used to interpolate velocity within a subcell are relatively ``low-order" in that they are linear, and the normal velocity is assumed to be constant along the subcell edge that coincides with a cell edge. However, they allow the normal velocity to vary along subcell edges that are internal to the cell, making them more flexible than the lowest-order (``RT\textsubscript{0} space") functions described by \cite{zhang2012}.
 \item Interpolated velocity vectors are discontinuous across boundaries between subcells in the sense that the component along the subcell boundary can be different on each side of the boundary. However, the normal component of velocity is continuous across subcell boundaries, as required by continuity of flow. The interpolated velocity is also continuous in the sense used by \cite{zhang2012}, which is that the subcells share a single cell-centroid velocity. Based on their testing of related methods, \cite{zhang2012} note that ``velocity continuity is not as important as local conservation for the purpose of streamline applications," i.e., particle tracking.
 \item The new method is locally conservative. In mathematical terms, this means that the divergence of the velocity field is uniform through each subcell and the same in all subcells, thereby honoring the divergence for the cell as a whole. In practical terms, it means that subcells will not appear to contain external sources or sinks of water or storage effects if the cell as a whole does not contain external sources, sinks, or storage effects. This helps avoid artificial convergence or spreading apart of closely adjacent particle tracks.
@@ -43,7 +43,7 @@ \subsubsection{Calculation of vertex velocities}
 
 For a cell with $N$ polygon faces, the degrees of freedom are $4N$ cell-vertex velocity components and $2$ cell-centroid velocity components, for a total of $4N + 2$ degrees of freedom. The constraints are $2N$ known cell-face normal velocities, $N$ continuity conditions (equality of normal components) along boundaries between subcells, $N - 1$ independent local conservation conditions (equality of divergences), $2$ constraints from setting the cell-centroid velocity to a weighted average of the cell-vertex velocities, and $1$ ``closure constraint," for a total of $4N + 2$ constraints, which matches the number of degrees of freedom.
 
-The closure constraint corresponds to Option 1 described by \cite{zhang2012} and is conceptually equivalent to the closure constraint used in Pollock's subcell method for quad refinement \citep{pollock2015}. It assumes that the velocity field in the cell is irrotational. This is a rigorous assumption provided the porosity and hydraulic conductivity are homogeneous throughout the cell, the hydraulic conductivity is isotropic, and density is uniform. In \mf, the porosity and hydraulic conductivity are homogeneous throughout a cell, but the hydraulic conductivity can be anisotropic and density can allowed to vary. Strictly speaking it is the head gradient, and not the velocity, that is irrotational in the general case. With 10:1 anisotropy, tests by \cite{zhang2012} of a low-order, locally conservative method with quadrilateral subcells showed ``negligible differences between the different closure constraints" they describe, which represent varying degrees of mathematical rigor. It is possible to incorporate a more mathematically rigorous closure constraint such as Option 2 presented by \cite{zhang2012}, which is due to \cite{cordes1992}, into the calculation of vertex velocities, but that has not been done in the version of the PRT Model described here.
+The closure constraint corresponds to Option 1 described by \cite{zhang2012} and is conceptually equivalent to the closure constraint used in Pollock's subcell method for quad refinement \citep{pollock2015}. It assumes that the velocity field in the cell is irrotational. This is a rigorous assumption provided the porosity and hydraulic conductivity are homogeneous throughout the cell, the hydraulic conductivity is isotropic, and density is uniform. In \mf, the porosity and hydraulic conductivity are homogeneous throughout a cell, but the hydraulic conductivity can be anisotropic and density can be allowed to vary. Strictly speaking it is the head gradient, and not the velocity, that is irrotational in the general case. With 10:1 anisotropy, tests by \cite{zhang2012} of a low-order, locally conservative method with quadrilateral subcells showed ``negligible differences between the different closure constraints" they describe, which represent varying degrees of mathematical rigor. It is possible to incorporate a more mathematically rigorous closure constraint such as Option 2 presented by \cite{zhang2012}, which is due to \cite{cordes1992}, into the calculation of vertex velocities, but that has not been done in the version of the PRT Model described here.
 
 \subsubsection{Velocity Interpolation and Solution of the Particle Trajectory}
 
@@ -61,7 +61,7 @@ \subsubsection{Velocity Interpolation and Solution of the Particle Trajectory}
 \label{eqn:vylinear}
 \end{equation}
 
-\noindent where $t$ is time ($T$), subscripted $v$ indicates a velocity component ($L/T$) in the $\tilde{x}$ or $\tilde{y}$ direction, and $\tilde{x}$ and $\tilde{y}$ are Cartesian coordinates ($L$) with their origin at one of the subcell vertices that corresponds to a cell vertex, and rotated such that the $\tilde{y} = 0$ axis coincides with the cell face. Because the normal velocity is assumed to be uniform along the cell face, the coefficient $w_{\tilde{y} \tilde{x}}$ in equation \ref{eqn:vylinear} is zero, eliminating the dependence of equation \ref{eqn:vylinear} on $\tilde{x}$ and allowing it to be integrated analytically, independently of equation \ref{eqn:vxlinear}. The resulting expression for $\tilde{y}$ as a function of time can then be substituted into equation \ref{eqn:vxlinear} to yield an equation that can, once again, be integrated analytically to yield an expression for $\tilde{x}$ as a function of time. The particular forms of the analytical expressions depend on the values of and relationships between the coefficients in equations \ref{eqn:vxlinear} and \ref{eqn:vylinear}, and there are a number of special cases to consider. In any case, equipped with the appropriate expressions, one could guess a value of $t$, analytically calculate the corresponding values of $\tilde{x}$ and $\tilde{y}$, and continue updating one's guesses for $t$ until the point ($\tilde{x}$, $\tilde{y}$) landed on a subcell boundary, at which point the exit time and location of the particle would be known.
+\noindent where $t$ is time ($T$), subscripted $v$ indicates a velocity component ($L/T$) in the $\tilde{x}$ or $\tilde{y}$ direction, and subscripted $w$ indicates the constant rate ($1/T$) at which a velocity component changes with $\tilde{x}$ or $\tilde{y}$. The coordinates $\tilde{x}$ and $\tilde{y}$ are Cartesian coordinates ($L$) with their origin at one of the subcell vertices that corresponds to a cell vertex, and rotated such that the $\tilde{y} = 0$ axis coincides with the cell face. Because the normal velocity is assumed to be uniform along the cell face, the coefficient $w_{\tilde{y} \tilde{x}}$ in equation \ref{eqn:vylinear} is zero, eliminating the dependence of equation \ref{eqn:vylinear} on $\tilde{x}$ and allowing it to be integrated analytically, independently of equation \ref{eqn:vxlinear}. The resulting expression for $\tilde{y}$ as a function of time can then be substituted into equation \ref{eqn:vxlinear} to yield an equation that can, once again, be integrated analytically to yield an expression for $\tilde{x}$ as a function of time. The particular forms of the analytical expressions depend on the values of and relationships between the coefficients in equations \ref{eqn:vxlinear} and \ref{eqn:vylinear}, and there are a number of special cases to consider. In any case, equipped with the appropriate expressions, one could guess a value of $t$, analytically calculate the corresponding values of $\tilde{x}$ and $\tilde{y}$, and continue updating one's guesses for $t$ until the point ($\tilde{x}$, $\tilde{y}$) landed on a subcell boundary, at which point the exit time and location of the particle would be known.
 
 To facilitate the identification of possible exit faces and bounding of possible solutions, and to streamline the solution procedure, the generalized Pollock's method transforms equations \ref{eqn:vxlinear} and \ref{eqn:vylinear} into nonorthogonal coordinates $\alpha$ and $\beta$ that correspond to barycentric coordinates on the triangular subcell. The coordinate transformation is linear and results in a pair of ordinary differential equations of the same general form as equations \ref{eqn:vxlinear} and \ref{eqn:vylinear}, with $\alpha$ and $\beta$ taking the places of $\tilde{x}$ and $\tilde{y}$, respectively. The equation for $\beta$ is effectively decoupled from the equation for $\alpha$, which ultimately allows $\alpha$ and $\beta$ to be expressed explicitly as functions of time.
 
@@ -94,7 +94,7 @@ \subsection{Summary of Main Assumptions and Limitations}
 \item ``Backward tracking" is not supported explicitly. However, backward tracking can be performed by preprocessing the binary budget and head files to reverse the chronological order of the flow and head information. FloPy \citep{bakker2016flopy, hughes2023flopy, flopysoftware} provides utilities for this.
 \item For non-rectagular cells, a generalization of Pollock's method is used. This method is based on refinement of a polygon cell (in plan view) into triangular subcells, with a linear variation of velocity within each subcell that honors the cell-face normal flows generated by the associated GWF Model. Exit location and time are generally solved using an iterative, one-dimensional root-finding procedure, making the method semi-analytical.
 \item Cell-face normal velocities are assumed to be uniform along each cell face, as in Pollock's original method. This ``low-order" approximation has the advantage of preventing particles from crossing no-flow faces at the outer boundary of a model.
-\item As in the methods described by \cite{zheng2012}, variations in fluid density are not accounted for when enforcing local conservation.
+\item As in the methods described by \cite{zhang2012}, variations in fluid density are not accounted for when enforcing local conservation.
 \item The ``closure constraint" used in calculating cell-vertex velocities is approximate in the case of anisotropic hydraulic conductivity and/or variable density. It is conceptually equivalent to the closure closure constraint used in Pollock's subcell method for quad refinement \citep{pollock2015}.
 \item For a triangular cell, the generalized Pollock's method could be used to calculate the particle trajectory without refinement into triangular subcells. This would improve computational efficiency for triangular cells and will likely be programmed in a future version of the PRT Model.
 \end{itemize}

doc/mf6io/framework/binaryoutput.tex

@@ -771,7 +771,7 @@ \subsection{Particle Track File}
 \noindent Field 12: \texttt{`X'} {\color{red} \footnotesize{DOUBLE}} \\
 \noindent Field 13: \texttt{`Y'} {\color{red} \footnotesize{DOUBLE}} \\
 \noindent Field 14: \texttt{`Z'} {\color{red} \footnotesize{DOUBLE}} \\
-\noindent Field 15: \texttt{`NAME'} {\color{red} \footnotesize{CHARACTER(LEN=LENTXT)}} \\
+\noindent Field 15: \texttt{`NAME'} {\color{red} \footnotesize{CHARACTER(LEN=40)}} \\
 
 \vspace{4mm}
 \noindent The ``NAME'' field may be empty.

doc/mf6io/mf6ivar/dfn/prt-mip.dfn

@@ -38,4 +38,4 @@ reader readarray
 layered true
 optional true
 longname zone number
-description is an integer zone number assigned to each cell.  IZONE may be positive, negative, or zero.  The current cell's zone number is recorded with each particle track datum.  If the ISTOPZONE option is set to any value other than zero in a PRP Package, particles released by that PRP Package terminate if they enter a cell whose IZONE value matches ISTOPZONE.  If ISTOPZONE is not specified or is set to zero in a PRP Package, IZONE has no effect on the termination of particles released by that PRP Package.
+description is an integer zone number assigned to each cell.  IZONE may be positive, negative, or zero.  The current cell's zone number is recorded with each particle track datum.  If a PRP package's ISTOPZONE option is set to any value other than zero, particles released by that PRP Package terminate if they enter a cell whose IZONE value matches ISTOPZONE.  If ISTOPZONE is not specified or is set to zero in a PRP Package, IZONE has no effect on the termination of particles released by that PRP Package.  Each PRP Package may configure a single ISTOPZONE value.

doc/mf6io/mf6ivar/dfn/prt-oc.dfn

@@ -96,7 +96,7 @@ reader urword
 tagged true
 optional false
 longname track keyword
-description keyword to specify that record corresponds to a binary track file.
+description keyword to specify that record corresponds to a binary track file.  Each PRT Model's OC Package may have only one binary track output file.
 
 block options
 name trackfile
@@ -108,7 +108,7 @@ reader urword
 tagged false
 optional false
 longname file keyword
-description name of the output file to write tracking information.
+description name of the binary output file to write tracking information.
 
 block options
 name trackcsv_filerecord
@@ -129,7 +129,7 @@ reader urword
 tagged true
 optional false
 longname track keyword
-description keyword to specify that record corresponds to a CSV track file.
+description keyword to specify that record corresponds to a CSV track file.  Each PRT Model's OC Package may have only one CSV track file.
 
 block options
 name trackcsvfile