In gopalpenny/anem: Analytical element models for simple aquifers
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
# devtools::install_github("https://github.com/gopalpenny/anem")
library(tidyverse)
# library(anem)
devtools::load_all("~/Projects/R_packages/anem")
We need to validate the anem package for some simple use cases. In this validation we will consider:
- hydraulic head in a confined aquifer
- discharge potential in an unconfined aquifer
- flow direction
- behavior along boundaries
Hydraulic head in a confined aquifer
We assume that all flow within the aquifer follows the Depuit-Forschheimer assumption such that all flow is 2-Dimensional and parallel to the aquifer bottom, which is horizontal. Flow is governed by Darcy's law, which can be integrated to obtain the Thiem solution describing steady state hydraulic head and flow in the vicinity ($r<R$) of a pumping well:
\begin{align}
h - h_0 = -\frac{Q}{2 \pi T} \log \left( \frac{r}{R} \right)
\label{eq:thiemunconfined}
\end{align}
where $h$ is hydraulic head at a distance $r$ from a fully penetrating well that pumps at rate $Q$ (positive for injection, negative for abstraction), $h_0$ is the undisturbed hydraulic head when $Q=0$, $T$ is the transmissivity of the confined aquifer defined as saturated hydraulic conductivity multiplied by the thickness of the aquifer ($T= K_{sat} z_0$), and $R$ is the radius of influence of the pumping well, outside of which the hydraulic head is undisturbed by pumping ($h=h_0$). This solution is valid for steady state flow in confined aquifers.
To test the behavior, we define a confined aquifer with the following parameters and rectangular no-flow boundaries:
aquifer <- define_aquifer(
aquifer_type = "confined",
Ksat=1e-4,
z0=10,
bounds=data.frame(bound_type=rep("NF",4),m=c(Inf,0,Inf,0),b=c(0,0,1000,1000))
)
print(aquifer)
We define a single pumping well in the center of the aquifer, with a pumping radius just less than the width of the aquifer.
well <- define_wells(
Q=-1,
R=1200,
diam=1,
x=500,
y=500
) %>%
generate_image_wells(aquifer)
print(well)
And finally, we select three locations for observation, including one at the well:
loc <- tribble(~pID,~x,~y,
1, 750, 500,
2, 625, 750,
3, 500, 500)
The well and its images and radii of influence, along with observation locations, look like this:
p_full <- ggplot() +
geom_segment(data=aquifer$bounds,aes(x1,y1,xend=x2,yend=y2)) +
geom_polygon(data=gen_circles(well),aes(x,y,group=id),alpha=0.25) +
geom_point(data=well,aes(x,y,color=well_image)) +
geom_point(data=loc,aes(x,y,shape="obs location")) +
scale_shape_manual(values=4) +
coord_equal()
p_aquifer <- p_full + coord_equal(xlim=c(0,1000),ylim=c(0,1000))
gridExtra::grid.arrange(p_full,p_aquifer,nrow=1)
The image wells are labeled "Image (+Q)", indicating they have the same sign as the "Actual" well. If the "Actual" well is negative, the images will also have negative pumping rates. The first point (pID=1) is within the radius of influence of four wells. We can manually calculate $h-h_0$ using the Thiem solution for the influence of each well
loc_radii <- loc %>% crossing(well %>% rename(xw=x,yw=y)) %>%
mutate(r=sqrt((xw-x)^2+(yw-y)^2),r=if_else(r<diam/2,diam/2,r)) %>%
filter(r<R) %>%
mutate(`h-h0`=-Q/(2*pi*aquifer$Ksat*aquifer$z0)*log(r/R))
loc_radii %>% filter(pID==1,r<R) %>% select(pID,wID,Q,well_image,r,`h-h0`) %>%
knitr::kable("html",1)
loc_radii_net <- loc_radii %>% group_by(pID) %>% summarize(net=sum(`h-h0`))
The sum of the $h-h_0$ column is r round(loc_radii_net$net[loc_radii_net$pID==1],4). We can compare these results with those from anem and get_potential_differential():
get_potential_differential(loc[loc$pID==1,],well,aquifer)
We can compare the same results for the second observation point:
loc_radii %>% filter(pID==2,r<R) %>% select(pID,wID,Q,well_image,r,`h-h0`) %>%
knitr::kable("html",1)
The sum of the $h-h_0$ column is r round(loc_radii_net$net[loc_radii_net$pID==2],4).
get_potential_differential(loc[loc$pID==2,],well,aquifer)
And the results for the third observation point.
loc_radii %>%
filter(pID==3,r<R) %>% mutate(r=format(r, scientific=F)) %>%
select(pID,wID,Q,well_image,r,`h-h0`) %>%
knitr::kable("html",1,align=c('r','r','r','l','r','r'))
The sum of the $h-h_0$ column is r round(loc_radii_net$net[loc_radii_net$pID==3],3).
get_potential_differential(loc[loc$pID==3,],well,aquifer)
In this last case, the observation location is within the radius of the well casing, which means that the effect of the co-located well is measured a distance of $r=diam/2=0.5$. In this case, the results are also valid.
Hydraulic head in an unconfined aquifer
The solution can also be written for unconfined aquifers in terms of discharge potential, $\phi$:
$$\phi - \phi_0 = -\frac{Q}{\pi K_{sat}} \log \left( \frac{r}{R} \right),$$
where discharge potential is defined as $\phi=h^2$, and $\phi_0$ is the undisturbed discharge potential when $Q=0$.
Following the same analysis as above but redefining the aquifer as "unconfined", we further validate the package.
aquifer_uconfined <- define_aquifer(
aquifer_type="unconfined",
Ksat=1e-4,
bounds=aquifer$bounds
)
print(aquifer_uconfined)
The first point (pID=1) is within the radius of influence of four wells. We can manually calculate $h-h_0$ using the Thiem solution for the influence of each well
loc_radii <- loc %>% crossing(well %>% rename(xw=x,yw=y)) %>%
mutate(r=sqrt((xw-x)^2+(yw-y)^2),r=if_else(r<diam/2,diam/2,r)) %>%
filter(r<R) %>%
mutate(`phi-phi0`=-Q/(pi*aquifer_uconfined$Ksat)*log(r/R))
loc_radii %>% filter(pID==1,r<R) %>% select(pID,wID,Q,well_image,r,`phi-phi0`) %>%
knitr::kable("html",1)
loc_radii_net <- loc_radii %>% group_by(pID) %>% summarize(net=sum(`phi-phi0`))
The sum of the $\phi-\phi_0$ column is r round(loc_radii_net$net[loc_radii_net$pID==1],3). We can compare these results with those from anem and get_potential_differential():
get_potential_differential(loc[loc$pID==1,],well,aquifer_uconfined)
We can compare the same results for the second observation point:
loc_radii %>% filter(pID==2,r<R) %>% select(pID,wID,Q,well_image,r,`phi-phi0`) %>%
knitr::kable("html",1)
The sum of the $\phi-\phi_0$ column is r round(loc_radii_net$net[loc_radii_net$pID==2],3).
get_potential_differential(loc[loc$pID==2,],well,aquifer_uconfined)
And the results for the third observation point.
loc_radii %>%
filter(pID==3,r<R) %>% mutate(r=format(r, scientific=F)) %>%
select(pID,wID,Q,well_image,r,`phi-phi0`) %>%
knitr::kable("html",1,align=c('r','r','r','l','r','r'))
The sum of the $\phi-\phi_0$ column is r format(loc_radii_net$net[loc_radii_net$pID==3],scientific=F,nsmall=3).
get_potential_differential(loc[loc$pID==3,],well,aquifer_uconfined)
In this last case, the observation location is within the radius of the well casing, which means that the effect of the co-located well is measured a distance of $r=diam/2=0.5$. In this case, the results are also valid.
Particle tracking in a confined aquifer
Particle tracking performs better when particles are far away from wells. Close to wells the velocity changes quickly such that the numerical tracking routine does relatively poorly in this region.
Track a particle starting at $(x,y)=(10,0)$ moving away from an injection well at $(x,y)=(0,0)$. The grid cell size (x-direction) is 2.2 m with a step distance of 2 m.
wells <- data.frame(x=0,y=0,Q=1,diam=1,R=1000) %>% define_wells()
bounds_df <- data.frame(bound_type=rep("PB",4),m=c(Inf,0,Inf,0),b=c(-100,500,1000,-500))
aquifer <- define_aquifer(Ksat=1,n=0.3,z0=10,h0=0,aquifer_type="confined",bounds=bounds_df)
x0 <- 10
particle <- track_particles(c(x0,0),wells,aquifer,step_dist=2,grid_length = 500) %>% rename(x_num=x)
particle <- particle %>%
mutate(x_eqn=sqrt(wells$Q*time_days*3600*24/(pi*aquifer$n*aquifer$z0)+x0^2),
pct_diff=(x_num-x_eqn)/x_eqn * 100)
particle %>% filter(x_num<1000) %>%
gather(var,x,x_num,x_eqn,pct_diff) %>%
mutate(type=ifelse(grepl("x",var),"x","pct_diff")) %>%
ggplot() +
geom_line(aes(time_days,x,color=var,linetype=var),size=1.5) +
facet_wrap(~type,ncol=1,scales="free_y",strip.position = "left") +
theme(strip.placement = "outside",strip.background = element_rect(fill=NA),
axis.title.y=element_blank())
Track a particle starting at $(x,y)=(10,0)$ moving away from an injection well at $(x,y)=(0,0)$. The grid size is 11 m with a step distance of 20 m.
x0 <- 10
particle2 <- track_particles(c(x0,0),wells,aquifer,step_dist=20,grid_length = 100) %>% rename(x_num=x)
particle2 <- particle2%>%
mutate(x_eqn=sqrt(wells$Q*time_days*3600*24/(pi*aquifer$n*aquifer$z0)+x0^2),
pct_diff=(x_num-x_eqn)/x_eqn * 100)
particle2 %>%
filter(x_num<1000) %>%
gather(var,x,x_num,x_eqn,pct_diff) %>%
mutate(type=ifelse(grepl("x",var),"x","pct_diff")) %>%
ggplot() +
geom_line(aes(time_days,x,color=var,linetype=var),size=1) +
facet_wrap(~type,ncol=1,scales="free_y",strip.position = "left") +
theme(strip.placement = "outside",strip.background = element_rect(fill=NA),
axis.title.y=element_blank())
Particle tracking in an unconfined aquifer
I was unable to derive an equation for particle location in an unconfined aquifer. The integral of $dx/dt$ was quite gnarly. Instead, I compare the results of an unconfined aquifer with the theoretical results from a confined aquifer.
In this case with an initial undisturbed aquifer depth of 100 m, large $K_{sat}=1$ m s^-1^, and small pumping ($Q=0.001$ cumec), the comparison looks quite good when simulated over 5 years.
wells2 <- data.frame(x=0,y=0,Q=0.001,diam=1,R=1000) %>% define_wells()
aquifer_un <- define_aquifer(Ksat=1,n=0.3,h0=100,aquifer_type="unconfined",bounds=bounds_df)
x0 <- 10
particle_un <- track_particles(c(x0,0),wells2,aquifer_un,step_dist=2,grid_length = 500,t_max = 365*5) %>% rename(x_num=x)
particle_un <- particle_un %>%
mutate(x_eqn=sqrt(wells2$Q*time_days*3600*24/(pi*aquifer_un$n*aquifer_un$h0)+x0^2),
pct_diff=(x_num-x_eqn)/x_eqn * 100)
particle_un %>% filter(x_num<1000) %>%
gather(var,x,x_num,x_eqn,pct_diff) %>%
mutate(type=ifelse(grepl("x",var),"x","pct_diff")) %>%
ggplot() +
geom_line(aes(time_days,x,color=var,linetype=var),size=1.5) +
facet_wrap(~type,ncol=1,scales="free_y",strip.position = "left") +
theme(strip.placement = "outside",strip.background = element_rect(fill=NA),
axis.title.y=element_blank())
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 = "#>" )
# devtools::install_github("https://github.com/gopalpenny/anem") library(tidyverse) # library(anem) devtools::load_all("~/Projects/R_packages/anem")
We need to validate the anem package for some simple use cases. In this validation we will consider:
- hydraulic head in a confined aquifer
- discharge potential in an unconfined aquifer
- flow direction
- behavior along boundaries
Hydraulic head in a confined aquifer
We assume that all flow within the aquifer follows the Depuit-Forschheimer assumption such that all flow is 2-Dimensional and parallel to the aquifer bottom, which is horizontal. Flow is governed by Darcy's law, which can be integrated to obtain the Thiem solution describing steady state hydraulic head and flow in the vicinity ($r<R$) of a pumping well:
\begin{align} h - h_0 = -\frac{Q}{2 \pi T} \log \left( \frac{r}{R} \right) \label{eq:thiemunconfined} \end{align}
where $h$ is hydraulic head at a distance $r$ from a fully penetrating well that pumps at rate $Q$ (positive for injection, negative for abstraction), $h_0$ is the undisturbed hydraulic head when $Q=0$, $T$ is the transmissivity of the confined aquifer defined as saturated hydraulic conductivity multiplied by the thickness of the aquifer ($T= K_{sat} z_0$), and $R$ is the radius of influence of the pumping well, outside of which the hydraulic head is undisturbed by pumping ($h=h_0$). This solution is valid for steady state flow in confined aquifers.
To test the behavior, we define a confined aquifer with the following parameters and rectangular no-flow boundaries:
aquifer <- define_aquifer( aquifer_type = "confined", Ksat=1e-4, z0=10, bounds=data.frame(bound_type=rep("NF",4),m=c(Inf,0,Inf,0),b=c(0,0,1000,1000)) ) print(aquifer)
We define a single pumping well in the center of the aquifer, with a pumping radius just less than the width of the aquifer.
well <- define_wells( Q=-1, R=1200, diam=1, x=500, y=500 ) %>% generate_image_wells(aquifer) print(well)
And finally, we select three locations for observation, including one at the well:
loc <- tribble(~pID,~x,~y, 1, 750, 500, 2, 625, 750, 3, 500, 500)
The well and its images and radii of influence, along with observation locations, look like this:
p_full <- ggplot() + geom_segment(data=aquifer$bounds,aes(x1,y1,xend=x2,yend=y2)) + geom_polygon(data=gen_circles(well),aes(x,y,group=id),alpha=0.25) + geom_point(data=well,aes(x,y,color=well_image)) + geom_point(data=loc,aes(x,y,shape="obs location")) + scale_shape_manual(values=4) + coord_equal() p_aquifer <- p_full + coord_equal(xlim=c(0,1000),ylim=c(0,1000)) gridExtra::grid.arrange(p_full,p_aquifer,nrow=1)
The image wells are labeled "Image (+Q)", indicating they have the same sign as the "Actual" well. If the "Actual" well is negative, the images will also have negative pumping rates. The first point (pID=1) is within the radius of influence of four wells. We can manually calculate $h-h_0$ using the Thiem solution for the influence of each well
loc_radii <- loc %>% crossing(well %>% rename(xw=x,yw=y)) %>% mutate(r=sqrt((xw-x)^2+(yw-y)^2),r=if_else(r<diam/2,diam/2,r)) %>% filter(r<R) %>% mutate(`h-h0`=-Q/(2*pi*aquifer$Ksat*aquifer$z0)*log(r/R)) loc_radii %>% filter(pID==1,r<R) %>% select(pID,wID,Q,well_image,r,`h-h0`) %>% knitr::kable("html",1) loc_radii_net <- loc_radii %>% group_by(pID) %>% summarize(net=sum(`h-h0`))
The sum of the $h-h_0$ column is r round(loc_radii_net$net[loc_radii_net$pID==1],4). We can compare these results with those from anem and get_potential_differential():
get_potential_differential(loc[loc$pID==1,],well,aquifer)
We can compare the same results for the second observation point:
loc_radii %>% filter(pID==2,r<R) %>% select(pID,wID,Q,well_image,r,`h-h0`) %>% knitr::kable("html",1)
The sum of the $h-h_0$ column is r round(loc_radii_net$net[loc_radii_net$pID==2],4).
get_potential_differential(loc[loc$pID==2,],well,aquifer)
And the results for the third observation point.
loc_radii %>% filter(pID==3,r<R) %>% mutate(r=format(r, scientific=F)) %>% select(pID,wID,Q,well_image,r,`h-h0`) %>% knitr::kable("html",1,align=c('r','r','r','l','r','r'))
The sum of the $h-h_0$ column is r round(loc_radii_net$net[loc_radii_net$pID==3],3).
get_potential_differential(loc[loc$pID==3,],well,aquifer)
In this last case, the observation location is within the radius of the well casing, which means that the effect of the co-located well is measured a distance of $r=diam/2=0.5$. In this case, the results are also valid.
Hydraulic head in an unconfined aquifer
The solution can also be written for unconfined aquifers in terms of discharge potential, $\phi$: $$\phi - \phi_0 = -\frac{Q}{\pi K_{sat}} \log \left( \frac{r}{R} \right),$$ where discharge potential is defined as $\phi=h^2$, and $\phi_0$ is the undisturbed discharge potential when $Q=0$.
Following the same analysis as above but redefining the aquifer as "unconfined", we further validate the package.
aquifer_uconfined <- define_aquifer( aquifer_type="unconfined", Ksat=1e-4, bounds=aquifer$bounds ) print(aquifer_uconfined)
The first point (pID=1) is within the radius of influence of four wells. We can manually calculate $h-h_0$ using the Thiem solution for the influence of each well
loc_radii <- loc %>% crossing(well %>% rename(xw=x,yw=y)) %>% mutate(r=sqrt((xw-x)^2+(yw-y)^2),r=if_else(r<diam/2,diam/2,r)) %>% filter(r<R) %>% mutate(`phi-phi0`=-Q/(pi*aquifer_uconfined$Ksat)*log(r/R)) loc_radii %>% filter(pID==1,r<R) %>% select(pID,wID,Q,well_image,r,`phi-phi0`) %>% knitr::kable("html",1) loc_radii_net <- loc_radii %>% group_by(pID) %>% summarize(net=sum(`phi-phi0`))
The sum of the $\phi-\phi_0$ column is r round(loc_radii_net$net[loc_radii_net$pID==1],3). We can compare these results with those from anem and get_potential_differential():
get_potential_differential(loc[loc$pID==1,],well,aquifer_uconfined)
We can compare the same results for the second observation point:
loc_radii %>% filter(pID==2,r<R) %>% select(pID,wID,Q,well_image,r,`phi-phi0`) %>% knitr::kable("html",1)
The sum of the $\phi-\phi_0$ column is r round(loc_radii_net$net[loc_radii_net$pID==2],3).
get_potential_differential(loc[loc$pID==2,],well,aquifer_uconfined)
And the results for the third observation point.
loc_radii %>% filter(pID==3,r<R) %>% mutate(r=format(r, scientific=F)) %>% select(pID,wID,Q,well_image,r,`phi-phi0`) %>% knitr::kable("html",1,align=c('r','r','r','l','r','r'))
The sum of the $\phi-\phi_0$ column is r format(loc_radii_net$net[loc_radii_net$pID==3],scientific=F,nsmall=3).
get_potential_differential(loc[loc$pID==3,],well,aquifer_uconfined)
In this last case, the observation location is within the radius of the well casing, which means that the effect of the co-located well is measured a distance of $r=diam/2=0.5$. In this case, the results are also valid.
Particle tracking in a confined aquifer
Particle tracking performs better when particles are far away from wells. Close to wells the velocity changes quickly such that the numerical tracking routine does relatively poorly in this region.
Track a particle starting at $(x,y)=(10,0)$ moving away from an injection well at $(x,y)=(0,0)$. The grid cell size (x-direction) is 2.2 m with a step distance of 2 m.
wells <- data.frame(x=0,y=0,Q=1,diam=1,R=1000) %>% define_wells() bounds_df <- data.frame(bound_type=rep("PB",4),m=c(Inf,0,Inf,0),b=c(-100,500,1000,-500)) aquifer <- define_aquifer(Ksat=1,n=0.3,z0=10,h0=0,aquifer_type="confined",bounds=bounds_df) x0 <- 10 particle <- track_particles(c(x0,0),wells,aquifer,step_dist=2,grid_length = 500) %>% rename(x_num=x) particle <- particle %>% mutate(x_eqn=sqrt(wells$Q*time_days*3600*24/(pi*aquifer$n*aquifer$z0)+x0^2), pct_diff=(x_num-x_eqn)/x_eqn * 100)
particle %>% filter(x_num<1000) %>% gather(var,x,x_num,x_eqn,pct_diff) %>% mutate(type=ifelse(grepl("x",var),"x","pct_diff")) %>% ggplot() + geom_line(aes(time_days,x,color=var,linetype=var),size=1.5) + facet_wrap(~type,ncol=1,scales="free_y",strip.position = "left") + theme(strip.placement = "outside",strip.background = element_rect(fill=NA), axis.title.y=element_blank())
Track a particle starting at $(x,y)=(10,0)$ moving away from an injection well at $(x,y)=(0,0)$. The grid size is 11 m with a step distance of 20 m.
x0 <- 10 particle2 <- track_particles(c(x0,0),wells,aquifer,step_dist=20,grid_length = 100) %>% rename(x_num=x) particle2 <- particle2%>% mutate(x_eqn=sqrt(wells$Q*time_days*3600*24/(pi*aquifer$n*aquifer$z0)+x0^2), pct_diff=(x_num-x_eqn)/x_eqn * 100)
particle2 %>% filter(x_num<1000) %>% gather(var,x,x_num,x_eqn,pct_diff) %>% mutate(type=ifelse(grepl("x",var),"x","pct_diff")) %>% ggplot() + geom_line(aes(time_days,x,color=var,linetype=var),size=1) + facet_wrap(~type,ncol=1,scales="free_y",strip.position = "left") + theme(strip.placement = "outside",strip.background = element_rect(fill=NA), axis.title.y=element_blank())
Particle tracking in an unconfined aquifer
I was unable to derive an equation for particle location in an unconfined aquifer. The integral of $dx/dt$ was quite gnarly. Instead, I compare the results of an unconfined aquifer with the theoretical results from a confined aquifer.
In this case with an initial undisturbed aquifer depth of 100 m, large $K_{sat}=1$ m s^-1^, and small pumping ($Q=0.001$ cumec), the comparison looks quite good when simulated over 5 years.
wells2 <- data.frame(x=0,y=0,Q=0.001,diam=1,R=1000) %>% define_wells() aquifer_un <- define_aquifer(Ksat=1,n=0.3,h0=100,aquifer_type="unconfined",bounds=bounds_df) x0 <- 10 particle_un <- track_particles(c(x0,0),wells2,aquifer_un,step_dist=2,grid_length = 500,t_max = 365*5) %>% rename(x_num=x) particle_un <- particle_un %>% mutate(x_eqn=sqrt(wells2$Q*time_days*3600*24/(pi*aquifer_un$n*aquifer_un$h0)+x0^2), pct_diff=(x_num-x_eqn)/x_eqn * 100)
particle_un %>% filter(x_num<1000) %>% gather(var,x,x_num,x_eqn,pct_diff) %>% mutate(type=ifelse(grepl("x",var),"x","pct_diff")) %>% ggplot() + geom_line(aes(time_days,x,color=var,linetype=var),size=1.5) + facet_wrap(~type,ncol=1,scales="free_y",strip.position = "left") + theme(strip.placement = "outside",strip.background = element_rect(fill=NA), axis.title.y=element_blank())
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.