In gopalpenny/anem: Analytical element models for simple aquifers

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width=4, 
  fig.height=2.5
)
knitr::opts_knit$set(
  eval.after = 'fig.cap'
)

library(anem)
library(tidyverse)
library(patchwork)

Abstract

Numerous tools exist for simulating groundwater flow via finite element models (e.g., MODFLOW) or analytic element models (e.g., gflow). These tools allow for analysis of complex groundwater behavior, but such functionality that may not be necessary in projects that demand simple groundwater models combined with (a) ease of use or (b) advanced statistical or theoretical analysis. We address each of these needs through the development of a web application and R package that simulate 2-dimensional, steady-state groundwater flow in confined and unconfined aquifers with rectangular boundaries using analytic elements and the method of well images. These platform-independent toosl offer functionality for defining aquifers and pumping wells, determining hydraulic head and lateral flow, calculating drawdown relationships among groups of wells, and evaluating particle transport. The ease of use of the web application should be especially valuabe in low-budget applications, such as small water management agencies or educational settings. Integration of groundwater models with R can provide numerous benefits for quickly simulating and analyzing groundwater hydrology problems. Simulation of idealized aquifer scenarios can serve as a quick-and-dirty first-cut analysis or be used to simulate synthetic scenarios as a way to explore variability and uncertainty. The functionality of the R package can be easily analyzed, mapped, and incorporated into more complex models in R. Indeed, the linear behavior of groundwater flow and principle of superposition means that physically plausible groundwater dynamics can be incorporated into more complex models of hydrology and coupled human-natural systems. We validate these tools and describe multiple scenarios in which the web application and R package might be used.

Introduction

Groundwater plays an essential role in environmental and water resources management. In various contexts it may serve as a resource, hazard, or transport medium. As such, groundwater modeling is ubiquitous for a variety of applications in industry, academia, and policy. The fundamental principles of groundwater dynamics are mathematically simple and commonly taught in civil and environmental engineering courses. However, there remains an uncessary learning curve to bridge the divide between the simplicity of groundwater equations and the available tools for modeling groundwater.

Groundwater models can be broadly categorized into numerical and analytical approaches. Numerical models simulate groundwater behavior using a finite-element representation of Richard's equation. Groundwater hydraulics are simulated on a discrete grid by stepping forward through time such that each time step solves the discrete version of Richard's equation. The most common numerical modeling platform, MODFLOW, is freely available and runf on both Windows and Unix operating systems.

Analytical groundwater models utilize the analytical element method, whereby two-dimensonal, steady-state versions of the groundwater equations are combined to reproduce groundwater behavior under simplifying assumptions. Groundwater dynamics depend on the sum of analytical solutions to groundwater scenarios containing individual elements.

The tools described here fill a missing gap in the available options for groundwater modeling. Existing modeling tools should suffice when the stakes are high or sufficient resources are available, but may not be appropriate for scenarios requiring inexpensive first-cut solutions. For instance, graphical setup in MODFLOW is currently only available on Windows systems. A recently deployed tool, MAGNET4Water offers online access and potential applications across platforms. Of the freely-available analytical models, gflow, XX, and ;lkj operate only on Windows operating systems. The web app provides an easy-to-use interface for quick-and-dirty modeling on any computing platform with web browsers and internet access.

Another requirement is the applications of complex analysis of aquifers. One exception is the python package by Bakker (timml). The R package we provide provides groundwater modeling in R coupled with the vast analytical tools within R.

We present two related platforms for groundwater modeling that offer additional flexibility and make groundwater modeling more readily accessible for simple scenarios. The first is an R package, anem, which models simple groundwater scenarios. 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. It provides speicific functionality 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.

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 package allows for constant background flow, any number of pumping or injection wells, as well as particle tracking. 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.

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 <- crossing(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=14),
        legend.text=element_text(size=12),
        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"))

R package

Analytic elements

Principle of superposition

The analytic element method represents the aquifer domain as the sum of analytical solutions of multiple elements added to the aquifer. We focus on 2-dimensional aquifers that reasonably adhere to the Depuit-Forscheimer assumption, which states that flow is parallel to the bottom of the aquifer. For a 2D confined aquifer with uniform thickness, the groundwater continuity equation can be written as

\begin{align} \nabla^2 h = 0, \end{align}

where $h$ is the hydraulic head. Because this equation is linear for $h$, it allows for the principle of superposition, where $f(a h_1 + b h_2) = a f(h_1) + b f(h_2)$. This means that the effect of pumping on hydraulic head from multiple wells can be calculated as the sum of the effects from individual wells, which can be estimated separately.

In 2D unconfined aquifers, the groundwater continuity equation can be written as

\begin{align} \nabla^2 \phi = 0, \end{align}

where $\phi = h^2$ is discharge potential and $h$ is the thickness of the water table. The principle of superposition once again applies because $f(a \phi_1 + b \phi_2) = a f(\phi_1) + b f(\phi_2)$. In this case, the the effect of each analytic element on discharge potential is calculated separately and the results are summed. In this way, the aquifer scenario is represented by the combined effect of recharge, well pumping, and aquifer boundaries.

Recharge

In steady steady state, recharge can be represented by constant background flow. In confined aquifers, this requires setting a constant flow defined using Darcy's law as $Q = -z K_{sat} \nabla h$, where $z$ is defined as the thickness of the aquifer $(z_0)$ for confined aquifers or the thickness of the water table $(h)$ for unconfined aquifers. Integrating this equation for confined aquifers in the x-direction yields

\begin{align} h - h_0 = \frac{Q_x}{K_{sat} z_0}(x_0 - x), \end{align}

where $h_0$ and $x_0$ refer to the head at some datum. The head gradient can be written as

\begin{align} \frac{\partial h}{\partial x} &= -\frac{Q_x}{ K_{sat} z_0} \end{align}

For unconfined aquifers, integration results in terms of discharge potential are

\begin{align} \phi - \phi_0 &= \frac{2 Q_x}{K_{sat}} (x_0 - x), \ \frac{\partial \phi}{\partial x} &= -\frac{2 Q_x}{K_{sat}}. \end{align}

Following the principle of superposition, the effect of recharge on hydraulic head or discharge potential can be calculated independently of other analytic elements.

Pumping from a single well

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. Generally speaking, this means that the thickness of the aquifer must be considerably less than the lateral extent of the aquifer, and wells must be screened from the bottom to the top of the wetted portion of the aquifer. Flow is governed by Darcy's law which, in the case of radial flow to a well, can be integrated to obtain the Thiem solution describing steady state hydraulic head in the vicinity ($r \le R$) of a pumping well:

\begin{align} h - h_0 &= -\frac{Q}{2 \pi T} \log \left( \frac{r}{R} \right) \ \frac{\partial h}{\partial r} &= -\frac{Q}{2 \pi T r} \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. The head gradient in this case, focusing on the x-direction, is:

\begin{align} \frac{\partial h}{\partial x} &= -\frac{Q}{2 \pi T} \frac{x-x_i}{r^2}, \end{align}

where $x$ is the observation location and $x_i$ is the location of the well. The solution can also be written for unconfined aquifers as:

\begin{align} \phi - \phi_0 &= -\frac{Q}{\pi K_{sat}} \log \left( \frac{r}{R} \right) \ \frac{\partial \phi}{\partial x} &= -\frac{Q}{\pi K_{sat}} \frac{x-x_i}{r^2} \end{align}

where $\phi_0$ is the undisturbed discharge potential when $Q=0$. Although the radius of influence does not affect the shape of the cone of depression, it determines how far it spreads from the well. Multiple approaches for calculating this radius have been proposed for confined and unconfined aquifers (CITES). The package implements the (CITE) approach for confined aquifers and (CITE) approach for unconfined aquifers.

Boundary conditions

In the analytic element method, the effect of boundary realized using strategic placement of wells that reproduce the conditions of the boundary. For instance, two wells with equivalent pumping mimic a no-flow boundary along the line of equal distance, because no flow occurs across this line (Figure 1c, horizontal boundary). Two wells pumping at an equivalent rate but opposite sign mimic a constant head boundary. In this case, the flow lines are perpendicular to the equal distance line meaning that no flow occurs along the boundary and it is of a constant head (Figure 1c, vertical boundary). This approach can be used to exactly reproduce no-flow and constant head boundaries by successively imaging wells across the boundaries.

cap <- "Figure 1: Basic functionality of anem R package. (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

Flow direction and velocity

According to Darcy's law, flow occurs in the direction opposite the gradient of the hydraulic head, which can be determined analytically using the analytical element method. Because both hydraulic head and discharge potential follow the principle of superposition, it follows that their derivatives also follows the principle of superposition. Furthermore, using the chain rule, $\partial h/\partial l = 1/(2h) \partial \phi/\partial l$, meaning that the hydraulic gradient can be calculated analytically in confined and unconfined aquifers.

Flow is then calculated from the hydraulic gratdient using Darcy's law, $q_l = - K_{sat} \partial \phi/\partial l$, and the flow velocity is equal to $v_l = q_l / n$, where $n$ is porosity of the aquifer.

Package description

Simulating groundwater hydraulics

The R package that we have developed applies the method of images in the aquifer scenarios described above. Aquifers are defined by their type (confined or unconfined), saturated hydraulic conductivity ($K_{sat}$), porosity ($n$), undisturbed head (i.e., the hydraulic head without pumping, $h_0$), and thickness ($z_0$, for confined aquifers only). Aquifers can also contain backround flow as a uniform flow field, or as two distince uniform fields on either side of a divide.

Wells are defined by a location, pumping rate, diameter, and radius of influence, and imaged across the aquifer borders contained within the radius of influence. The package allows users to calculate hydraulic head and flow direction at any number of locations with a single function call. The effect of each well on hydraulic head is calculated individually using the equations above, and the results are summed. Flow direction is calculated as the derivative of hydraulic head at any given location, proceeding analytically for both confined and unconfined aquifers.

Drawdown relationships

One of the most important advantages of this package is the ability to quickly calculate the relationship between pumping and drawdown relationships among groups of wells. Here we define the drawdown relationship as the average drawdown in a group of wells ($N \ge 1$) due to unit pumping from another group of wells. Results are given in terms of head (m) for confined aquifers or discharge potential (m^2^) for unconfined aquifer. The average drawdown and unit pumping can be weighted to favor certain wells within each group, if desired. This framing of drawdown into a single parameter between groups of wells allows for realistic framing of groundwater within much more complex models. For instance, such a framing was applied by Müller et al. (2018?) in using non-cooperative game theory to simulate cooperation in transboundary aquifers. Additional applications could include dynamical systems, where farmers change their abstraction ($Q_i$) and irrigated crop area ($I_i$) based on water levels within their wells ($h_i$). In this case, water level for user $i$ is the sum of drawdown relationships due to all the other users ($D_{ij}$) times their pumping ($Q_j$). If crops are similar, then we can set $Q_i = c I_i$ and write the following governing equation for the dynamical system:

\begin{align} \dot{I_i} = \sum a D_{ij} Q_j + b I_i \end{align}

In other words, this framework of drawdown relationships provides simple ways to incorporate groundater behavior into complex models in a realistic framework, and these relationships are easily estimated with the anem package.

Particle tracking

The package allows for particle tracking, either by initializing particles at a user-specified location or calculating capture zones around pumping wells. Particle tracking is simulated numerically by forward stepping in time using the Euler method. First the aquifer domain is discretized and fluid velocities are calculated on a grid. Instananeous particle velocity is calculated using bilinear interpolation of the grid at the particle location, unless the particle is inside a grid cell adjacent to a well. In this case, the fluid velocity is calculated at the particle location (i.e., without bilinear interpolation) to ensure precision in the vicinity of wells. Time stepping is variable, such that the particle travels a fixed distance each time step. Finally, well capture zones can be calculated by initializing a group of particles around each well and running particle tracking in reverse.

Web application

The web application incorporates the essential functionality of the R package, while providing an interactive environment to build and simulate a given scenario. There is a base map to assist with placing wells and aquifer components on the map. Particle tracking is also enabled through point tracking or well capture zones. The results can be interactively updated by changing well or aquifer characteristics. Finally, scenarios can downloaded from the application as R data files. These files can then be shared and uploaded to the application, allowing for reproducability and collaboration on the same scenarios.

Use cases

The following "use cases" demonstrate how the application might be used by planners, academics, and educators.

Case 1: Locating the source of urban contamination

Web application

A municipality receives complaints about the quality of its municipal water supply. Testing data suggests the possibility of contamination of the unconfined source aquifer below the city; the municipality wishes to establish, quickly and inexpensively, possible points of origin of the contamination. In this example, W1, W2, and W3 are municipal wells, each pumping 5 MGD (0.22 cumec); W1 is the contaminated well. Around these wells are three sources of discharge: D1, D2, and D3, which are injection wells which each operate at 0.45 MGD (0.02 cumec); and S1 and S2 represent points of known surface discharge. The area of interest is bounded by a river to the northeast, represented by a constant-head boundary. Notice that several of the discharge points are between the features of interest in the model (the municipal wells and the constant-head boundary) such that the municipality, at the outset, is plausibly uncertain as to the relative effects of pumping at the municipal abstraction wells and of the constant-head boundary, but has can make a reasonable estimate of the lateral flow toward the river.

We implemented this imaginary scenario in the web application as an aquifer bounded on one side with a constant head boundary, later recharge flowing towards the river, and both the pumping and injection wells (Figure XX). Well capture zones were simulated for a maximum of 10 years, after which all particle terminated at one of the boundaries or an injection well. The contaminated well draws water from all of the injection wells and one of the surface contamination sites, meaning there are four possible sources that could contaminate the well in question. However, the suspected sites can be reduced further given that the three injection wells are captured by multiple wells. No contamination has appeared at either of the other municipal pumping sites, which indicates that D1, D2, and D3 can be elimated from consideration. The municipality can begin by investigating the most likely source of contamination, S1.

Education opportunity

In addition to the practical application described above, this scenario could be used as an educational opportunity. For instance, consider the following game in which students must design a plan for capturing the contamination, using the minimum amount of pumping. This will involve adding pumping wells to the map and ensuring none of the lines touch the main city well. Each well costs \$1000, and the cost of pumping is \$3/cumec/well.

Case 2: Application for pumping well

A groundwater management agency has received an application to drill a new well at site W. The agency’s mandate prohibits new wells from withdrawing quantities of water that will “materially impair” the withdrawals of existing wells, which the agency interprets as a diminution of more than 5% from a known baseline. The agency wishes to determine roughly how much water can be withdrawn at W without impairing withdrawals at existing wells C1, C2, and C3. The area of interest is located near a river with the wells extracting a from a long but narrow confined aquifer connected to the stream.

Web application

The boundaries, aquifer, and wells are prepared in the web application (Figure Xa), with a constant-head boundary mimicking the river and no-flow parallel boundaries on either side. Hydraulic heads with and without the new well are calculated at each of the exisint wells, and the percent reduction is recorded. At the existing wells, the maximum reduction is 2.6%, and the well can be approved. If the agency has uncertainty about the parameters, they can adjust them manually in the app and check the results.

# MGD to cumec... L / s
q_mgd <- 1 # * 1M * m^3/gal / 24 hr/day / 3600 sec / hr
q_mgd * 5e5 * 0.00378541 / 24 / 3600

# L / s to cumec -- > /1000

case3 <- tibble::tibble(Well=c("E1","E2","E3","W"),
                        Baseline=c(48.451,48.27,47.931,49.066)) %>%
  dplyr::mutate(Minimum=Baseline* 0.95,
                # `Head (Q = -4)`=c(47.241,45.006,44.548,48.572),
                `Head (Q = -9)`=c(47.74,47.026,45.686,42.208),
                `Head (Q = -10)`=c(47.66,46.887,45.437,41.446)) %>%
  dplyr::mutate(#`% change (Q = -4)`=round((Baseline-`Head (Q = -4)`)/Baseline * 100,2),
                `% change (Q = -9)`=round((Baseline-`Head (Q = -9)`)/Baseline * 100,2),
                `% change (Q = -10)`=round((Baseline-`Head (Q = -10)`)/Baseline * 100,2)) %>%
  dplyr::arrange(Well) %>% dplyr::filter(!grepl("D",Well)) %>% dplyr::select(-dplyr::contains("Q = -4"),-dplyr::contains("Head"))
case3 %>% xtable::xtable(caption="Hydraulic head at existing pumping wells.",label="t:gw-agency") %>% print(include.rownames=FALSE,caption.placement="top")
case3 %>% knitr::kable() %>% kableExtra::kable_styling()

# urban contamination check
house_heads <- tibble::tibble(wID = c(14,15,16,17),
                          House = paste0("H",1:4),
                          `Head, with pumping`=c(28.018, 28.082, 27.542, 27.638),
                          `Head, no pumping`=c(31.734, 31.85, 31.522, 31.373)) %>% 
  dplyr::mutate(`Change in head, m` = `Head, no pumping`-`Head, with pumping`)
# house_heads %>% knitr::kable() %>% kableExtra::kable_styling()
house_heads %>% dplyr::select(-wID) %>% xtable::xtable(caption="Change in head at H1-H4 due to cessation of W2 pumping.",label="t:no-pumping") %>% print(include.rownames=FALSE,caption.placement="top")

Monte Carlo with R application

Some applications in industry or academia may demand a more complete analysis of uncertainty. In this case, let's assume there is uncertainty regarding saturated hyraulic conducitivity, $(K_{sat})$, aquifer thickness, $(z_0)$, and porosity $(n)$. In this synthetic scenario, we assume uniform distrubtions for each property, with $K_{sat}$ in the range $[1\times10^{-4},2\times10^{-3}]$ m/s, $z_0$ in $[10, 20]$ m, and $n$ in $[0.35,0.4]$.

The scenario prepared using the web application can be exported to Rdata file and imported into R for analysis. Monte carlo analysis was conducted with N = 1000 aquifers with randomly generated parameters. Uncertainty in the radius of influence of each well was accounted for by using the randomly generated parameters to re-calculate this radius using Cooper and Jacob () for each random realization. Pumping in existing wells was set to zero. For each random scenario, the maximum decline in head in existing wells was calculated and compared to the head without pumping from well D. The likelihood of of greater reduction in head increases with increasing pumping (Figure Xb), and the probability of compliance decreases (Figure Xc).

Case 3: Industrial plant going out of commission

Education opportunity

An industrial plant has been operating for over 50 years but now plans to close. The plant uses considerable amount of water. Because the plant is in a rural area, their water must be pumped directly from the ground. A local manager might reasonably wonder about the effect of closing the plant on the local water table. In some cases, closing of large industrial operations has caused flooding in nearby houses because the cessation of abstraction causes an increase in the local water table. To assess the situation, a planner wants to use anem-app to conduct an initial assessment of the sitation.

place an abstraction well to represent the plant and a number of observation wells (with zero pumping rate) in locations of the nearest houses.

Use scenario in "plant_closure.rds". What happens when the plant closes and how will this affect people in nearby houses? Is there anything you can do to keep this from happening, including the possibility of increasing pumping in the houses or drilling new wells?

gopalpenny/anem documentation built on Dec. 20, 2020, 5:27 a.m.