In dpphat/MFtools: This is a package for reading MODFLOW and MT3DMS files

The MFtools library is used for investigating the input and output files from MODFLOW simulations. The library is based on fixed-width file formats as are generated when using Groundwater Vistas. This library has been tested using Groundwater Vistas files. The library includes functions for the following tasks:

Reading the discritization file (readdis(rootname))
Reading the layer property flow file (readlpf(rootname))
Reading the basic transport package file (readbtn(rootname))
Reading the headsave file (readhds(rootname, NLAY, NTTS))
Reading the unformatted concentration file (readucn(ucrootname, NLAY, NTTS))
Reading the well package file (readwel(rootname, NSP))
Reading no-flow array file (readnoflow(FILE, rootname))
Writing the headsave file (readhds(hds, ofl))
Writing the unformatted concentration file based on a concentration dataframe (writeucn(ucn, ofl))

The library also includes a set of utility functions for the following tasks:

Calculating contaminant mass based on information from the headsave file (i.e., saturated thickness), porosity read from the basic transport package, and concentrations from the unformatted concentration file (masscalc(BTN, HDS, UCN))
Reading a noflow matrix file (readnoflow(FILE, rootname))
Coordinate rotation for rotated models. This utility calculates X and Y coordinates from a data-frame (i.e., the BOT dataframe in the discretization file) based on a user-specified rotation and X and Y off-sets (rot(X_off,\ Y_off,\ ROT,\ Xin\ =\ X,\ Yin\ =\ Y))
Coordinate rotation from global coordinates to model coordinates (rot2mod(X_off,\ Y_off,\ ROT,\ Xin\ =\ X,\ Yin\ =\ Y))
Calculating longitudinal dispersivity based on the approach developed by Xu and Eckstein (1995) (disp(PLUME_L,\ UNITS\ =\ "feet"))

The functions from the MFtools library use rootnames to read files. Rootnames are based on the MODFLOW rootname convention such that files consist of a rootname and an extension that describes the MODFLOW file.

(NOTE that by my personal convention my script headings include the libraries that I will be using):

ROOT_NM <- system.file("extdata", "Ogll_sv2.01_ss.dis", package = "MFtools") %>% 
tools::file_path_sans_ext()

library(MFtools)
library(tidyverse)
library(magrittr)
library(scales)

rnm <- ROOT_NM

With the rootname we can now read in the discretization file

d <- readdis(rnm)
d

The structure of the variable, d, is a list consisting of the information from the MODFLOW discretization file. For indexing purposes, some of the information is arranged in a dataframe, such as the TOP and BOT variables that describe the top and bottom elevations of each model layer.

MFtools includes a function to read the layer property flow file (.lpf).

p <- readlpf(rnm)
p

The .lpf file has lots of really useful information. We can isolate this information using the "$" indexing operation.

p$PROPS

p$PROPS %>% group_by(LAY) %>% summarise(MIN = min(HK), MEDIAN = median(HK), MAX = max(HK))

The model is a four layer model. Summarising does not give us a lot of information about the distribution of hydraulic conductivity values.

We can, instead, use ggplot to plot a histogram of hydraulic conductivity values by layer

p$PROPS %>% ggplot(aes(x = HK)) +
            geom_histogram(aes(y = (..count..)/tapply(..count.., ..PANEL.., sum)[..PANEL..]),
                colour = "black", 
                fill = "yellow", 
                alpha = 0.6) +
            scale_y_continuous(limits = c(0, 0.2), 
                breaks = seq(0, 0.2, 0.05), 
                labels = percent, 
                name = "FREQUENCY COUNT (PERCENT)") +   
            scale_x_log10(breaks = 10^(-10:10), 
                 minor_breaks = c(1:9 %o% 10^(-10:10)), 
                 labels = scales::trans_format("log10", math_format(10^.x))) +
            facet_wrap(~LAY, ncol = 1) +
            labs(x = "HYDRAULIC CONDUCTIVITY (FEET / DAY)") +
            theme_bw()

We can also investigate groundwater elevations by model layer

h <- readhds(rnm, NLAY = d$NLAY, 1)
h %>% filter(GWE > 2000) %>% filter(GWE < 10000) %>%
            ggplot(aes(x = GWE)) +
            geom_histogram(aes(y = (..count..)/tapply(..count.., ..PANEL.., sum)[..PANEL..]),
                colour = "black", 
                fill = "blue", 
                alpha = 0.6) +
            scale_y_continuous(limits = c(0, 0.2), 
                breaks = seq(0, 0.2, 0.05), 
                labels = percent, 
                name = "FREQUENCY COUNT (PERCENT)") +   
            facet_wrap(~LAY, ncol = 1) +
            labs(x = "GROUNDWATER ELEVATION (FEET AMSL)") +
            theme_bw()