This is a brief follow up to my last post on February 5, 2018 which described a python script for parsing a synthetic alluvial fan realization to support a MODFLOW simulation. With a few relatively minor additions, the script has now been extended to also support development of a reactive transport simulation using the U.S. Geological Survey’s PHAST model. Specifically, these additions include:

  • PHAST uses a simple x-y-z grid, as opposed to the layer model employed by MODFLOW, so the hydraulic conductivity values represented at each column, row, and layer in the MODFLOW input files are written to an x-y-z file which is read and interpolated by PHAST.
  • The ‘prism’ keyword for property zone definition within PHAST’s transport input file is employed to delineate the active and inactive cells across the model, utilizing (1) the alluvial fan surface as a top surface, and (2) the minimum fan thickness near the edges of the model to define the prism perimeter.
  • For reactive transport, the spatial distributions of mineral phase or reactive surface abundances can be correlated with hydraulic conductivity and binned (i.e., digitized) into a user-defined number of groups. These are read, in turn, as a separate input file and interpolated by PHAST.

The script, input files, and read-me file in my GitHub repository for this project have all been updated to reflect these additions.

The example synthetic alluvial fan developed in my December 27, 2018 post (which was parsed to enable a MODFLOW simulation, as described in my last post) was also used here to demonstrate a PHAST simulation. This demo scenario now entails the introduction of a user-specified groundwater composition at the two sediment input points (i.e., stream channels along left boundary). This composition that differs from that of the ambient groundwater present across the entire fan as an initial condition:


fan process summary
Visualization of flow and transport through synthetic alluvial fan (vertical exag. = 10X). Blue shading within alluvium indicates relatively permeable material. Backward particle tracks are color-coded to indicate relative groundwater velocity. Contour plot indicates the location of a conservative tracer (one horizontal slice) after 8 years of steady influx of water with differing composition at sediment source locations.

While the steady-state flow model component converged for the unconfined aquifer simulation that featured a number of dry cells (as was the case with the earlier MODFLOW run), the PHAST solute transport model developed numerical errors that I left un-diagnosed for now, owing to time constraints. As a workaround for this demonstration exercise, I instead ran the simulation assuming confined conditions across all active cells representing the fan. This eliminated the solute transport issue but produced a slightly different flow field than that generated by the earlier MODFLOW run.

Two groundwater chemistry scenarios were addressed: a conservative scenario with no water-rock interactions, and a reactive scenario with a cation exchanger present throughout the fan. The spatial distribution of the abundance of the cation exchanger, expressed as moles of exchange sites per kilogram of water in the pore space, was delineated into ten groups featuring a gradation from a relatively sparse to a relatively dense number of sites. These values were assumed to be inversely correlated with the local hydraulic conductivity (or, conceptually, correlated with the clay content). The cation exchanger was assumed to be in equilibrium with the ambient fan groundwater composition as an initial condition.

The PHAST transport and chemistry input files are also included in the project GitHub repository. Some sample results are provided here.

Conservative solute transport model: cation ratios resulting purely from mixing of end member waters.
Solute transport and mixing that also entails cation exchange: modeled cation ratios in groundwater after 8 years of elapsed simulation time.