Solution method

For steady-state, single-layer Dupuit-Forchheimer flow in the horizontal plane, the governing equation is:

\[ T\frac{\partial^2{h}}{\partial x^2} + T\frac{\partial^2{h}}{\partial y^2} + N = 0 \]

with \(T\) the aquifer transmissivity \([L^2/T]\), equal to \(kH\) with \(k\) the horizontal hydraulic conductivity \([L/T]\) and \(H\) the saturated thickness \([L]\), \(h\) the hydraulic head \([L]\) and \(N\) a sink/source term \([L/T]\).

When the aquifer is phreatic, the saturated thickness \(H\) is variable and depends on \(h\), which renders the flow equation non-linear when solving for \(h\). Therefore, the discharge potential \(\Phi~[L^3/T]\) is introduced so that the governing flow equation becomes:

\[ \frac{\partial^2{\Phi}}{\partial x^2} + \frac{\partial^2{\Phi}}{\partial y^2} + N = 0 \]

which is a linear equation that can be solved for \(\Phi\). The resulting hydraulic heads can then be computed from:

\[ \begin{aligned} \Phi &= kHh - (\frac{1}{2}kH^2 + kHz_b) &\qquad h \ge z_t \\ \Phi &= \frac{1}{2}k(h - z_b)^2 &\qquad z_b \lt h \lt z_t \end{aligned} \]

with \(z_t\) and \(z_b\) the aquifer top and bottom elevation \([L]\), respectively.

In addition to the discharge potential \(\Phi\), there exists a stream function \(\Psi~[L^3/T]\) which is defined as a function that is constant along streamlines and whose difference between two points is equal to the amount of water flowing between those two points. By convention, the stream function increases to the right when looking in the direction of flow. It can be used to visualize the amount of flow in a groundwater system, e.g. when constructing flow nets. \(\Psi\) is only defined for the Laplace equation, i.e. when \(N\) is zero. In modeling terms, this means that it is only specified outside area-sinks (see paragraph Area-sink below).

The discharge potential \(\Phi\) and the stream function \(\Psi\) fulfill the Cauchy-Riemann conditions and can therefore be combined into a complex potential \(\Omega\) whose real and imaginary terms are \(\Phi\) and \(\Psi\), respectively. By solving the flow problem in terms of \(\Omega\), \(\Phi\) and \(\Psi\) are solved simultaneously in the complex plane, which is described by the complex coordinate \(\zeta=x+iy\). The solution in terms of \(\Omega\) yields \(\Phi\) and eventually \(h\).

The negative derivative of \(\Omega\) with respect to \(\zeta\) is the complex discharge \(W\), whose real component is \(Q_x\) and whose negative imaginary component is \(Q_y\), i.e. the \(x\) and \(y\) components of the discharge vector \(Q~[L^2/T]\). \(Q\) represents the amount of flow in the aquifer integrated along the vertical (or, alternatively, the first derivative of \(\Phi\) in the \(x\) and \(y\) directions). The vertical component \(Q_z\) is determined from mass balance considerations as discussed above. The Darcy flux \(q~[L/T]\) (also called specific discharge or Darcy velocity) is computed as \(Q/H\). The average linear flow velocity \(v~[L/T]\) is computed as \(q/n\) with \(n\) the aquifer porosity, which should represent the effective porosity.

Each element gives a solution for the governing flow equation in terms of \(\Omega\), as well as for the complex discharge \(W\). An element has one or more free parameters, e.g. the discharge for a well element. These parameters may be user-specified or are a priori unknown, in which case they need to be computed from other conditions, e.g. a specified head. The parameters are specified at so called collocation points or control points. Since the solution for a given unknown element parameter depends on the results of all other elements, the complete system of equations needs to be solved simultaneously for all unknowns. This is done by setting up the flow equations at the collocation points in a matrix form of \(Ax=b\), which is then solved for the unknown parameter vector \(x\). Once the element parameters are known, the contributions of each element to the complex potential \(\Omega\) and the complex discharge \(W\) at any point can be calculated and superimposed to form the resulting \(\Omega\) and \(W\) values, which then yield the variables of interest such as \(h\) and \(Q\).

For a more complete overview of analytic element modeling, the reader is referred to Haitjema (1995), Strack (1989), Strack (2003) and Hunt (2006).

Well

The simplest element is a well with a specified discharge \(Q\), which can be created with the well() function. It requires the \(x\) and \(y\) coordinate of the well and the pumping rate \(Q\) (positive for withdrawal). The well radius \(r_w\) can also be specified (default is \(0.3\) length units). Inside the well annulus as defined by \(r_w,\) the solution is undefined. When an output variable is requested inside this circle, its location is projected onto the nearest point on the outer annulus. It is assumed that the well is screened across the full (saturated) aquifer.

As an example, we specify two extraction wells. We then create the model with a call to aem(). This function needs the aquifer parameters and the elements. You can supply the elements individually or put them into a single list(). The aem() function automatically solves the system of equations by calling solve.aem(). By default, the model has a variable saturated thickness which depends on the calculated head. You can change this by setting type = 'confined' to force the model to use a constant saturated thickness, set equal to the aquifer thickness. This may be useful when you want to force the system of equations to remain linear when elements with a hydraulic resistance are supplied (see below).

Finally, we create the contouring grid and plot the head contours in blue and the streamlines in orange using contours(). Note the ugly branch cuts (in AEM jargon) of the streamlines extending from the well centers towards the negative x-axis, which is an artifact of contouring the stream function.

# create well elements
w1 = well(x = 200, y = 0, Q = 300)   
w2 = well(x = -200, y = 0, Q = 1000)

# create the model. This automatically solves the system of equations.
m = aem(k = k, top = top, base = base, n = n, w1, w2)

# set up the contouring grid 
xg = seq(-600, 600, length = 100)
yg = seq(-300, 300, length = 100)

# plot head contours and streamlines
contours(m, xg, yg, col = 'dodgerblue', nlevels = 20, drawlabels = FALSE)
contours(m, xg, yg, variable = 'streamfunction', col = 'orange',
         nlevels = 20, drawlabels = TRUE, add = TRUE)

# create head-specified wells
hw1 = headwell(xw = 300, yw = 100, hc = 6)
hw2 = headwell(xw = -200, yw = -100, xc = 0, yc = 0, hc = 7)

# create reference point element
rf = constant(x = -1000, y = 0, h = 8)

# create and solve model
m = aem(k, top, base, n, hw1, hw2, rf) 

# plot head contours
contours(m, xg, yg, col = 'dodgerblue', nlevels = 20)

# computed discharge of hw1
m$elements$hw1$parameter
#> [1] 292.5061

# computed discharge of hw2
element_discharge(m, name = 'hw2')
#>      hw2 
#> 779.7601

Reference point

When solving for unknown element parameters, a so-called reference point (or reference head) needs to be added to the model in order to have enough equations in the matrix formulation to solve for all unknowns. This is a constant value which is added to all hydraulic head values. It can be created by the constant() function which needs an \(x\) and \(y\) coordinate and the reference head \(h_c\). It is important to understand that this reference point is not an actual aquifer feature, but a mathematical necessity. Therefore, its value and location should be chosen is such a way that it does not influence the results too much. In general, the reference point should be never be placed inside the zone of interest (called the “near-field” in AEM lingo). It should be placed in the far-field and given a “neutral” head value as to not induce unrealistic gradients towards the near-field.

Line-sink

Line-sinks are elements with a line geometry which can be used to simulate groundwater flow to and from surface water features such as streams or canals. The parameter of a line-sink is its strength \(\sigma\), which is defined as the total flow to/from the line-sink per unit length of the line \([L^2/T]\) (positive going into the line-sink). The parameter is defined at a collocation point at the center of the line-sink.

Here, a single line-sink is created using linesink(). In addition to the start- and endpoints of the line-sink \(x_0-y_0\) and \(x_1-y_1\), the strength parameter needs to be specified. Optionally, a width can also be given, which for a strength-specified line-sink is only used in particle tracking to determine when the particles have reached the line.

Note that the locations of all elements in an aem object with line and point geometries (except the reference point) can be plotted using the plot.aem() method. Individual element locations can always be plotted by calling plot.element() on the respective element object.

# create line-sink
ls = linesink(x0 = -200, y0 = 200, x1 = 200, y1 = -200, sigma = 5)

# create and solve the model
m = aem(k, top, base, n, ls)

# plot head contours and line-sink geometry
contours(m, xg, yg, col = 'dodgerblue', nlevels = 20, drawlabels = FALSE)
plot(m, add = TRUE)

el = list(rf = rf) # list of elements
nls = 10           # number of line-sinks
xls = seq(-700, 700, length = nls + 1) # x-coordinates of line-sinks
yls = c(rep(c(-25, 25), nls), -25)     # y-coordinates of line-sinks
hc = 7             # stream level, m
res = 2            # hydraulic resistance of streambed, days
width = 3          # stream width, m

# create head-specified line-sinks
for(i in seq(nls)) {
  hl = headlinesink(x0 = xls[i], 
                    x1 = xls[i + 1], 
                    y0 = yls[i],
                    y1 = yls[i + 1],
                    hc = hc,
                    resistance = res,
                    width = width)
  el[[paste0('hls_', i)]] = hl
}

# create and solve model
m = aem(k, top, base, n, el, verbose = TRUE)
#> Solving analytic element model ... 
#>  Non-linear model with 11 elements and 11 unknowns 
#>  Entering outer iteration loop ... 
#>   Iteration 1 
#>   Iteration 2 
#>   Iteration 3 
#>   Iteration 4 
#>   Iteration 5 
#>   Iteration 6 
#>   Iteration 7 
#>   Iteration 8 
#>   Iteration 9 
#>   Iteration 10 
#> Model solved

# plot head contours
contours(m, xg, yg, col = 'dodgerblue', nlevels = 10)
plot(m, add = TRUE)

Area-sink

Area-sinks are circular elements that can be added at the top or bottom of the model. They provide a constant flux and can be used to simulate areal recharge or leakage through the aquifer base. The circular geometry is sufficient as long as the disc covers the entire area of interest and the heads inside are fixed by head-specified elements. The parameter is the specified flux value N \([L/T]\), defined at the center of the area-sink and applied uniformly across the entire disc.

Use of an area-sink introduces vertical flow components, besides the ones that occur due to the curvature of the water-table in a phreatic aquifer when variably saturated thickness is used in the aem() model. Note that the stream function does not give reasonable results inside the area covered by an area-sink, as discussed above.

The areasink() function requires the center coordinates of the circular area-sink, as well as the flux value N (positive is into the aquifer, note that this differs from other elements) and a radius R. The value of R is not important as long as its large enough so that the disc covers the modeled region. By default, the area-sink is placed at the aquifer top, which can be overridden by setting location = 'base'.

To illustrate to effect an area-sink providing areal recharge at the aquifer top, the example of the head-specified line-sinks above is used and an area-sink is added. Note the increased head values compared to the case without recharge.

# add area-sink to list of elements
el$as = areasink(x = 0, y = 0, N = 0.001, R = 1000) 

# create and solve the model
m = aem(k, top, base, n, el)

# plot head contours
contours(m, xg, yg, col = 'dodgerblue', nlevels = 10)
plot(m, add = TRUE)

# create head-specified area-sink
has = headareasink(x = 0, 
                   y = 200, 
                   h = 5, 
                   resistance = 1, 
                   R = 100)

# create and solve model
m = aem(k, top, base, n, w1, w2, rf, has)

# plot head contours and area-sink geometry
contours(m, xg, yg, col = 'dodgerblue', nlevels = 20)
plot(has, add = TRUE, col = adjustcolor('grey50', alpha = 0.5))

Uniform flow

An element can be created which models a specified, uniform background flow in the aquifer. In general, the regional gradient should be modeled by including head-specified surface water elements which will induce a background flow in the model. For certain applications however, it may be useful to specify this uniform background flow explicitly, e.g. when the regional system is not fully included in the model. This can be done by using uniformflow(). It requires a transmissivity value which remains constant in the entire domain (i.e. has a constant saturated thickness independent of the simulated heads), a gradient which is positive in the direction of flow, and an angle in degrees counterclockwise from the positive x-axis:

# create uniform flow element
uf = uniformflow(TR = k * (top - base), gradient = 0.001, angle = -30)

# create model
m = aem(k, top, base, n, uf, rf)

# plot head contours
contours(m, xg, yg, col = 'dodgerblue', nlevels = 20)

# add stream lines
contours(m, xg, yg, variable = 'streamfunction', 
         col = 'orange', nlevels = 20, add = TRUE)

Capture zone

A special case of particle tracking is the capture zone delineation for a well. This can be done using capzone(), a thin wrapper around tracelines() which places npar evenly distributed particles around the well’s radius and tracks them back in time using time steps of length dt until the requested capture time is reached. By plotting the tracelines, the capture zone can be be delineated. Initial particle locations are placed at the aquifer base (zstart argument). Different results may be obtained in models with vertical flow components by placing initial particle locations at different zstart values. By increasing the amount of particles or decreasing the time step length, a sharper delineation of the capture zone envelope may be obtained.

# 5-year capture zone of well A
cpA_5 = capzone(m, wA, time = 5 * 365, npar = 10) 

# plot head contours and element geometries
xg = seq(-800, 300, length = 100)
yg = seq(-600, 300, length = 100)

contours(m, xg, yg, col = 'dodgerblue', nlevels = 20)
plot(m, add = TRUE)

# plot capture zone output
plot(cpA_5, add = TRUE)