diff --git a/doc/mf6io/gwe/gwe.tex b/doc/mf6io/gwe/gwe.tex index b314eef36fc..8c96986c8dc 100644 --- a/doc/mf6io/gwe/gwe.tex +++ b/doc/mf6io/gwe/gwe.tex @@ -1,153 +1,144 @@ -Like GWT \citep{modflow6gwt}, the GWE Model simulates three-dimensional transport in flowing groundwater. The primary difference between GWT and GWE is that heat (i.e., temperature), instead of concentration, is the simulated ``species.'' As such, the GWE Model solves the heat transport equation using numerical methods and a generalized control-volume finite-difference approach, which can be used with regular MODFLOW grids (DIS Package) or with unstructured grids (DISV and DISU Packages). The GWE Model is designed to work with most of the new capabilities released with the GWF Model, including the Newton flow formulation, XT3D \citep{modflow6xt3d}, unstructured grids, advanced packages, and the movement of water between packages. The GWF and GWE (and, if active, GWT) models operate simultaneously during a \mf simulation to represent coupled groundwater flow and heat transport. The GWE Model can also run separately from a GWF Model by reading the heads and flows saved by a previously run GWF Model. The GWE model is also capable of working with the flows from another groundwater flow model as long as the cell-by-cell and boundary flows and groundwater heads are written to ``linker'' files in the correct format. -The purpose of the GWE Model is to calculate changes in groundwater temperature in both space and time. Groundwater temperature within an aquifer can change in response to different energy transport processes. These processes include (1) convective (advective) transport of heat with flowing groundwater, (2) the combined hydrodynamic dispersion processes of velocity-dependent mechanical dispersion and conduction (analogous to chemical diffusion), (3) thermal equilibrium with the aquifer matrix, (4) mixing of fluids from groundwater sources and sinks, and (5) direct addition of thermal energy. +The Groundwater Energy (GWE) Model for \mf simulates three-dimensional transport of thermal energy in flowing groundwater based on a generalized control-volume finite-difference approach. The GWE Model is designed to work with the Groundwater Flow (GWF) Model \citep{modflow6gwf} for \mf, which simulates transient, three-dimensional groundwater flow. The version of the GWE model documented here must use the same spatial discretization used by the GWF Model; however, that spatial discretization can be represented by regular MODFLOW grids consisting of layers, rows, and columns, or by more general unstructured grids. The GWE Model simulates (1) advective transport, (2) the combined hydrodynamic dispersion processes of velocity-dependent mechanical dispersion and thermal conduction in groundwater, (3) thermal conduction in the solid aquifer material, (4) storage of thermal energy in the groundwater and solid aquifer material, (5) thermal equilibration between the groundwater and solid aquifer material, (5) zero-order decay or production of thermal energy in the groundwater and the solid, (6) mixing from groundwater sources and sinks, and (7) direct addition or removal of thermal energy. The GWE Model can also represent advective energy transport through advanced package features, such as streams, lakes, multi-aquifer wells, and the unsaturated zone. If the GWF Model that provides the flow field for a GWE Model uses the Water Mover (MVR) Package to connect flow packages, then energy transport between these packages can also be represented by activating the Mover Energy Transport (MVE) Package. The transport processes described here have been implemented in a fully implicit manner and are solved in a system of equations using iterative numerical methods. The present version of the GWE Model for \mf does not have an option to calculate steady-state transport solutions; if a steady-state solution is required, then transient evolution of the temperature must be represented using multiple time steps until no further significant changes in temperature are detected. Finally, temperatures calculated by a GWE Model may be used in the Buoyancy (BUY) and Viscosity (VSC) Packages. Instructions provided in the \mf Description of Input and Output document explain how pass GWE-calculated temperatures to BUY and VSC. -For GWE, the temperature at any point in the aquifer is assumed to instantaneously equilibrate between the aqueous and solid phase domains. For example, a pulse of heat convecting through an aquifer will be retarded through thermal equilibration with the aquifer material. Conversely, the introduction of cold groundwater into a previously warm region of the aquifer will warm up, at least in part, as energy within the aquifer matrix transfers to the aqueous phase. Unlike GWT, the GWE Model type does not support an immobile domain. The energy that is transferred between the aqeous and solid phases of the groundwater system are tracked in the GWE Model budget. +Transport of thermal energy is sometimes simulated using a solute-transport model. In this approach, which takes advantage of the close analogy between the governing equations for thermal-energy transport and solute transport, the solute-transport model is run with values of input parameters that have been scaled to render the governing equation representative of thermal-energy transport \citep{thorne2006, hechtmendez, mazheng2010}. The GWT Model for \mf can be used in this way. However, the new GWE Model offers the following advantages: model input and output in terms of parameters and in units that relate directly to thermal-energy transport; variation of bulk thermal conductivity with changing saturation; and simulation of thermal conduction through the solid matrix, even in dry cells. -This section describes the data files for a \mf Groundwater Energy Transport (GWE) Model. A GWE Model is added to the simulation by including a GWE entry in the MODELS block of the simulation name file. There are three types of spatial discretization approaches that can be used with the GWE Model: DIS, DISV, and DISU. The input instructions for these three packages are not described here in this section on GWE Model input; input instructions for these three packages are described in the section on GWF Model input. +\subsection{Mathematical Model} \label{sec:mathmodel} -The GWE Model is designed to permit input to be gathered, as it is needed, from many different files. Likewise, results from the model calculations can be written to a number of output files. The GWE Model Listing File is a key file to which the GWE model output is written. As \mf runs, information about the GWE Model is written to the GWE Model Listing File, including much of the input data (as a record of the simulation) and calculated results. Details about the files used by each package are provided in this section. +The Groundwater Energy (GWE) Model solves the following governing equation for thermal energy transport: -The GWE Model reads a file called the Name File, which specifies most of the files that will be used in a groundwater energy transport simulation. Several files are always required whereas other files are optional depending on the question(s) being addressed by the model. The Output Control Package receives instructions from the user to control the amount and frequency of output. Details about the Name File and the Output Control Package are described in this section. +\begin{equation} +\label{eqn:pdegwe} +\begin{aligned} +\rho_w C_{pw} \frac {\partial \left ( S_w \theta T \right )}{\partial t} ++ \rho_s C_{ps} \frac {\partial \left [ \left (1 - \theta \right ) T \right ]}{\partial t} += - \rho_w C_{pw} \nabla \cdot \left ( \matr{q} T \right ) ++ \rho_w C_{pw} \nabla \cdot \left ( \matr{D}_T \nabla T \right )\\ ++ \rho_w C_{pw} q'_s T_s + E_s + E_a - \gamma_{1w} \theta S_w +- \gamma_{1s} \left (1 - \theta \right ) \rho_s, +\end{aligned} +\end{equation} -For the GWE Model, ``flows'' (unless stated otherwise) represent the ``flow'' of energy, often expressed in units of energy (e.g., joules) per time, rather than groundwater flow. +\noindent The two terms on the left-hand side of equation \ref{eqn:pdegwe} represent the rates of change in thermal energy storage in the water and the solid matrix, respectively. The six terms on the right-hand side of equation \ref{eqn:pdegwe} represent the rates at which thermal energy is advected, dispersed and/or conducted, added or removed by groundwater inflows or outflows, added or removed directly to or from the water, produced/decayed in the water, and produced/decayed in the solid matrix, respectively. The parameters and variables in equation \ref{eqn:pdegwe} are defined as follows: $S_w$ is the water saturation (dimensionless) defined as the volume of water per volume of voids; $\theta$ is the effective porosity, defined as volume of voids participating in transport per unit volume of aquifer; $\rho_w$ and $\rho_s$ are the densities ($M/L^3$) of the water and solid-matrix material, respectively; $C_{pw}$ and $C_{ps}$ are the specific heat capacities ($E/(M \, deg)$) of the water and solid-matrix material, respectively; $T$ is temperature ($deg$); $t$ is time ($T$); $\matr{q}$ is the vector of specific discharge ($L/T$); $\matr{D}_T$ is the second-order tensor of hydrodynamic dispersion coefficients for thermal energy transport ($L^2/T$); $q'_s$ is the volumetric flow rate per unit volume of aquifer (defined as positive for flow into the aquifer) for fluid sources and sinks ($1/T$), $T_s$ is the temperature of the source or sink fluid ($deg$), $E_s$ is rate of energy loading per unit volume of aquifer ($M/L^3T$), $E_a$ is rate of energy exchange with advanced stress packages ($M/L^3T$), $\gamma_{1w}$ is the zero-order energy decay rate coefficient in the water ($E/L^3T$), and $\gamma_{1s}$ is the zero-order energy decay rate coefficient in the solid ($E/MT$). Note that $\gamma_{1w}$ is defined on a per-volume-of-water basis, whereas $\gamma_{1s}$ is defined on a per-mass-of-solid basis. Note that $\rho_w$, $\rho_s$, $C_{pw}$, and $C_{ps}$ are assumed to remain constant with time, although the solid properties can vary spatially from cell to cell. -\subsection{Information for Existing Heat Transport Modelers} -An important goal of the \mf GWE Model is to alleviate the need for ``parameter equivalents'' when simulating heat transport in groundwater systems. In the past, codes like HST3D \citep{kipp1987} or VS2DH \citep{healy1996} simulated energy transport directly by supporting the use of native heat transport units. For example, users could directly specify thermal conductivity of the fluid and solid phases, as well as the heat capacity of both phases. Alternatively, codes like MT3DMS \citep{zheng1999mt3dms}, MT3D-USGS \citep{mt3dusgs}, and MODFLOW-USG \citep{modflowusg} could be used to simulate the movement of heat in groundwater, but required users to leverage existing variables as surrogates for heat transport. For example, the molecular diffusion parameter may be used as a surrogate for simulating thermal conduction in an aquifer \citep{mazheng2010, hechtmendez}. +Equation \ref{eqn:pdegwe} is closely analogous to the equation solved by the Groundwater Transport (GWT) Model for solute transport \citep{modflow6gwt}, with the following differences: -The following list summarizes important aspects of GWE for simulating heat transport with \mf: +\begin{itemize} +\item The GWT solute-transport equation describes the movement and conservation of solute mass, with the ``solute-mass density" (solute mass per volume of aquifer) represented by the volumetric concentration $C$ [$M/L^3$]. Equation \ref{eqn:pdegwe} describes the movement and conservation of thermal energy, with the ``thermal-energy density" (energy per volume of aquifer) represented by $\rho C_{p} T$ [$E/L^3$]. +\item Equation \ref{eqn:pdegwe} is based on the assumption that at any point in the model, the groundwater and solid matrix through which it moves are at the same temperature, i.e., they are at thermal equilibrium. The second term on the left-hand side of equation \ref{eqn:pdegwe} represents the rate of change in thermal energy stored in the solid. This term is analogous to the equilibrium-sorption term in the solute-transport equation, which represents the rate of change in solute mass sorbed onto the solid. There is no separate ``immobile domain" represented in the GWE Model, since all the material within a cell is assumed to be at the same temperature. +\item The GWT Model allows both zero-order (independent of concentration) and first-order (proportional to concentration) production/decay of solute and sorbate mass in the water to represent simple chemical reactions. The GWE model allows only zero-order (independent of temperature) production/decay of thermal energy, but it can occur in both the water and the solid. +\item Like the GWT solute-transport equation, equation \ref{eqn:pdegwe} includes a dispersion term. In the GWT solute-transport equation, the dispersion term represents mechanical dispersion and molecular diffusion of solute mass in the water. In equation \ref{eqn:pdegwe}, the dispersion term represents mechanical dispersion and conduction of thermal energy in the water, as well as conduction of thermal energy in the solid. Conduction in the solid has no analog in the GWT solute-transport equation. +\end{itemize} -\begin{enumerate} +\noindent The GWE Model gives \mf a thermal-energy-transport simulation capability that is comparable to that in the finite-difference-based simulation codes HST3D \citep{kipp1987}, VS2DH \citep{healy1996}, SEAWAT 4 \citep{langevin2008seawat}, the Block-Centered Transport (BCT) Process for MODFLOW-USG \citep{modflowusg, panday2019bct}, and MT3D-USGS \citep{mt3dusgs}. -\item The GWE Model uses parameters that are native to heat transport, including thermal conductivity of water, heat capacity of water, thermal conductivity of the aquifer material, heat capacity of of the aquifer material, and latent heat of vaporization. Therefore, users do not need to pre-calculate solute-transport ``parameter equivalents'' when generating GWE model input; users can instead enter native parameter values that are readily available. +\subsection{Energy-Transport Packages} \label{sec:packages} -\item Thermal energy transport budgets written to the \mf list file are reported in units of energy (e.g., joules). Previously, using a program like MT3D-USGS \citep{mt3dusgs} to emulate heat transport using solute transport, units in the list file budget did not correspond to thermal energy, but were reported in units of $\frac{m^{3 \;\circ}C}{d}$. To convert to thermal energy units, values in the list file had to be post-processed by multiplying each line item by the density of water ($\rho_w$) and the heat capacity of water ($C_p$) \citep{langevin2008seawat}. +\noindent Equation~\ref{eqn:pdegwe} can be rewritten in the form -\item Thermal equilibrium between the aqueous and solid phases is assumed. Thus, simulated temperatures are representative of both phases. As a result, thermal conduction between adjacent cells may still occur even in the absence of convection. +\begin{equation} +\label{eqn:pdegwefterms} +\begin{aligned} +f^{EST} + f^{ADV} + f^{CND} + f^{SSM} + f^{ESL} + f^{APT} = 0, +\end{aligned} +\end{equation} -\item In GWE, dry cells (devoid of groundwater) remain active for simulating thermal conduction. For example, energy (heat) transfer will be simulated between a partially saturated cell (i.e., ``water-table'' cell) and an overlying dry cell. In this way, a more full accounting of various heat transport processes is represented in the subsurface. Moreover, this approach readily supports heat transport in the unsaturated-zone when the UZE (unsaturated-zone energy transport) Package is active. +\noindent where -\item Heat transport is supported for all five of the advanced GWF packages using the following packages in GWE: (1) streamflow energy transport, SFE Package; (2) lake energy transport, LKE Package; (3) multi-aquifer well energy transport, MWE Package; (4) unsaturated zone energy transport, UZE Package; and the (5) Water Mover Package, MVE. Similar to GWT, GWE will simulate heat transfer between an advanced package and the groundwater system via groundwater surface-water exchange; however, GWE also simulates a conductive transfer of heat between an advanced package feature and the aquifer. To take advantage of this functionality, users must specify the thermal conductivity of the material separating a stream from the aquifer, for example, the thermal conductivity of the streambed (or lakebed), as well as the thickness of the streambed (or lakebed). As with the advanced GWT packages, GWE simulates thermal convection between package features, such as between two stream reaches for example. Also, dispersive heat transport among among advanced package features is not represented, similar to GWT. +\begin{equation} +\label{eqn:fterms} +\begin{aligned} +& f^{EST} = - \rho_w C_{pw} \frac {\partial \left ( S_w \theta T \right )}{\partial t} +- \rho_s C_{ps} \frac {\partial \left [ \left (1 - \theta \right ) T \right ]}{\partial t} +- \gamma_{1w} \theta S_w - \gamma_{1s} \left (1 - \theta \right ) \rho_s \\ +& f^{ADV} = - \rho_w C_{pw} \nabla \cdot \left ( \matr{q} T \right ) \\ +& f^{CND} = \rho_w C_{pw} \nabla \cdot \left ( \matr{D}_T \nabla T \right ) \\ +& f^{SSM} = \rho_w C_{pw} q'_s T_s \\ +& f^{ESL} = E_s \\ +& f^{APT} = E_a, +\end{aligned} +\end{equation} -\item Where the GWF model simulates evaporation from an open body of water, for example from the surface of a stream or lake, the latent heat of vaporization may be used to simulate evaporative cooling. As water is converted from liquid to gas, the energy required by the phase change is drawn from the remaining body of water and the resulting cool down is calculated. +\noindent and the three-letter acronyms in superscript denote the GWE Model packages that handle the corresponding terms of the thermal-energy equation. The term $ f^{APT}$ represents energy transfer between the groundwater model and advanced stress packages. Depending on the specific advanced package, mechanisms of energy transfer can include advection, conduction, and thermal equilibration. (Note that an analogous term should appear in equation 2--4 of \cite{modflow6gwt} to represent transfer of solute mass between a groundwater model and advanced stress packages of a GWT model, but was omitted from that equation.) -\end{enumerate} +Each of the GWE Model packages described here is analogous to a package associated with the GWT Model \citep{modflow6gwt}. Differences are highlighted below. -Many of the same considerations listed for the GWT model should be kept in mind when developing a GWE model. For convenience, many of those considerations are adapted for GWE and repeated here. +\subsubsection{GWE Packages Directly Analogous to GWT Packages} -\begin{enumerate} +The following GWE Model packages are directly analogous to packages of the same name in the GWT Model and will not be described in detailed here: Advection (ADV), Source-Sink Mixing (SSM), Mover Transport (MVT), Initial Conditions (IC), Flow Model Interface (FMI), Model Observations (OBS), and Output Control (OC). In the \mf source code, each of these packages is implemented as a single module that can be used for both GWT and GWE models. From the user's perspective, the only difference is that the GWE Model packages deal with transport of thermal energy instead of solute mass, and their input data are accordingly in units representative of energy transport. In addition, the GWE Model uses the Spatial Discretization (DIS, DISV, or DISU) Packages, which are also used by the GWT and Groundwater Flow (GWF) Models. -\item A GWE Model can access flows calculated by a GWF Model that is running in the same simulation as the GWE Model. Alternatively, a GWE Model can read binary head and budget files created from a previous GWF Model simulation (provided these files contain all of the required information for all time steps); there is no specialized flow and transport link file \citep{zheng2001modflow} as there is for MT3D. Details on these two different use cases are provided in the chapter on the FMI Package. +The Constant Temperature (CTP) and Energy Source Loading (ESL) Packages in the GWE Model are directly analogous to the Constant Concentration (CNC) and Mass Source Loading (SRC) Packages in the GWT Model. Although also they function very similarly to their GWT Model counterparts, the CTP and ESL Packages have been assigned distinct names that better represent their connections with temperature and energy rather than concentration and solute mass. -\item The GWE Model is based on a generalized control-volume finite-difference method, which means that heat transport can be simulated using regular MODFLOW grids consisting of layers, rows, and columns, or heat transport can be simulated using unstructured grids. +\subsubsection{Energy Storage and Transfer (EST) Package} -\item GWE and GWT use the same advection package source code. As a result, advection can be simulated using central-in-space weighting, upstream weighting, or an implicit second-order TVD scheme. Currently, neither the GWE or GWT models can use a Method of Characteristics (particle-based approaches) or an explicit TVD scheme to simulate convective (or advective) transport. Consequently, the GWE Model may require a higher level of spatial discretization than other transport models that use higher order terms for advection dominated systems. This can be an important limitation in problems involving sharp heat fronts. +The Energy Storage and Transfer (EST) Package in the GWE Model is analogous to the Mobile Storage and Transfer (MST) Package in the GWT Model. The EST Package accounts for changes in thermal energy storage in the groundwater and the solid matrix under the assumption of thermal equilibrium between water and solid, as well as zero-order production/decay of thermal energy in the water and the solid. Thermal equilibration is analogous to equilibrium sorption. In the MST package, sorption is optional, whereas thermal equilibrium is integral to the formulation of the GWE Model. -\item The Viscosity Package can use temperatures from the GWE model to adjust the viscosities in the flow model. +Because the temperature is assumed to be the same in all material within a cell, there is no need to simulate the effects of temperature changes in a separate ``immobile domain.'' Consequently, there is no analog in the GWE Model to the Immobile Storage and Transfer (IST) Package in the GWT Model. -\item GWE and GWT use the same Source and Sink Mixing (SSM) Package for representing the effects of GWF stress package inflows and outflows on simulated temperatures and concentrations. In a GWE simulation, there are two ways in which users can assign concentrations to the individual features in these stress package. The first way is to activate a temperature auxiliary variable in the corresponding GWF stress package. In the SSM input file, the user provides the name of the auxiliary variable to be used for temperature. The second way is to create a special SPC file, which contains user-assigned time-varying temperatures for stress package features. +\subsubsection{Conduction and Dispersion (CND) Package} -\item The GWE model includes an Energy Storage and Transfer (EST) Package that is analogous to the MST Package in the GWT Model. The GWE Model does not simulate immobile domains and therefore does not include an analog of the IST Package in the GWT Model. +The Conduction and Dispersion (CND) Package in the GWE Model is analogous the Dispersion (DSP) Package in the GWT Model. As mentioned earlier in this chapter, the DSP Package accounts only for molecular diffusion and mechanical dispersion (together referred to simply as ``dispersion") in the water, since diffusion of solute into the solid matrix material is typically not of practical interest. By contrast, the CND Package accounts for conduction (thermal diffusion) in both the water and the solid, as well as mechanical dispersion of thermal energy in the water. -\item A GWE-GWE Exchange (introduced in version 6.5.0) can be used to tightly couple multiple heat transport models, as might be done in a nested grid configuration. +The dispersion term handled by the CND Package is given by the third expression listed in equation \ref{eqn:fterms}. The term in parentheses is the negative of the dispersive flux, which includes mechanical dispersion in the water and conduction in the water and solid. Accordingly, the dispersion tensor, $\matr{D}_T$, is expressed as the sum of contributions from mechanical dispersion and conduction: -\item There is no option to automatically run the GWE Model to steady state using a single time step. This is an option available in MT3DMS \citep{zheng2010supplemental}. Steady state conditions must be determined by running the transport model under transient conditions until temperatures stabilize. +\begin{equation} +\label{eqn:dt_tensor} +\begin{aligned} +\matr{D}_T = S_w \theta \matr{D}^{mech} + D^{mol} \matr{I}, +\end{aligned} +\end{equation} -\item As is the case with GWT, the GWE Model has not yet been programmed to work with the Skeletal Storage, Compaction, and Subsidence (CSUB) Package for the GWF Model. +\noindent where $\matr{D}^{mech}$ is the mechanical dispersion tensor, -\end{enumerate} +\begin{equation} +\label{eqn:dmolforheat} +\begin{aligned} +D^{mol} = \frac {{k_T}_{bulk}}{\rho_w C_{pw}} +\end{aligned} +\end{equation} -\subsection{Units of Length and Time} -The GWE Model formulates the groundwater energy transport equation without using prescribed length and time units. Any consistent units of length and time can be used when specifying the input data for a simulation. This capability gives a certain amount of freedom to the user, but care must be exercised to avoid mixing units. The program cannot detect the use of inconsistent units. +\noindent is a molecular diffusion coefficient ($L^2/T$) that combines the conductive capacity of the water and the solid, and -\subsection{Thermal Energy Budget} -A summary of all inflow (sources) and outflow (sinks) of thermal energy is referred to as an energy budget. \mf calculates an energy budget for the overall model as a check on the acceptability of the solution, and to provide a summary of the sources and sinks of energy to the flow system. The energy budget is printed to the GWE Model Listing File for specified time steps. +\begin{equation} +\label{eqn:ktbulk} +\begin{aligned} +{k_T}_{bulk} = S_w \theta {k_T}_w + \left ( 1 - \theta \right ) {k_T}_s +\end{aligned} +\end{equation} -\subsection{Time Stepping} +\noindent is a bulk thermal conductivity ($E/(T L \, deg)$) that depends on the thermal conductivities of the water and the solid, ${k_T}_w$ and ${k_T}_s$ ($E/(T L \, deg)$, respectively. The particular form of the bulk conductivity in equation \ref{eqn:ktbulk} corresponds to the \textit{ad hoc} assumption that the total conductive flux is the sum of independent contributions from the water and the solid, which are proportional to the volumes of water and solid, respectively. One could describe different conductive behavior using a nonlinear expression for bulk conductivity {\citep{campbell1994, markle2006}. However, the bulk conductivity expression in equation \ref{eqn:ktbulk} is the only one available in the version of the GWE Model described here. Note that the CND Package does not account explicitly for tortuosity. The user should adjust the input values of ${k_T}_w$ and ${k_T}_s$ to account for the effects of tortuosity if desired. -For the present implementation of the GWE Model, all terms in the heat transport equation are solved implicitly. With the implicit approach applied to the transport equation, it is possible to take relatively large time steps and efficiently obtain a stable solution. If the time steps are too large, however, accuracy of the model results will suffer, so there is usually some compromise required between the desired level of accuracy and length of the time step. An assessment of accuracy can be performed by simply running simulations with shorter time steps and comparing results. +Once the dispersion tensor is calculated, the CND Package handles dispersion in a way that is directly analogous to the method used in the DSP Package. As in the DSP Package, the XT3D option is the default method for representing dispersion in the CND Package, and XT3D can be turned off by the user in package input. -In \mf time step lengths are controlled by the user and specified in the Temporal Discretization (TDIS) input file. When the flow model and heat transport model are included in the same simulation, then the length of the time step specified in TDIS is used for both models. If the GWE Model runs in a separate simulation from the GWE Model, then the time steps used for the heat transport model can be different, and likely shorter, than the time steps used for the flow solution. Instructions for specifying time steps are described in the TDIS section of this user guide; additional information on GWF and GWE configurations are in the Flow Model Interface section. +\subsubsection{Advanced Stress Packages} +The version of the GWE Model described here offers four advanced stress packages, each of which is analogous and functions similarly to an advanced package in the GWT Model: Lake Energy Transport (LKE), Multi-Aquifer Well Energy Transport (MWE), Streamflow Energy Transport (SFE), and Unsaturated Zone Energy Transport (UZE). The GWE advanced stress packages will not be described in detail here except to note differences from the corresponding GWT packages. The primary difference, of course, is that GWE advanced stress packages explicitly simulate thermal energy transport within their respective features, not solute mass. Optionally, the SSM package may be used to represent the effects of thermal energy exchange between an advanced transport package and the subsurface for cases where the temperature of the feature is known. +The LKE, SFE, and MWE Packages offer a mechanism for energy to conduct between the lake or stream reach and the aquifer through a thermally conductive layer that can represent, for example, a lake or stream bed or a well casing. The LKT, SFT, and MWT packages do not offer analogous diffusion of solute through a conductive layer. The LKE and SFE packages also account for removal of thermal energy from surface water by evaporation. -\newpage -\subsection{GWE Model Name File} -\input{gwe/namefile.tex} +The UZE Package simulates thermal energy transport through the unsaturated zone. It calculates cell temperatures based on cell-cell flows generated by the Unsaturated Zone Flow (UZF) Package in the GWF Model and the assumption of thermal equilibrium between the water and the solid matrix material. Equilibrium sorption, which can be considered a solute-transport analog of thermal equilibrium, is not currently simulated by the corresponding Unsaturated Zone Transport (UZT) Package in the GWT Model. -%\newpage -%\subsection{Structured Discretization (DIS) Input File} -%\input{gwf/dis} +\subsection{Exchanges} -%\newpage -%\subsection{Discretization with Vertices (DISV) Input File} -%\input{gwf/disv} +The GWF-GWE Exchange establishes the link between a GWE model and its corresponding GWF model when the two models are run in the same simulation. The GWE-GWE Exchange establishes the link between two GWE models that have cell-cell connections between the two models. -%\newpage -%\subsection{Unstructured Discretization (DISU) Input File} -%\input{gwf/disu} +\subsection{Transport in Dry Cells} -\newpage -\subsection{Initial Conditions (IC) Package} -\input{gwe/ic} +Because thermal conduction can occur in a dry cell (a cell in which the calculated head is below the cell bottom), model cells are never deactivated in a GWE Model due to going dry, regardless of whether or not the Newton-Raphson formulation is being used. If the Newton-Raphson formulation is \textit{not} used, the only energy-transport process simulated in a dry cell is thermal conduction through the solid matrix. If the Newton-Raphson formulation \textit{is} used, dry cells are assumed to transmit water without storing it, as in the GWT Model. In that case, the energy-transport processes simulated in a dry cell are thermal equilibration between the transmitted water and the solid matrix through which it is passing, thermal conduction through the solid matrix, and advection of thermal energy by the transmitted water. -\newpage -\subsection{Output Control (OC) Option} -\input{gwe/oc} +\subsection{Summary of Main Assumptions and Limitations} -\newpage -\subsection{Observation (OBS) Utility for a GWE Model} -\input{gwe/gwe-obs} +The following is a list of the main assumptions and limitations of the GWE Model for \mf: -\newpage -\subsection{Advection (ADV) Package} -\input{gwe/adv} - -\newpage -\subsection{Conduction and Dispersion (CND) Package} -\input{gwe/cnd} - -\newpage -\subsection{Energy Storage and Transfer (EST) Package} -\input{gwe/est} - -\newpage -\subsection{Source and Sink Mixing (SSM) Package} -\input{gwe/ssm} - -\newpage -\subsection{Constant Temperature (CTP) Package} -\input{gwe/ctp} - -\newpage -\subsection{Energy Source Loading (ESL) Package} -\input{gwe/esl} - -\newpage -\subsection{Streamflow Energy Transport (SFE) Package} -\input{gwe/sfe} - -\newpage -\subsection{Lake Energy Transport (LKE) Package} -\input{gwe/lke} - -\newpage -\subsection{Multi-Aquifer Well Energy Transport (MWE) Package} -\input{gwe/mwe} - -\newpage -\subsection{Unsaturated-Zone Energy Transport (UZE) Package} -\input{gwe/uze} - -\newpage -\subsection{Flow Model Interface (FMI) Package} -\input{gwe/fmi} - -\newpage -\subsection{Mover Energy Transport (MVE) Package} -\input{gwe/mve} - -\newpage -\subsection{Groundwater Energy Transport (GWE) Exchange} -\input{gwe/gwe-gwe} +\begin{itemize} +\item Groundwater is assumed to be in thermal equilibrium with the solid matrix material through which it flows. +\item Densities and specific heats of the water and solid are assumed to be constant in time. However, the solid properties can vary spatially. +\item Like the GWT Model, the GWE Model does not account for changes in porosity with time. Changes in porosity are accounted for implicitly in the storage term of the GWF Model, but \mf does not explicitly track the evolution of porosity due to storage changes. +\item The GWE Model does not simulate freezing and does not check for or prevent calculated temperatures below freezing. +\item In partially saturated cells, the GWE Model does not represent thermal conduction separately in the saturated and dry portions of the cell. A single bulk thermal conductivity is calculated based on the overall water and solid contents of the cell, and conduction between cells is approximated as occurring over the entire shared cell face based on this single conductivity. +\item Dry cells remain active in a GWE model so that thermal conduction through the solid matrix material can be simulated. +\item Terms in the thermal energy transport equation that are associated with groundwater flows are expressed relative to a reference temperature of zero (in whatever temperature unit is being used in the simulation). This choice of reference temperature is straightforward but arbitrary. Choosing a different reference temperature would change the numerical values of all such terms, but the overall energy balance would be unaffected, since the sum of those terms would remain the same as a consequence of groundwater mass conservation. This should be borne in mind when interpreting the reported values of budget terms associated with groundwater flows. +\item A GWE model must use the same spatial discretization as the GWF model from which it obtains groundwater flow information. +\item The present version of the GWE Model does not have an option to directly calculate steady-state transport solutions; a steady-state solution must be approached using a sufficiently long transient simulation. +\end{itemize} diff --git a/doc/mf6io/prt/prt.tex b/doc/mf6io/prt/prt.tex index 385f45f4260..47005382aca 100644 --- a/doc/mf6io/prt/prt.tex +++ b/doc/mf6io/prt/prt.tex @@ -1,168 +1,104 @@ -The PRT Model calculates three-dimensional, advective particle trajectories in flowing groundwater. The PRT Model is designed to work with the Groundwater Flow (GWF) Model \citep{modflow6gwf} and uses the same spatial discretization, which may be represented using either a structured (DIS) or an unstructured (DISV) grid. The PRT Model replicates much of the functionality of MODPATH 7 \citep{pollock2016modpath7} and offers support for a much broader class of unstructured grids. The PRT Model can be run in the same simulation as the associated GWF Model or in a separate simulation that reads previously calculated flows from a binary budget file. The version of the PRT Model documented here does not support grids of DISU type, tracking of particles through advanced stress package features such as lakes or streams reaches, or exchange of particles between PRT models. -This section describes the data files for a \mf Particle Tracking (PRT) Model. A PRT Model is added to the simulation by including a PRT entry in the MODELS block of the simulation name file. There are currently two types of spatial discretization approaches that can be used with the PRT Model: DIS and DISV. The input instructions for these three packages are not described here in this section on PRT Model input; input instructions for these three packages are described in the section on GWF Model input. Note that for a PRT Model, the maximum number of vertices for a cell in a DISV grid is limited to 8. +The Particle Tracking (PRT) Model for \mf calculates three-dimensional, advective particle trajectories in flowing groundwater. The PRT Model is designed to work with the Groundwater Flow (GWF) Model \citep{modflow6gwf} for \mf, which simulates transient, three-dimensional groundwater flow. The PRT Model replicates much of the functionality of MODPATH 7 \citep{pollock2016modpath7} and offers support for a much broader class of unstructured grids. The PRT model uses the same spatial discretization used by the GWF Model, which may be represented using either a structured (DIS) or an unstructured (DISV) grid. The PRT Model can be run in the same simulation as the associated GWF Model or in a separate simulation that reads previously calculated flows from a binary budget file. The version of the PRT Model documented here does not support grids of DISU type, tracking of particles through advanced stress package features such as lakes or streams reaches, or exchange of particles between PRT models. The PRT Model simulates the forward tracking of particles; however, backward particle tracking is possible using a procedure described here for reversing simulated flows from a previous GWF simulation and using those flows as input to PRT. -The PRT Model is designed to permit input to be gathered, as it is needed, from many different files. Likewise, results from the model calculations can be written to a number of output files. Details about the files used by each package are provided in this section on the PRT Model Instructions. +\subsection{Tracking Approach} \label{sec:trackingapproach} -The PRT Model reads a file called the Name File, which specifies most of the files that will be used in a simulation. Several files are always required whereas other files are optional depending on the simulation. The Output Control Package receives instructions from the user to control the amount and frequency of output. Details about the Name File and the Output Control Package are described in this section. +The Particle Tracking Model (PRT) Model solves for advective particle trajectories based on cell-cell flows calculated by the Groundwater Flow (GWF) Model. The overall approach to tracking a particle through a model cell is similar to that in MODPATH 7. -For the PRT Model, ``flows'' (unless stated otherwise) represent particle mass ``flow'' in mass per time, rather than groundwater flow. Each particle is currently assigned unit mass, and the numerical value of the flow can be interpreted as particles per time. +\begin{itemize} +\item On each time step of the simulation, convert cell-cell flows on the faces of the cell that contains the particle to normal velocities (velocity components normal to the cell faces). The normal velocity is assumed to be uniform along each cell face. +\item Interpolate the normal-velocity information to enable calculation of a particle velocity at any point within the cell. +\item Solve analytically or semi-analytically for the particle trajectory through the cell. +\item Continue tracking the particle from cell to cell as necessary until the end of the time step or until the particle exits the model or terminates for another reason. +\end{itemize} -\subsection{Units of Length and Time} -The PRT Model formulates the particle trajectory equations without using prescribed length and time units. Any consistent units of length and time can be used when specifying the input data for a simulation. This capability gives a certain amount of freedom to the user, but care must be exercised to avoid mixing units. The program cannot detect the use of inconsistent units. +The interpolation and solution method used by the PRT Model in a given cell depends on the geometry of the cell in plan view. -\subsection{Time Stepping} -In \mf time step lengths are controlled by the user and specified in the Temporal Discretization (TDIS) input file. When the flow model and particle-tracking model run in the same simulation, the time step length specified in TDIS is used for both models. If the PRT Model runs in a separate simulation, the time discretization may differ. Instructions for specifying time steps are described in the TDIS section of this user guide; additional information on GWF and PRT configurations are in the Flow Model Interface section. By default, the PRT model will terminate particles at the end of the simulation's final time step; PRT can be configured to extend tracking until particles terminate naturally (i.e. at boundary conditions, sinks, or specified stop times) via Particle Release Point (PRP) package settings. +\begin{itemize} +\item For rectangular cells, Pollock's method \citep{pollock2016modpath7} is used. +\item For cells adjacent to quad-refined cells subject to the constraints described in \citep{pollock2016modpath7}, Pollock's subcell method for quad refinement \citep{pollock2015} is used. +\item For all other \mf cells, a generalization of Pollock's method is used. The generalized method is described below. For rectangular cells, it is mathematically equivalent to Pollock's original method. +\end{itemize} -\subsection{Specifying Cell Face Flows using IFLOWFACE} +\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. 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. Assignment of external flows is done by setting the value of an auxiliary package variable called IFLOWFACE to a value that corresponds to one of the cell faces. To accommodate non-rectangular cells, the IFLOWFACE numbering scheme in the PRT Model is based on clockwise ordering of the lateral cell faces. For a DIS-grid cell, IFLOWFACE = 1 corresponds to the ``western'' face, i.e., the face parallel to the y axis that has the lesser x coordinate. For a DISV-grid cell, IFLOWFACE = 1 corresponds to the face that extends in the clockwise direction from the first vertex listed for that cell in the CELL2D input block. IFLOWFACE numbering of lateral cell faces proceeds in clockwise direction. IFLOWFACE = -2 corresponds to the cell bottom, and IFLOWFACE = -1 corresponds to the cell top. +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, as described in the \mf input instructions. -\subsection{Particle Mass Budget} -A summary of all inflow (sources) and outflow (sinks) of particle mass is called a mass budget. The particle mass budget is printed to the PRT Model Listing File for selected time steps. In the current implementation, each particle is assigned unit mass, and the numerical value of the flow can be interpreted as particles per time. +The motion of a particle is determined by the groundwater velocity field in which the particle is immersed. In a fully saturated cell, or the saturated portion of a partially saturated cell, the velocity field is calculated from the flows entering and exiting the cell. In a completely dry cell, or the dry portion of a partially saturated cell, the fate of a particle depends on whether the cell is an active part of the flow simulation, whether the particle is in a dry or wet part of the cell, and a user-selected option. -\subsection{Vertical Tracking} +A cell can be inactive either because it has been removed from the active simulation using the IDOMAIN array or because it is completely dry, i.e., the head in the cell is below the bottom elevation of the cell. Deactivation of completely dry cells is the default behavior in MODFLOW 6. However, when the Newton-Raphson formulation is used to solve for groundwater flow, completely dry cells remain active. Particle tracking through inactive and dry-but-active cells is discussed in detail in the PRT chapter of the \mf Description of Input and Output document. -When a particle is in the flow field, vertical motion can be solved in the same way as lateral motion. Special handling is necessary above the water table. +\subsection{Generalized Pollock's Method For Non-Rectangular Cells} \label{sec:genpollockmethod} -A dry cell is either an inactive cell or an active but dry cell, as can occur with the Newton formulation. +\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: -Normally, an inactive cell might be dry or explicitly disabled (idomain). With Newton, dry cells remain active. +\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. 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. +\end{itemize} -Release-time and tracking-time considerations are described (and implemented) separately. +As in MODPATH 7, movement of a particle in the $z$ (vertical) direction is calculated independently from its movement in the $(x, y)$ (horizontal) plane \citep{pollock2016modpath7} . Therefore, the discussions of vertex velocities and the trajectory solution below focus on two-dimensional motion within the horizontal plane. -\subsubsection{Release} +\subsubsection{Calculation of vertex velocities} -At release time, PRT decides whether to release each particle or to terminate it unreleased. +A single velocity vector ($x$ and $y$ velocity components) could be calculated at each cell vertex such that each vertex velocity respects the uniform normal velocity along each cell faces. A centroid velocity could then be calculated by averaging the vertex velocities. However, this approach would not be locally conservative. The additional degrees of freedom (velocity components to be calculated) that are needed to allow local conservation can be introduced by using higher-order basis functions, but that approach typically requires expensive numerical integration of the particle trajectory \citep{zhang2012}. An alternative is to introduce the requisite additional degrees of freedom by allowing each subcell to have its own distinct velocity vector at the cell centroid, but that has the obvious disadvantage that the velocity is discontinuous at the cell centroid. The method used in the PRT Model is novel in that it uses relatively low-order basis functions that allow for semi-analytical solution of the particle trajectory, together with two velocities at each cell vertex (one for each of the two subcells that share the vertex) and a single cell-centroid velocity to provide sufficient degrees of freedom for local conservation. This approach can also accommodate the discontinuity in cell-face normal velocity that occurs when two cell faces meet at a $180^\circ$ angle, as occurs in quad-refined grids. -If the release cell is active, the particle will be released at the user-specified location. +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. -If the release cell is inactive, behavior is determined by the DRAPE option. If the DRAPE option is enabled, the particle will be released from the top-most active cell beneath it, if any. If there is no active cell underneath the particle in any layer, or if DRAPE is not enabled, the particle will terminate unreleased (with status code 8). +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. -Since under the Newton formulation dry cells can remain active, the DRAPE option has no effect when Newton is on (assuming particles are not released into disabled grid regions). Vertical tracking behavior with Newton can be configured using the DRY\_TRACKING\_METHOD option discussed below. +\subsubsection{Velocity Interpolation and Solution of the Particle Trajectory} -\subsubsection{Tracking} +As in Pollock's method \citep{pollock2016modpath7} for a rectangular cell, the velocity is assumed to vary linearly in space within a triangular subcell: -With the Newton formulation, particles can be released into dry-but-active cells. +\begin{equation} +v_{\tilde{x}} = \frac{d \tilde{x}}{dt} = v_{\tilde{x}}^{(0)} + w_{\tilde{x} \tilde{x}} \tilde{x} + + w_{\tilde{x} \tilde{y}} \tilde{y} +\label{eqn:vxlinear} +\end{equation} -Particle trajectories are solved over the same time discretization used by the flow model. A particle may be immersed in the flow field in one time step, and find that the water table has dropped below it in the next. +\noindent and -A particle which finds itself in an inactive cell will terminate with status code 7. This is consistent with MODPATH 7's behavior. +\begin{equation} +v_{\tilde{y}} = \frac{d \tilde{y}}{dt} = v_{\tilde{y}}^{(0)} + w_{\tilde{y} \tilde{x}} \tilde{x} + + w_{\tilde{y} \tilde{y}} \tilde{y} , +\label{eqn:vylinear} +\end{equation} -A particle in a dry but active cell, or above the water table in a partially saturated cell, need not terminate. We call such a particle dry. The PRP package provides an option DRY\_TRACKING\_METHOD determining how dry particles should behave. Supported values are DROP (default), STOP, and STAY. +\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. -If DROP is selected, or if a DRY\_TRACKING\_METHOD is unspecified, a dry particle is passed vertically and instantaneously to the water table (if the cell is partially saturated) or to the bottom of the cell (if the cell is dry). This repeats (i.e. the particle may drop through multiple cells) until it reaches the water table. Tracking then proceeds as usual. If the vertical column containing the particle is entirely dry, the particle will terminate upon reaching the bottom of the model grid. +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. -If STOP is selected, dry particles will be terminated. +Because of its relatively simple form, the expression for $\beta$ can be inverted to yield an expression for $t$ as a function of $\beta$. Evaluating the subcell face that corresponds to $\beta = 0$ (the cell face) as a possible exit face for the particle is then simply a matter of calculating the exit time by substituting $\beta = 0$ into the expression for $t$. (A negative exit time indicates that that is not an exit face.) Evaluating the subcell face that corresponds to $\alpha = 0$ (one of the subcell faces internal to the cell) is essentially a one-dimensional root-finding exercise: for a given value of $\beta$, one evaluates $t$ and then $\alpha$, and iteration on $\beta$ continues until $\alpha = 0$ is solved. Evaluating the second subcell face internal to the cell as a possible exit face can be formulated as a similar root-finding problem. Because finding the root involves numerical iteration, the method is considered semi-analytical. -If STAY is selected, a dry particle will remain stationary until a) the water table rises and tracking can continue, b) the particle terminates due to reaching its STOPTIME or STOPTRAVELTIME, or c) the simulation ends. +\subsection{Packages} \label{sec:packages} -\subsection{Particle Track Output} +The Model Input (MIP) Package is used to specify properties of the model cells: porosity and, optionally, a retardation factor that modifies particle velocities and a zone number that can be used to define particle termination criteria. -The PRT Model supports both binary and CSV particle track output files. A particle track CSV file contains the output data in tabular format. The first line of the CSV file contains column names. Each subsequent line in the file contains a row of data for a single particle track record, with the following fields: +Particles are released into the model using the Particle Release Point (PRP) Package. Multiple PRP packages can be defined for a given PRT model. Each PRP package, which corresponds roughly to a ``particle group'' in MODPATH 7 \citep{pollock2016modpath7}, specifies tracking parameters and release locations and times for the particles it releases. -\vspace{5mm} -\noindent Column 0: \texttt{`KPER'} {\color{red} \footnotesize{INTEGER}} \\ -\noindent Column 1: \texttt{`KSTP'} {\color{red} \footnotesize{INTEGER}} \\ -\noindent Column 2: \texttt{`IMDL'} {\color{red} \footnotesize{INTEGER}} \\ -\noindent Column 3: \texttt{`IPRP'} {\color{red} \footnotesize{INTEGER}} \\ -\noindent Column 4: \texttt{`IRPT'} {\color{red} \footnotesize{INTEGER}} \\ -\noindent Column 5: \texttt{`ILAY'} {\color{red} \footnotesize{INTEGER}} \\ -\noindent Column 6: \texttt{`ICELL'} {\color{red} \footnotesize{INTEGER}} \\ -\noindent Column 7: \texttt{`IZONE'} {\color{red} \footnotesize{INTEGER}} \\ -\noindent Column 8: \texttt{`ISTATUS'} {\color{red} \footnotesize{INTEGER}} \\ -\noindent Column 9: \texttt{`IREASON'} {\color{red} \footnotesize{INTEGER}} \\ -\noindent Column 10: \texttt{`TRELEASE'} {\color{red} \footnotesize{DOUBLE}} \\ -\noindent Column 11: \texttt{`T'} {\color{red} \footnotesize{DOUBLE}} \\ -\noindent Column 12: \texttt{`X'} {\color{red} \footnotesize{DOUBLE}} \\ -\noindent Column 13: \texttt{`Y'} {\color{red} \footnotesize{DOUBLE}} \\ -\noindent Column 14: \texttt{`Z'} {\color{red} \footnotesize{DOUBLE}} \\ -\noindent Column 15: \texttt{`NAME'} {\color{red} \footnotesize{CHARACTER(LEN=LENBOUNDNAME)}} \\ +The Flow Model Interface (FMI) Package provides groundwater-flow information to the PRT Model. If a PRT model is run in the same simulation as the corresponding GWF model, the FMI Package provides updated flow information from the GWF model to the PRT model on each time step. If the PRT model is run is a separate simulation, the FMI Package reads flow information from the binary budget and head files written during the flow simulation. -\vspace{2mm} -\noindent where +\subsection{Exchanges} \label{sec:exchanges} -\begin{description} \itemsep0pt \parskip0pt \parsep0pt -\item \texttt{KPER} is the stress period number -\item \texttt{KSTP} is the time step number -\item \texttt{IMDL} is the number of the model the particle originated in -\item \texttt{IPRP} is the number of the particle release point (PRP) package the particle originated in -\item \texttt{IRPT} is the release point number -\item \texttt{ILAY} is the layer number -\item \texttt{ICELL} is the cell number -\item \texttt{IZONE} is the zone number -\item \texttt{ISTATUS} is the particle status code -\item \texttt{IREASON} is the reporting reason code -\item \texttt{TRELEASE} is the particle release time -\item \texttt{T} is the particle tracking time -\item \texttt{X} is the particle x coordinate -\item \texttt{Y} is the particle y coordinate -\item \texttt{Z} is the particle z coordinate -\item \texttt{NAME} is the name of the particle release point -\end{description} +The GWF-PRT Exchange establishes the link between a PRT model and its corresponding GWF model when the two models are run in the same simulation. As the version of the PRT Model described here does not support exchange of particles from one model to another, a PRT-PRT Exchange has not been implemented. -The ISTATUS field indicates the status of the particle: +\subsection{Backward Tracking} \label{sec:backwardtracking} -\begin{description} \itemsep0pt \parskip0pt \parsep0pt -\item \texttt{0}: particle was released -\item \texttt{1}: particle is being actively tracked -\item \texttt{2}: particle terminated at a boundary face -\item \texttt{3}: particle terminated in a weak sink cell -\item \texttt{4}: \textit{unused} -\item \texttt{5}: particle terminated in a cell with no exit face -\item \texttt{6}: particle terminated in a stop zone -\item \texttt{7}: particle terminated in an inactive cell -\item \texttt{8}: particle terminated immediately upon release into a dry cell -\item \texttt{9}: particle terminated in a subcell with no exit face -\end{description} +Unlike MODPATH 7, the PRT Model does not include an option for ``backward tracking,'' i.e., tracking of particles backward in time. However, backward tracking can be performed by running a GWF Model flow simulation, reversing the order of the flow budget and head information written to the binary budget and head files, and then running a PRT Model simulation based on the reversed budget and head files. To facilitate the process, FloPy \citep{bakker2016flopy, hughes2023flopy, flopysoftware}, a collection of python scripts that offers full support for preprocessing and postpostcessing \mf simulations, offers a utility to reverse binary budget and head files. -The IREASON field indicates the reason the particle track record was saved: +\subsection{Summary of Main Assumptions and Limitations} -\begin{description} \itemsep0pt \parskip0pt \parsep0pt -\item \texttt{0}: particle was released -\item \texttt{1}: particle exited a cell -\item \texttt{2}: time step ended -\item \texttt{3}: particle terminated -\item \texttt{4}: particle exited a weak sink cell -\item \texttt{5}: user-specified tracking time -\end{description} +The following is a list of the main assumptions and limitations of the PRT Model for \mf documented here: -By default, the PRT Model reports all particle events. The user can optionally select which events are reported, as explained in the Output Control (OC) subsection. Because multiple tracking events may coincide (e.g. exiting a cell and exiting a weak sink cell), particle track records may be duplicates except for the ISTATUS and/or IREASON codes. - -The binary particle track file contains the same particle track data in a record-based binary format explained in the Particle Track File subsection of the Description of Binary Output Files section. - - - - - -\newpage -\subsection{PRT Model Name File} -\input{prt/namefile.tex} - -%\newpage -%\subsection{Structured Discretization (DIS) Input File} -%\input{gwf/dis} - -%\newpage -%\subsection{Discretization with Vertices (DISV) Input File} -%\input{gwf/disv} - -%\newpage -%\subsection{Unstructured Discretization (DISU) Input File} -%\input{gwf/disu} - -\newpage -\subsection{Model Input (MIP) Package} -\input{prt/mip} - -\newpage -\subsection{Particle Release Point (PRP) Package} -\input{prt/prp} - -\newpage -\subsection{Output Control (OC) Option} -\input{prt/oc} - -\newpage -\subsection{Flow Model Interface (FMI) Package} -\input{prt/fmi} +\begin{itemize} +\item Particle tracks calculated by the PRT Model take into account advective motion only. Diffusive or dispersive transport of particles is not simulated. +\item The PRT Model uses the same spatial discretization used by the GWF Model, which may be represented using either a structured (DIS) or an unstructured (DISV) grid. Grids of DISU type are not supported. +\item Tracking of particles through advanced stress package features such as lakes and streams reaches and exchange of particles between PRT models is not supported. +\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{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 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}