Scenario

Level 1 data

These are processing instructions for creating level 1 modified datasets.

eval(parse(text = app.c.chunks[["comb_sw_irr_1"]]))

eval(parse(text = app.c.chunks[["pod_gw_1"]]))

eval(parse(text = unlist(app.c.chunks[c("pod_wells_1", "pod_wells_2")])))

eval(parse(text = app.c.chunks[["irr_lands_1"]]))

eval(parse(text = unlist(app.c.chunks[c("et_1", "et_2")])))

Level 2 data

These are processing instructions for creating level 2 modified datasets. Recall that level 2 datasets are dependent on level 1 data.

eval(parse(text = unlist(app.c.chunks[paste0("entity_components_", 1:3)])))

eval(parse(text = app.c.chunks[["rs_entities_1"]]))

eval(parse(text = app.c.chunks[["rs_rech_non_irr_1"]]))

Copy calibrated files

Copy all the calibrated model files from the “archive” folder to a new folder named “model” in our working directory.

dir.create("model", showWarnings = FALSE)
files <- list.files("archive", full.names = TRUE)
file.copy(files, "model", overwrite = TRUE)

Specify the identifier for the model run.

id <- "wrv_mfusg"

Load R objects from the “model/model.rda” file into your R environment, these objects include d.in.mv.ave, misc, reduction, rs, ss.stress.periods, tr.stress.periods, and trib. A description of each R object is given in the help documentation of the UpdateWaterBudget function in the wrv package.

load("model/model.rda")

The rs object in the appendix D code is named rs.model; use the following command to address this issue:

rs.model <- rs

Recalculate Specified Flows

Recalculate the specified-flow boundary conditions using our modified datasets—this results in a new MODFLOW Well file (“model/wrv_mfusg.wel”) and quality assurance tables (“model/qa-*.csv”). Note that all modified datasets need to be passed as arguments in the UpdateWaterBudget function.

UpdateWaterBudget("model", id, qa.tables = "english", pod.wells = pod.wells,
                  comb.sw.irr = comb.sw.irr, et = et,  pod.gw = pod.gw,
                  entity.components = entity.components, rs.entities = rs.entities)

Run MODFLOW

Run MODFLOW-USG using model scenario conditions.

wd <- setwd("model")
system2(file.path(getwd(), "RunModflow.bat"), stdout = FALSE, stderr = FALSE)
setwd(wd)

Read water budgets

Read water-budget output for the archived calibrated simulation, located in the “archive” folder, and save the results for later use. To do so, run the code chunks named “read_budget_1” and “read_budget_2” in appendix D using the following commands:

dir.run <- "archive"  # code below is dependent on this object
eval(parse(text = unlist(app.d.chunks[paste0("read_budget_", 1:2)])))

And save the average volumetric flow rates to an R object named budget1 using

budget1 <- c("Water-table recharge"             = mean(d.rech$flow.in),
             "Streamflow losses"                = mean(d.river$flow.in),
             "Tributary basin underflow"        = mean(d.trib$flow),
             "Water-table discharge"            = mean(d.rech$flow.out),
             "Streamflow gains"                 = mean(d.river$flow.out),
             "Production well pumping"          = mean(d.well$flow),
             "Stanton Crossing outlet boundary" = mean(d.drain.1$flow),
             "Silver Creek outlet boundary"     = mean(d.drain.2$flow))
budget1 <- setNames(as.integer(abs(budget1)) * 0.296106669, names(budget1))

Next, read the water-budget output for the scenario simulation, located in the “model” folder.

dir.run <- "model"
eval(parse(text = unlist(app.d.chunks[paste0("read_budget_", 1:2)])))

And save the average volumetric flow rates to an R object named budget2 using

budget2 <- c("Water-table recharge"             = mean(d.rech$flow.in),
             "Streamflow losses"                = mean(d.river$flow.in),
             "Tributary basin underflow"        = mean(d.trib$flow),
             "Water-table discharge"            = mean(d.rech$flow.out),
             "Streamflow gains"                 = mean(d.river$flow.out),
             "Production well pumping"          = mean(d.well$flow),
             "Stanton Crossing outlet boundary" = mean(d.drain.1$flow),
             "Silver Creek outlet boundary"     = mean(d.drain.2$flow))
budget2 <- setNames(as.integer(abs(budget2)) * 0.296106669, names(budget2))

Compare water budgets

A comparison between water budgets is made by constructing the following table.

Table: Water budgets for modified and archived simulation, specified as volumetric flow rates average over 1998–2010. [Inflow: water entering the aquifer system. Outflow: water leaving the aquifer system. Component: a water budget component in the groundwater-flow model. Archived Rate: is the mean volumetric flow rate for the archived simulation. Modified Rate: is the mean volumetric flow rate for the modified simulation. Difference: is the archived rate subtracted from the modified rate. Abbreviations: acre-ft/yr, acre-feet per year]

d <- data.frame(rate1 = budget1, rate2 = budget2)
d$direction <- c("**Inflow**", "", "", "**Outflow**", "", "", "", "")
d$component <- names(budget1)
d$diff <- budget2 - budget1
d$pchange <- ((budget2 - budget1) / budget1) * 100
d <- rbind(d[1:3, ], NA, d[4:8, ], NA, NA)
d$component[nrow(d)] <- "**Inflow - Outflow**"
d$rate1[nrow(d)] <- sum(budget1[1:3]) - sum(budget1[5:8])
d$rate2[nrow(d)] <- sum(budget2[1:3]) - sum(budget2[5:8])
rownames(d) <- NULL
d <- d[, c("direction", "component", "rate1", "rate2", "diff", "pchange")]
colnames(d) <- c("", "Component", "Archived rate<br>(acre-ft/yr)",
                 "Scenario rate<br>(acre-ft/yr)", "Difference<br>(arce-ft/yr)",
                 "Percent<br>change")
d[, 3] <- formatC(d[, 3], format = "f", digits = 0, big.mark = ",")
d[, 4] <- formatC(d[, 4], format = "f", digits = 0, big.mark = ",")
d[, 5] <- formatC(d[, 5], format = "f", digits = 0, big.mark = ",")
d[, 6] <- formatC(d[, 6], format = "f", digits = 1, big.mark = ",")
d[, 1:2][is.na(d[, 1:2])]  <- ""
d[, 3:6][d[, 3:6] == "NA"] <- ""
knitr::kable(d, format = "markdown", padding = 0, booktabs = TRUE, digits = 0,
             row.names = FALSE, col.names = colnames(d), align = "llrrrr")

Read hydraulic heads

Read simulated hydraulic heads for archived simulation, located in “archive” folder, and save results.

Read the simulated hydraulic head (head) output for the archived calibrated model, located in the “archive” folder, and save the results for later use. To do so, run the code chunk named “read_head” in appendix D using the following commands:

dir.run <- "archive"
eval(parse(text = app.d.chunks[["read_head"]]))
head1 <- rs.heads.lay1

Next, read the head output for the scenario simulation, located in the “model” folder.

dir.run <- "model"
eval(parse(text = app.d.chunks[["read_head"]]))
head2 <- rs.heads.lay1

Compare head distributions

Draw a map the simulated water table difference, defined as the archived water table subtracted from the modified water table, during August 16, 2007. Add an observation point located at a longitude and latitude of -114.222664 and 43.383047 decimal degrees, respectively.

ext <- spTransform(alluvium.extent, CRS("+init=epsg:4326"))
r <- head2[["2007-08-16"]] - head1[["2007-08-16"]]
r[] <- round(r[] * 3.28084, digits = 3)
r <- projectRaster(r, crs = CRS("+init=epsg:4326"), method = "ngb")
Pal <- colorNumeric("Spectral", r[], na.color = "transparent")
map <- CreateWebMap("Topo")
ll  <- coordinates(ext)
map <- setView(map, lng = ll[1], lat = ll[2], zoom = 11)
map <- addRasterImage(map, r, colors = Pal, opacity = 0.8)
map <- addPolylines(map, data = ext, weight = 3, color = "#000000")
map <- addLegend(map, pal = Pal, values = r[], opacity = 0.8,
                 title = "Head Diff.<br/>(feet)")
ll <- c(-114.222664, 43.383047)
txt <- sprintf("<b>Longitude:</b> %s<br/><b>Latitude:</b> %s", ll[1], ll[2])
icon <- icons(iconUrl = "markers/marker-blue.png",
              iconWidth = 34, iconHeight = 34, iconAnchorX = 17, iconAnchorY = 34)
map <- addMarkers(map, lng = ll[1], lat = ll[2], popup = txt, icon = icon)
map

Finally, at the observation point, plot the head difference in model layer 1 over time.

pnt <- data.frame(lon = ll[1], lat = ll[2])
coordinates(pnt) <- c("lon", "lat")
proj4string(pnt) <- CRS("+init=epsg:4326")
pnt <- spTransform(pnt, crs(hill.shading))
d <- t(extract(head2, pnt) - extract(head1, pnt))
d <- data.frame(Date = as.Date(rownames(d)), difference = d)
ylab <- paste("Hydraulic head difference, in", c("meters", "feet"))
PlotGraph(d, ylab = ylab, col = "#025D8C", conversion.factor = 3.28084,
          center.date.labels = TRUE, seq.date.by = "year")