In gopalpenny/anem: Analytical element models for simple aquifers
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
fig.width=4,
fig.height=2.5
)
library(ggplot2)
library(dplyr)
library(anem)
library(dplyr)
This package was created to evaluate the hydraulic relationships among wells, in order to estimate the effect of a group of wells on drawdown at these wells or other wells. It it based on method of images from the analytical element modeling approach, in which the 2-dimensional characteristics of aquifers are reproduced by strategically placing wells within the domain.
This package models simple aquifer and well configurations. The boundaries of the aquifer can be specified as no flow or constant head boundaries, and the corners of the aquifer must be right angles. For fully bounded aquifers, this means that the aquifer must be a rectangle. The constant head boundaries take the head of the undisturbed aquifer, h0.
Units are designed for length dimensions given in meters (or m^2^), and time in seconds.
As you will see in this vignette, the package has a number of functions that are used to set up the aquifer boundaries and wells. For practical purposes, the package was designed for simple but real aquifers, and wells and boundaries of the aquifers can be imported as shapefiles and then prepared. The basic functionality of the package, described in more detail in sections below, can be seen in the following figure.
bounds_df <- define_bounds(tibble(bound_type=c("PB","PB","CH","NF"),m=c(Inf,0,Inf,0),b=c(0,5000,5000,0)))
recharge_params <- list(recharge_type="F", recharge_vector=c(0,0,1,0),
flow=5e-6, x0=5000, y0=5000)
aquifer <- define_aquifer("confined", 1e-3, h0=0, z0=10,
bounds=bounds_df,recharge=recharge_params)
wells_orig <- define_wells(x=c(1000,1700,3500,4000),y=c(4000,1700,3500,1000),
R=rep(4000,4),diam=rep(1,4),Q=c(-0.02,-0.02,0.02,-0.02))
wells <- wells_orig %>%
anem::generate_image_wells(aquifer)
# wells <- define_wells(x=2500,y=2500,
# R=6000,diam=1,Q=-0.01) %>%
# anem::generate_image_wells(aquifer)
# anem::get_gridded_hydrodynamics()
loc <- expand.grid(x=seq(0,10000,length.out=201),y=seq(-5000,5000,length.out=201))
thm <- theme(panel.background = element_rect(fill=NA,color=NA),
panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
legend.position="none",axis.text=element_blank(),
axis.ticks=element_blank(),axis.title=element_blank(),
legend.title=element_text(size=9),
legend.text=element_text(size=8),
legend.spacing = unit(1,"mm"),
legend.margin = margin(0.25,0,0,0.25,unit = "mm"),
legend.box.margin = margin(0,0,0,0,unit="mm"))
p_base <- ggplot() +
geom_segment(data=aquifer$bounds %>%
mutate(bt=factor(bound_type,levels=c("CH","NF","PB"),
labels=c("River","No flow","Open"))) %>%
filter(bt !="Open"),
aes(x1,y1,xend=x2,yend=y2,color=bt),size=1) +
scale_shape_manual("Well type",values = c(16,1,4,3)) +
scale_color_manual("Boundary",values = c("blue","#555555","#CCCCCC"))
well_color <- "#AA3333"
constant_flow <- tibble(y=seq(500,5000-500,length.out = 6),x=0)
f1_a <- p_base +
geom_segment(data=constant_flow,aes(x,y,xend=x+500,yend=y),arrow = arrow(length=unit(1,"mm"),type = "closed")) +
geom_point(data=wells %>% filter(well_image=="Actual"),
aes(x,y,shape=well_type),size=3,stroke=1) +
# geom_point(data=wells %>% filter(well_image=="Actual",orig_wID==4),
# aes(x,y,shape=well_type),size=3,stroke=1,color=well_color,show.legend = FALSE) +
coord_equal() +
labs(subtitle = "(a)") +
thm %+replace% theme(legend.position = "left")
wells_roi <- wells %>% filter(orig_wID==4) %>% anem::gen_circles()
f1_b <- p_base +
geom_polygon(data=wells_roi,aes(x,y,group=id),fill=well_color,alpha=0.2,color=NA,size=0.25) +
geom_point(data=wells %>% filter(well_image=="Actual",orig_wID!=4),
aes(x,y,shape=well_type),size=3,stroke=1,alpha=0.5) +
# geom_point(data=wells %>% filter(well_image=="Actual",orig_wID==4),
# aes(x,y,shape=well_type),size=3,stroke=1,color=well_color) +
geom_point(data=wells %>% filter(well_image=="Actual",orig_wID==4),aes(x,y),
size=3,stroke=1,color=well_color,shape=16) +
geom_point(data=wells %>% filter(grepl("\\+Q",well_image),orig_wID==4),aes(x,y),
size=2.5,stroke=1,color=well_color,shape=15) +
geom_point(data=wells %>% filter(grepl("\\-Q",well_image),orig_wID==4),aes(x,y),
size=2,stroke=1,color=well_color,shape=0) +
annotate("rect",xmin=3500,xmax=6500,ymin=-1500,ymax=1500,color="#992222",linetype="dashed",fill=NA,size=1) +
annotate("point",500,-4500,shape=0,size=2,color=well_color,stroke=1) +
annotate("text",1000,-4500,shape=0,size=4.5,label="Image wells",
color="black",hjust=0,vjust=0.45) +
coord_equal() +
labs(subtitle = "(b)") +
thm #%+replace% theme(legend.position="right")
streamfunction <- loc %>%
dplyr::bind_cols(streamfunction=get_stream_function(loc,wells %>% filter(orig_wID==4),aquifer)) %>%
dplyr::bind_cols(head=get_hydraulic_head(loc,wells,aquifer))
flowlines <- anem::get_contourlines(streamfunction %>% rename(z=streamfunction),nlevels=40)
f1_c <- p_base +
geom_path(data=flowlines,aes(x,y,group=line),color="#66CCEE",alpha=1) +
geom_point(data=wells %>% filter(well_image=="Actual"),aes(x,y,shape=well_image),
size=2,stroke=1,color="#AA3333") +
geom_point(data=wells %>% filter(well_image!="Actual"),aes(x,y,shape=well_image),
size=2,stroke=1,color="#AA3333") +
annotate("rect",xmin=3500,xmax=6500,ymin=-1500,ymax=1500,color="#AA3333",linetype="dashed",fill=NA,size=1) +
coord_equal(xlim=c(3500,6500),ylim=c(-1500,1500)) +
scale_shape_manual(values=c(16,0,15))+
labs(subtitle = "(c)") +
thm #%+replace% theme(legend.position = "right")
wells_actual <- wells %>% filter(well_image=="Actual")
gridded <- get_gridded_hydrodynamics(wells,aquifer,head_dim=c(200,200),flow_dim=c(10,10))
missing_points <-
do.call(rbind,lapply(split(wells_actual,seq(nrow(wells_actual))),
function(loc,df) df[which.min(sqrt((df$x - loc$x)^2 +
(df$y - loc$y)^2))[1],c("x","y")], df=gridded$flow))
x_multiplier <- 2e6
gridded$flow <- gridded$flow %>%
mutate(mag=sqrt(dx^2 + dy^2),
multiplier = (mag)^(1/4)*x_multiplier,
x2 = x + dx * multiplier,
y2 = y + dy * multiplier) %>%
anti_join(missing_points,by=c("x","y"))
# contours <- get_contourlines(gridded$head %>% rename(z=head_m),nlevels=20)
with(gridded$flow, max(abs(x2-x)))
f1_d <- ggplot() +
geom_raster(data=gridded$head,aes(x,y,fill=head_m)) +
# geom_path(data=contours,aes(x,y,group=line)) +
geom_segment(data=aquifer$bounds %>%
mutate(bt=factor(bound_type,levels=c("CH","NF","PB"),
labels=c("River","No flow","Open (no boundary)"))) %>%
filter(bound_type!="PB"),
aes(x1,y1,xend=x2,yend=y2,color=bt),size=2,show.legend = FALSE) +
geom_segment(data=gridded$flow,aes(x,y,xend=x2,yend=y2),
arrow = arrow(ends="last",type="closed",length=unit(0.75,"mm")),color="black") +
geom_point(data=wells %>% filter(well_image=="Actual"),
aes(x,y,shape=well_type),size=3,stroke=1,show.legend = FALSE) +
scale_shape_manual("Well type",values = c(16,1,4,3)) +
scale_color_manual("Boundary",values = c("blue","#555555","#CCCCCC")) +
scale_fill_viridis_c("Hydraulic\nhead, m",breaks=-1:2,labels=c("-1","0 (river)","1","2")) +
labs(subtitle = "(d)") +
coord_equal() + thm %+replace% theme(legend.position = "right",legend.key.height=unit(0.55,"cm"))
cap <- "(a) Aquifer scenario with four wells and arrows indicating constant background flow towards the river. (b) Well images and radii of influence for bottom-right well (red). (c) Reproduction of the no-flow and constant-head boundaries for single well using well images. (d) Hydraulic head and flow field for scenario. Note that the stream, flowing from top to bottom, goes from gaining to losing."
# ggpubr::ggarrange(f1_a,f1_b,f1_c,f1_d,nrow=1,labels=c("(a)","(b)","(c)","(d)"),widths=c(0.32,0.2,0.2,0.32)) %>%
# ggpubr::annotate_figure(bottom = "This is the caption.")
# f1_a | f1_b | f1_c | f1_d
gridExtra::grid.arrange(f1_a , f1_b , f1_c , f1_d,nrow=1,widths=c(0.33,0.2,0.2,0.35))
Understanding the method of images
The package relies on the principle of superposition and the method of images to generate groundwater hydrodynamics, illustrated by the two figures below (see ?get_segments_behavior and ?get_stream_function, respectively, which give examples for generating similar plots). The principle of superposition states that, because the groundwater equations are linear, the solution to the equations can be calculated as the sum of solutions. In the example below, the drawdown of the water table is calculated as the sum of drawdown due to three individual wells.
# define the aquifer
aquifer_super <- aquifer_confined_example
aquifer_super$Ksat <- 5e-4
# define wells and images
wells_super_actual <- define_wells(data.frame(x=c(400,500,900),y=500,R=1500,Q=-0.1,diam=1))
wells_super <- generate_image_wells(wells_super_actual,aquifer_super)
# set the segment
seg <- c(0,500,1000,500)
# get segment for each well and all wells
seg_w1 <- get_segments_behavior(seg,wells_super[wells_super$orig_wID==1,],aquifer_super,1000)
seg_w1 <- data.frame(seg_w1[,c("x","y","head")],Wells="W1")
seg_w2 <- get_segments_behavior(seg,wells_super[wells_super$orig_wID==2,],aquifer_super,1000)
seg_w2 <- data.frame(seg_w2[,c("x","y","head")],Wells="W2")
seg_w3 <- get_segments_behavior(seg,wells_super[wells_super$orig_wID==3,],aquifer_super,1000)
seg_w3 <- data.frame(seg_w3[,c("x","y","head")],Wells="W3")
seg_all <- get_segments_behavior(seg,wells_super,aquifer_super,1000)
seg_all <- data.frame(seg_all[,c("x","y","head")],Wells="All")
# plot results
seg_behavior <- rbind(seg_all,seg_w1,seg_w2,seg_w3)
ggplot() +
geom_line(data=seg_behavior,aes(x,head,linetype=Wells)) +
geom_vline(data=data.frame(x=c(0,1000)),aes(xintercept=x),color="gray") +
annotate(geom="text",x=40,y=86,label="Constant head boundary",angle=-90) +
annotate(geom="text",x=1040,y=86,label="No flow boundary",angle=-90)
In the method of images, wells are mirrored across aquifer boundaries to reproduce characterics of the boundary, such as no flow or constant head. The following examples demonstrate (a) how a no-flow boundary is created by mirroring a well across the boundary, and (b) how a constant head boundary is recreated by mirroring the well and inverting the pumping.
# Create a grid of locations
loc <- expand.grid(x=seq(-200,200,length.out=201),y=seq(-200,200,length.out=201))
# Constant head boundary
wells_constant_head <- define_wells(x=c(-100,100),y=c(-0,0),Q=c(1e-2,-1e-2),diam=c(0.1,0.1),R=c(500,500))
aquifer <- define_aquifer("confined",1e-4,z0=20,h0=0,n=0.35)
constant_head_boundary <- loc %>%
bind_cols(streamfunction=get_stream_function(loc,wells_constant_head,aquifer)) %>%
bind_cols(head=get_hydraulic_head(loc,wells_constant_head,aquifer))
# No flow boundary
wells_no_flow <- define_wells(x=c(-100,100),y=c(-0,0),Q=c(-1e-2,-1e-2),diam=c(0.1,0.1),R=c(500,500))
no_flow_boundary <- loc %>%
bind_cols(streamfunction=get_stream_function(loc,wells_no_flow,aquifer)) %>%
bind_cols(head=get_hydraulic_head(loc,wells_no_flow,aquifer))
thm <- theme(axis.text=element_blank(),axis.title=element_blank(),axis.ticks=element_blank(),
panel.background = element_rect(fill=NA,color="black"),
panel.border = element_rect(fill=NA,color="black"),
plot.title = element_text(size=10))
p_CH <- ggplot() +
geom_contour(data=constant_head_boundary,aes(x,y,z=head),bins=50,linetype="dashed",color="black") +
geom_contour(data=constant_head_boundary,aes(x,y,z=streamfunction),bins=25,color="darkblue") +
geom_point(data=wells_constant_head,aes(x,y,fill=well_type),size=3,shape=21) +
theme(legend.position=c(0.8,0.15),legend.title=element_blank()) +
coord_equal() + thm + ggtitle("Constant head boundary at x = 0")
p_NF <- ggplot() +
geom_contour(data=no_flow_boundary,aes(x,y,z=head),bins=25,linetype="dashed",color="black") +
geom_contour(data=no_flow_boundary,aes(x,y,z=streamfunction),bins=50,color="darkblue") +
geom_point(data=wells_no_flow,aes(x,y,fill=well_type),size=3,shape=21) +
theme(legend.position=c(0.8,0.1),legend.title=element_blank()) +
coord_equal() + thm + ggtitle("No flow boundary at x = 0")
gridExtra::grid.arrange(p_NF,p_CH,nrow=1)
A simple example of drawdown relationships
Load packages for the vignette
library(ggplot2)
library(dplyr)
library(anem)
Aquifers in this package are characterized by:
- Aquifer type, confined or unconfined
- Saturated hydraulic conductivity, Ksat
- Resting hydraulic head, h0
- Thickness of the aquifer, z0 (confined aquifers only)
- Porosity of the aquifer, n
- Aquifer boundaries, characterized by:
- Boundary ID, or bID
- Coordinates in a cartesian plane, given as slope, m, and intercept, b, and coordinates (x1, y1, x2, y2)
- Boundary type, NF, CH, or PB
- Additional user-defined parameters
Use the define_aquifer() function to create a new aquifer, which has the class aquifer. This class is essentially a list with names parameters and a separate print method.
# define aquifer
bounds_df <- data.frame(bound_type=c("CH","NF","NF","NF"),m=c(Inf,0,Inf,0),b=c(0,1000,1000,0))
aquifer_unconfined <- define_aquifer("unconfined",1e-3,bounds=bounds_df,h0=100,n=0.35)
print(aquifer_unconfined)
Wells are characterized by:
- Well ID, or wID
- Location x, y
- Pumping rate (+ for injection, - for pumping), Q
- Radius of influence, R
- Diameter, diam
- Well images (which are noted in the
well_image column)
- Additional user-defined columns, including groups or weights.
The code below Defines 8 pumping wells using random locations and arbitrarily divides wells into countries "A" and "B" using a threshold at y = 500 m. The define_wells() function ensures that the wells have all appropriate columns, and the class of the returned object is a tibble (which functions like a data.frame, but has a couple bells and whistles). Also note that the radius of influence (R) is defined arbitrarily, but there is also a function to calculate this manually -- see ?get_ROI for more details.
The generate_image_wells() function generates well images to recreate the bounderies defined by the aquifer, aquifer_unconfined.
# define wells and well images
set.seed(15)
wells_actual <- define_wells(x = runif(8,0,1000),
y = c(runif(4,0,500),runif(4,500,1000)),
Q = -1/4,
diam = 1,
R = 1000,
weights = 1,
country = rep(c("A","B"),each=4))
wells <- wells_actual %>% generate_image_wells(aquifer_unconfined)
print(wells)
Plot the aquifer and wells.
ggplot() +
geom_segment(data=aquifer_unconfined$bounds,aes(x1,y1,xend=x2,yend=y2,linetype=bound_type)) +
geom_abline(slope=0,intercept=500,linetype="dashed") +
geom_point(data=wells_actual,aes(x,y,fill=country),shape=21) +
coord_equal()
Get drawdown relationships using get_drawdown_relationships(). This function calculates the average drawdown at wells in each group defined by the column group_column. The average is taken as the weighted mean, determined by the weights_column. The weights were previously set equal for all wells so that the result here is a simple mean. The results show PHIii and PHIij.
drawdown_relationships <- get_drawdown_relationships(wells, aquifer_unconfined, group_column=country, weights_column=weights)
drawdown_relationships
drawdown_relationships %>%
knitr::kable("html") #%>% kableExtra::kable_styling()
Plot the head and flow
The hydrodynamics of the aquifer can be mapped by obtaining gridded heand flow using the get_gridded_hydrodynamics() function. The function takes as input the wells, aquifer, and grid dimensions for head and flow. It returns a list object with data.frames for head and flow, which can then be plotted.
gridded_1 <- get_gridded_hydrodynamics(wells,aquifer_unconfined,c(20,20),c(8,8))
ggplot() +
geom_raster(data=gridded_1$head,aes(x,y,fill=head_m)) +
geom_contour(data=gridded_1$head,aes(x,y,z=head_m),color="white",alpha=0.3) +
geom_segment(data=gridded_1$flow,aes(x,y,xend=x2,yend=y2),
arrow = arrow(ends="last",type="closed",length=unit(1,"mm")),color="black") +
geom_segment(data=aquifer_unconfined$bounds,aes(x1,y1,xend=x2,yend=y2,linetype=bound_type)) +
geom_point(data=wells %>% filter(wID==orig_wID),aes(x,y),shape=21) +
coord_equal()
Check the head and flow along the boundaries
The function get_bounds_behavior() is a helper function to generate hydraulic properties at the boundaries. It obtains hydraulic head and the flow normal to the boundaries by setting boundaries to segments and using get_segments_behavior(). Normal flow is defined such that positive flow has some component in the x-direction (the y-direction depends on the normal vector for the boundary). The function plot_bounds_behavior() is a wrapper around get_bounds_behavior(), and it can be used to quickly compare the hydraulic head along the boundaries and flow across the boundaries under two scenarios:
- no images (ie, homogeneous, uniform aquifer -- all well images removed)
- images (ie, aquifer with boundaries -- all well images intact)
bounds_behavior <- plot_bounds_behavior(wells,aquifer_unconfined)
gridExtra::grid.arrange(bounds_behavior$p_h,bounds_behavior$p_f,nrow=1)
Finally, the function plot_bounds_behavior also includes the raw data of head and flow used to create the above plots (bounds_behavior$bounds_behavior), as well as a summary of these values which includes mean head on each of the boundaries and mean flow as the mean of the absolute value of flow normal to each boundary (bounds_behavior$table). To numerically check that the boundaries are working as expected, we can print bounds_behavior$table:
bounds_behavior$table %>%
knitr::kable("html") #%>% kableExtra::kable_styling()
Defining recharge / background flow
Background flow is parameterized by defining "recharge" in the aquifer. Recharge can be defined as a constant flow (F) or a recharge divide (D). In both cases, flow is defined as m^3^ / m-s, in the direction of the recharge vector. Recharge is therefore defined by specifying:
- Type of recharge: constant flow (F) or recharge divide (D)
- Recharge vector, pointing in the direction of flow
- Flow magnitude, in units of m^3^ / m-s
- Coordinates where the undisturbed water table (no wells) has a hydraulic head equal to aquifer$h0.
See ?define_recharge for more details.
# Define recharge for type "D"
recharge_params <- list(recharge_type="D",recharge_vector=c(500,500,501,501),
flow_main=1e-3,flow_opp=2e-3,x0=0,y0=0)
# Define bounds as "PB" -- pervious bounds, which have no effect (other than to define output grid)
bounds_recharge <- define_bounds(data.frame(bound_type=rep("PB",4),m=c(Inf,0,Inf,0),b=c(0,0,1000,1000)))
aquifer_recharge <- define_aquifer("confined",1,h0=0,z0=1,n=0.35,recharge=recharge_params,bounds=bounds_recharge)
gridded_2 <- get_gridded_hydrodynamics(wells_actual,aquifer_recharge,c(20,20),c(8,8))
The results can then be plotted, just as above:
ggplot() +
geom_raster(data=gridded_2$head,aes(x,y,fill=head_m)) +
geom_contour(data=gridded_2$head,aes(x,y,z=head_m),color="white",alpha=0.3) +
geom_segment(data=gridded_2$flow,aes(x,y,xend=x2,yend=y2),
arrow = arrow(ends="last",type="closed",length=unit(1,"mm")),color="black") +
geom_segment(data=aquifer_recharge$bounds,aes(x1,y1,xend=x2,yend=y2,linetype=bound_type)) +
geom_point(data=wells_actual,aes(x,y),shape=21) +
coord_equal()
Particle tracking
The package implements particle tracking in two ways: (1) tracking individual particles using track_particles and (2) estimating capture zones of wells using get_capture_zone, which initializes particles in the vicinity of wells and tracks them in reverse. Here we use get_capture_zone on the previous aquifer to generate particle tracking away from wells.
# wells <- wells
# wells[4,"Q"] <- 0.25
well_particles <- get_capture_zone(wells,aquifer_unconfined,t_max=365,n_particles=4)
gridded_3 <- get_gridded_hydrodynamics(wells,aquifer_unconfined)
ggplot() +
geom_raster(data=gridded_1$head,aes(x,y,fill=head_m)) +
geom_contour(data=gridded_1$head,aes(x,y,z=head_m),color="white",alpha=0.3) +
geom_segment(data=aquifer_unconfined$bounds,aes(x1,y1,xend=x2,yend=y2,linetype=bound_type)) +
geom_path(data=well_particles,aes(x,y,group=i,color=as.factor(wID)),show.legend = FALSE) +
geom_point(data=wells_actual,aes(x,y,color=as.factor(wID)),shape=21,show.legend = FALSE) +
coord_equal()
Integrating with anem-app
The basic functionality of this package is implemented in a Shiny application, which can be run by calling anem_app() or accessed online at shinyapps.io. The Shiny application provides an GUI for setting aquifer properties, drawing boundaries, defining recharge, adding wells, and tracking particles. It also displays results directly in the app. You can also download scenarios created on anem-app and import them into R. To do so, check out ?import_app_rds from this package to learn more.
gopalpenny/anem documentation built on Dec. 20, 2020, 5:27 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width=4, fig.height=2.5 ) library(ggplot2) library(dplyr) library(anem) library(dplyr)
This package was created to evaluate the hydraulic relationships among wells, in order to estimate the effect of a group of wells on drawdown at these wells or other wells. It it based on method of images from the analytical element modeling approach, in which the 2-dimensional characteristics of aquifers are reproduced by strategically placing wells within the domain.
This package models simple aquifer and well configurations. The boundaries of the aquifer can be specified as no flow or constant head boundaries, and the corners of the aquifer must be right angles. For fully bounded aquifers, this means that the aquifer must be a rectangle. The constant head boundaries take the head of the undisturbed aquifer, h0.
Units are designed for length dimensions given in meters (or m^2^), and time in seconds.
As you will see in this vignette, the package has a number of functions that are used to set up the aquifer boundaries and wells. For practical purposes, the package was designed for simple but real aquifers, and wells and boundaries of the aquifers can be imported as shapefiles and then prepared. The basic functionality of the package, described in more detail in sections below, can be seen in the following figure.
bounds_df <- define_bounds(tibble(bound_type=c("PB","PB","CH","NF"),m=c(Inf,0,Inf,0),b=c(0,5000,5000,0))) recharge_params <- list(recharge_type="F", recharge_vector=c(0,0,1,0), flow=5e-6, x0=5000, y0=5000) aquifer <- define_aquifer("confined", 1e-3, h0=0, z0=10, bounds=bounds_df,recharge=recharge_params) wells_orig <- define_wells(x=c(1000,1700,3500,4000),y=c(4000,1700,3500,1000), R=rep(4000,4),diam=rep(1,4),Q=c(-0.02,-0.02,0.02,-0.02)) wells <- wells_orig %>% anem::generate_image_wells(aquifer) # wells <- define_wells(x=2500,y=2500, # R=6000,diam=1,Q=-0.01) %>% # anem::generate_image_wells(aquifer) # anem::get_gridded_hydrodynamics() loc <- expand.grid(x=seq(0,10000,length.out=201),y=seq(-5000,5000,length.out=201)) thm <- theme(panel.background = element_rect(fill=NA,color=NA), panel.grid.major=element_blank(), panel.grid.minor=element_blank(), legend.position="none",axis.text=element_blank(), axis.ticks=element_blank(),axis.title=element_blank(), legend.title=element_text(size=9), legend.text=element_text(size=8), legend.spacing = unit(1,"mm"), legend.margin = margin(0.25,0,0,0.25,unit = "mm"), legend.box.margin = margin(0,0,0,0,unit="mm")) p_base <- ggplot() + geom_segment(data=aquifer$bounds %>% mutate(bt=factor(bound_type,levels=c("CH","NF","PB"), labels=c("River","No flow","Open"))) %>% filter(bt !="Open"), aes(x1,y1,xend=x2,yend=y2,color=bt),size=1) + scale_shape_manual("Well type",values = c(16,1,4,3)) + scale_color_manual("Boundary",values = c("blue","#555555","#CCCCCC"))
well_color <- "#AA3333" constant_flow <- tibble(y=seq(500,5000-500,length.out = 6),x=0) f1_a <- p_base + geom_segment(data=constant_flow,aes(x,y,xend=x+500,yend=y),arrow = arrow(length=unit(1,"mm"),type = "closed")) + geom_point(data=wells %>% filter(well_image=="Actual"), aes(x,y,shape=well_type),size=3,stroke=1) + # geom_point(data=wells %>% filter(well_image=="Actual",orig_wID==4), # aes(x,y,shape=well_type),size=3,stroke=1,color=well_color,show.legend = FALSE) + coord_equal() + labs(subtitle = "(a)") + thm %+replace% theme(legend.position = "left")
wells_roi <- wells %>% filter(orig_wID==4) %>% anem::gen_circles() f1_b <- p_base + geom_polygon(data=wells_roi,aes(x,y,group=id),fill=well_color,alpha=0.2,color=NA,size=0.25) + geom_point(data=wells %>% filter(well_image=="Actual",orig_wID!=4), aes(x,y,shape=well_type),size=3,stroke=1,alpha=0.5) + # geom_point(data=wells %>% filter(well_image=="Actual",orig_wID==4), # aes(x,y,shape=well_type),size=3,stroke=1,color=well_color) + geom_point(data=wells %>% filter(well_image=="Actual",orig_wID==4),aes(x,y), size=3,stroke=1,color=well_color,shape=16) + geom_point(data=wells %>% filter(grepl("\\+Q",well_image),orig_wID==4),aes(x,y), size=2.5,stroke=1,color=well_color,shape=15) + geom_point(data=wells %>% filter(grepl("\\-Q",well_image),orig_wID==4),aes(x,y), size=2,stroke=1,color=well_color,shape=0) + annotate("rect",xmin=3500,xmax=6500,ymin=-1500,ymax=1500,color="#992222",linetype="dashed",fill=NA,size=1) + annotate("point",500,-4500,shape=0,size=2,color=well_color,stroke=1) + annotate("text",1000,-4500,shape=0,size=4.5,label="Image wells", color="black",hjust=0,vjust=0.45) + coord_equal() + labs(subtitle = "(b)") + thm #%+replace% theme(legend.position="right")
streamfunction <- loc %>% dplyr::bind_cols(streamfunction=get_stream_function(loc,wells %>% filter(orig_wID==4),aquifer)) %>% dplyr::bind_cols(head=get_hydraulic_head(loc,wells,aquifer)) flowlines <- anem::get_contourlines(streamfunction %>% rename(z=streamfunction),nlevels=40) f1_c <- p_base + geom_path(data=flowlines,aes(x,y,group=line),color="#66CCEE",alpha=1) + geom_point(data=wells %>% filter(well_image=="Actual"),aes(x,y,shape=well_image), size=2,stroke=1,color="#AA3333") + geom_point(data=wells %>% filter(well_image!="Actual"),aes(x,y,shape=well_image), size=2,stroke=1,color="#AA3333") + annotate("rect",xmin=3500,xmax=6500,ymin=-1500,ymax=1500,color="#AA3333",linetype="dashed",fill=NA,size=1) + coord_equal(xlim=c(3500,6500),ylim=c(-1500,1500)) + scale_shape_manual(values=c(16,0,15))+ labs(subtitle = "(c)") + thm #%+replace% theme(legend.position = "right")
wells_actual <- wells %>% filter(well_image=="Actual") gridded <- get_gridded_hydrodynamics(wells,aquifer,head_dim=c(200,200),flow_dim=c(10,10)) missing_points <- do.call(rbind,lapply(split(wells_actual,seq(nrow(wells_actual))), function(loc,df) df[which.min(sqrt((df$x - loc$x)^2 + (df$y - loc$y)^2))[1],c("x","y")], df=gridded$flow)) x_multiplier <- 2e6 gridded$flow <- gridded$flow %>% mutate(mag=sqrt(dx^2 + dy^2), multiplier = (mag)^(1/4)*x_multiplier, x2 = x + dx * multiplier, y2 = y + dy * multiplier) %>% anti_join(missing_points,by=c("x","y")) # contours <- get_contourlines(gridded$head %>% rename(z=head_m),nlevels=20) with(gridded$flow, max(abs(x2-x))) f1_d <- ggplot() + geom_raster(data=gridded$head,aes(x,y,fill=head_m)) + # geom_path(data=contours,aes(x,y,group=line)) + geom_segment(data=aquifer$bounds %>% mutate(bt=factor(bound_type,levels=c("CH","NF","PB"), labels=c("River","No flow","Open (no boundary)"))) %>% filter(bound_type!="PB"), aes(x1,y1,xend=x2,yend=y2,color=bt),size=2,show.legend = FALSE) + geom_segment(data=gridded$flow,aes(x,y,xend=x2,yend=y2), arrow = arrow(ends="last",type="closed",length=unit(0.75,"mm")),color="black") + geom_point(data=wells %>% filter(well_image=="Actual"), aes(x,y,shape=well_type),size=3,stroke=1,show.legend = FALSE) + scale_shape_manual("Well type",values = c(16,1,4,3)) + scale_color_manual("Boundary",values = c("blue","#555555","#CCCCCC")) + scale_fill_viridis_c("Hydraulic\nhead, m",breaks=-1:2,labels=c("-1","0 (river)","1","2")) + labs(subtitle = "(d)") + coord_equal() + thm %+replace% theme(legend.position = "right",legend.key.height=unit(0.55,"cm"))
cap <- "(a) Aquifer scenario with four wells and arrows indicating constant background flow towards the river. (b) Well images and radii of influence for bottom-right well (red). (c) Reproduction of the no-flow and constant-head boundaries for single well using well images. (d) Hydraulic head and flow field for scenario. Note that the stream, flowing from top to bottom, goes from gaining to losing." # ggpubr::ggarrange(f1_a,f1_b,f1_c,f1_d,nrow=1,labels=c("(a)","(b)","(c)","(d)"),widths=c(0.32,0.2,0.2,0.32)) %>% # ggpubr::annotate_figure(bottom = "This is the caption.") # f1_a | f1_b | f1_c | f1_d gridExtra::grid.arrange(f1_a , f1_b , f1_c , f1_d,nrow=1,widths=c(0.33,0.2,0.2,0.35))
Understanding the method of images
The package relies on the principle of superposition and the method of images to generate groundwater hydrodynamics, illustrated by the two figures below (see ?get_segments_behavior and ?get_stream_function, respectively, which give examples for generating similar plots). The principle of superposition states that, because the groundwater equations are linear, the solution to the equations can be calculated as the sum of solutions. In the example below, the drawdown of the water table is calculated as the sum of drawdown due to three individual wells.
# define the aquifer aquifer_super <- aquifer_confined_example aquifer_super$Ksat <- 5e-4 # define wells and images wells_super_actual <- define_wells(data.frame(x=c(400,500,900),y=500,R=1500,Q=-0.1,diam=1)) wells_super <- generate_image_wells(wells_super_actual,aquifer_super) # set the segment seg <- c(0,500,1000,500) # get segment for each well and all wells seg_w1 <- get_segments_behavior(seg,wells_super[wells_super$orig_wID==1,],aquifer_super,1000) seg_w1 <- data.frame(seg_w1[,c("x","y","head")],Wells="W1") seg_w2 <- get_segments_behavior(seg,wells_super[wells_super$orig_wID==2,],aquifer_super,1000) seg_w2 <- data.frame(seg_w2[,c("x","y","head")],Wells="W2") seg_w3 <- get_segments_behavior(seg,wells_super[wells_super$orig_wID==3,],aquifer_super,1000) seg_w3 <- data.frame(seg_w3[,c("x","y","head")],Wells="W3") seg_all <- get_segments_behavior(seg,wells_super,aquifer_super,1000) seg_all <- data.frame(seg_all[,c("x","y","head")],Wells="All") # plot results seg_behavior <- rbind(seg_all,seg_w1,seg_w2,seg_w3) ggplot() + geom_line(data=seg_behavior,aes(x,head,linetype=Wells)) + geom_vline(data=data.frame(x=c(0,1000)),aes(xintercept=x),color="gray") + annotate(geom="text",x=40,y=86,label="Constant head boundary",angle=-90) + annotate(geom="text",x=1040,y=86,label="No flow boundary",angle=-90)
In the method of images, wells are mirrored across aquifer boundaries to reproduce characterics of the boundary, such as no flow or constant head. The following examples demonstrate (a) how a no-flow boundary is created by mirroring a well across the boundary, and (b) how a constant head boundary is recreated by mirroring the well and inverting the pumping.
# Create a grid of locations loc <- expand.grid(x=seq(-200,200,length.out=201),y=seq(-200,200,length.out=201)) # Constant head boundary wells_constant_head <- define_wells(x=c(-100,100),y=c(-0,0),Q=c(1e-2,-1e-2),diam=c(0.1,0.1),R=c(500,500)) aquifer <- define_aquifer("confined",1e-4,z0=20,h0=0,n=0.35) constant_head_boundary <- loc %>% bind_cols(streamfunction=get_stream_function(loc,wells_constant_head,aquifer)) %>% bind_cols(head=get_hydraulic_head(loc,wells_constant_head,aquifer)) # No flow boundary wells_no_flow <- define_wells(x=c(-100,100),y=c(-0,0),Q=c(-1e-2,-1e-2),diam=c(0.1,0.1),R=c(500,500)) no_flow_boundary <- loc %>% bind_cols(streamfunction=get_stream_function(loc,wells_no_flow,aquifer)) %>% bind_cols(head=get_hydraulic_head(loc,wells_no_flow,aquifer)) thm <- theme(axis.text=element_blank(),axis.title=element_blank(),axis.ticks=element_blank(), panel.background = element_rect(fill=NA,color="black"), panel.border = element_rect(fill=NA,color="black"), plot.title = element_text(size=10)) p_CH <- ggplot() + geom_contour(data=constant_head_boundary,aes(x,y,z=head),bins=50,linetype="dashed",color="black") + geom_contour(data=constant_head_boundary,aes(x,y,z=streamfunction),bins=25,color="darkblue") + geom_point(data=wells_constant_head,aes(x,y,fill=well_type),size=3,shape=21) + theme(legend.position=c(0.8,0.15),legend.title=element_blank()) + coord_equal() + thm + ggtitle("Constant head boundary at x = 0") p_NF <- ggplot() + geom_contour(data=no_flow_boundary,aes(x,y,z=head),bins=25,linetype="dashed",color="black") + geom_contour(data=no_flow_boundary,aes(x,y,z=streamfunction),bins=50,color="darkblue") + geom_point(data=wells_no_flow,aes(x,y,fill=well_type),size=3,shape=21) + theme(legend.position=c(0.8,0.1),legend.title=element_blank()) + coord_equal() + thm + ggtitle("No flow boundary at x = 0") gridExtra::grid.arrange(p_NF,p_CH,nrow=1)
A simple example of drawdown relationships
Load packages for the vignette
library(ggplot2) library(dplyr) library(anem)
Aquifers in this package are characterized by:
- Aquifer type, confined or unconfined
- Saturated hydraulic conductivity, Ksat
- Resting hydraulic head, h0
- Thickness of the aquifer, z0 (confined aquifers only)
- Porosity of the aquifer, n
- Aquifer boundaries, characterized by:
- Boundary ID, or bID
- Coordinates in a cartesian plane, given as slope, m, and intercept, b, and coordinates (x1, y1, x2, y2)
- Boundary type, NF, CH, or PB
- Additional user-defined parameters
Use the define_aquifer() function to create a new aquifer, which has the class aquifer. This class is essentially a list with names parameters and a separate print method.
# define aquifer bounds_df <- data.frame(bound_type=c("CH","NF","NF","NF"),m=c(Inf,0,Inf,0),b=c(0,1000,1000,0)) aquifer_unconfined <- define_aquifer("unconfined",1e-3,bounds=bounds_df,h0=100,n=0.35) print(aquifer_unconfined)
Wells are characterized by:
- Well ID, or wID
- Location x, y
- Pumping rate (+ for injection, - for pumping), Q
- Radius of influence, R
- Diameter, diam
- Well images (which are noted in the
well_imagecolumn) - Additional user-defined columns, including groups or weights.
The code below Defines 8 pumping wells using random locations and arbitrarily divides wells into countries "A" and "B" using a threshold at y = 500 m. The define_wells() function ensures that the wells have all appropriate columns, and the class of the returned object is a tibble (which functions like a data.frame, but has a couple bells and whistles). Also note that the radius of influence (R) is defined arbitrarily, but there is also a function to calculate this manually -- see ?get_ROI for more details.
The generate_image_wells() function generates well images to recreate the bounderies defined by the aquifer, aquifer_unconfined.
# define wells and well images set.seed(15) wells_actual <- define_wells(x = runif(8,0,1000), y = c(runif(4,0,500),runif(4,500,1000)), Q = -1/4, diam = 1, R = 1000, weights = 1, country = rep(c("A","B"),each=4)) wells <- wells_actual %>% generate_image_wells(aquifer_unconfined) print(wells)
Plot the aquifer and wells.
ggplot() + geom_segment(data=aquifer_unconfined$bounds,aes(x1,y1,xend=x2,yend=y2,linetype=bound_type)) + geom_abline(slope=0,intercept=500,linetype="dashed") + geom_point(data=wells_actual,aes(x,y,fill=country),shape=21) + coord_equal()
Get drawdown relationships using get_drawdown_relationships(). This function calculates the average drawdown at wells in each group defined by the column group_column. The average is taken as the weighted mean, determined by the weights_column. The weights were previously set equal for all wells so that the result here is a simple mean. The results show PHIii and PHIij.
drawdown_relationships <- get_drawdown_relationships(wells, aquifer_unconfined, group_column=country, weights_column=weights) drawdown_relationships
drawdown_relationships %>% knitr::kable("html") #%>% kableExtra::kable_styling()
Plot the head and flow
The hydrodynamics of the aquifer can be mapped by obtaining gridded heand flow using the get_gridded_hydrodynamics() function. The function takes as input the wells, aquifer, and grid dimensions for head and flow. It returns a list object with data.frames for head and flow, which can then be plotted.
gridded_1 <- get_gridded_hydrodynamics(wells,aquifer_unconfined,c(20,20),c(8,8)) ggplot() + geom_raster(data=gridded_1$head,aes(x,y,fill=head_m)) + geom_contour(data=gridded_1$head,aes(x,y,z=head_m),color="white",alpha=0.3) + geom_segment(data=gridded_1$flow,aes(x,y,xend=x2,yend=y2), arrow = arrow(ends="last",type="closed",length=unit(1,"mm")),color="black") + geom_segment(data=aquifer_unconfined$bounds,aes(x1,y1,xend=x2,yend=y2,linetype=bound_type)) + geom_point(data=wells %>% filter(wID==orig_wID),aes(x,y),shape=21) + coord_equal()
Check the head and flow along the boundaries
The function get_bounds_behavior() is a helper function to generate hydraulic properties at the boundaries. It obtains hydraulic head and the flow normal to the boundaries by setting boundaries to segments and using get_segments_behavior(). Normal flow is defined such that positive flow has some component in the x-direction (the y-direction depends on the normal vector for the boundary). The function plot_bounds_behavior() is a wrapper around get_bounds_behavior(), and it can be used to quickly compare the hydraulic head along the boundaries and flow across the boundaries under two scenarios:
- no images (ie, homogeneous, uniform aquifer -- all well images removed)
- images (ie, aquifer with boundaries -- all well images intact)
bounds_behavior <- plot_bounds_behavior(wells,aquifer_unconfined) gridExtra::grid.arrange(bounds_behavior$p_h,bounds_behavior$p_f,nrow=1)
Finally, the function plot_bounds_behavior also includes the raw data of head and flow used to create the above plots (bounds_behavior$bounds_behavior), as well as a summary of these values which includes mean head on each of the boundaries and mean flow as the mean of the absolute value of flow normal to each boundary (bounds_behavior$table). To numerically check that the boundaries are working as expected, we can print bounds_behavior$table:
bounds_behavior$table %>% knitr::kable("html") #%>% kableExtra::kable_styling()
Defining recharge / background flow
Background flow is parameterized by defining "recharge" in the aquifer. Recharge can be defined as a constant flow (F) or a recharge divide (D). In both cases, flow is defined as m^3^ / m-s, in the direction of the recharge vector. Recharge is therefore defined by specifying:
- Type of recharge: constant flow (F) or recharge divide (D)
- Recharge vector, pointing in the direction of flow
- Flow magnitude, in units of m^3^ / m-s
- Coordinates where the undisturbed water table (no wells) has a hydraulic head equal to aquifer$h0.
See ?define_recharge for more details.
# Define recharge for type "D" recharge_params <- list(recharge_type="D",recharge_vector=c(500,500,501,501), flow_main=1e-3,flow_opp=2e-3,x0=0,y0=0) # Define bounds as "PB" -- pervious bounds, which have no effect (other than to define output grid) bounds_recharge <- define_bounds(data.frame(bound_type=rep("PB",4),m=c(Inf,0,Inf,0),b=c(0,0,1000,1000))) aquifer_recharge <- define_aquifer("confined",1,h0=0,z0=1,n=0.35,recharge=recharge_params,bounds=bounds_recharge) gridded_2 <- get_gridded_hydrodynamics(wells_actual,aquifer_recharge,c(20,20),c(8,8))
The results can then be plotted, just as above:
ggplot() + geom_raster(data=gridded_2$head,aes(x,y,fill=head_m)) + geom_contour(data=gridded_2$head,aes(x,y,z=head_m),color="white",alpha=0.3) + geom_segment(data=gridded_2$flow,aes(x,y,xend=x2,yend=y2), arrow = arrow(ends="last",type="closed",length=unit(1,"mm")),color="black") + geom_segment(data=aquifer_recharge$bounds,aes(x1,y1,xend=x2,yend=y2,linetype=bound_type)) + geom_point(data=wells_actual,aes(x,y),shape=21) + coord_equal()
Particle tracking
The package implements particle tracking in two ways: (1) tracking individual particles using track_particles and (2) estimating capture zones of wells using get_capture_zone, which initializes particles in the vicinity of wells and tracks them in reverse. Here we use get_capture_zone on the previous aquifer to generate particle tracking away from wells.
# wells <- wells # wells[4,"Q"] <- 0.25 well_particles <- get_capture_zone(wells,aquifer_unconfined,t_max=365,n_particles=4)
gridded_3 <- get_gridded_hydrodynamics(wells,aquifer_unconfined) ggplot() + geom_raster(data=gridded_1$head,aes(x,y,fill=head_m)) + geom_contour(data=gridded_1$head,aes(x,y,z=head_m),color="white",alpha=0.3) + geom_segment(data=aquifer_unconfined$bounds,aes(x1,y1,xend=x2,yend=y2,linetype=bound_type)) + geom_path(data=well_particles,aes(x,y,group=i,color=as.factor(wID)),show.legend = FALSE) + geom_point(data=wells_actual,aes(x,y,color=as.factor(wID)),shape=21,show.legend = FALSE) + coord_equal()
Integrating with anem-app
The basic functionality of this package is implemented in a Shiny application, which can be run by calling anem_app() or accessed online at shinyapps.io. The Shiny application provides an GUI for setting aquifer properties, drawing boundaries, defining recharge, adding wells, and tracking particles. It also displays results directly in the app. You can also download scenarios created on anem-app and import them into R. To do so, check out ?import_app_rds from this package to learn more.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.