## H process: Theis solution (Pumping well)

This benchmark is available as a Jupyter notebook: Parabolic/LiquidFlow/AxiSymTheis/axisym_theis.ipynb.

#modules
from IPython.display import Image
import os
import pyvista as pv
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import exp1
import vtk
from vtk.util.numpy_support import vtk_to_numpy
import matplotlib.tri as tri
import time
#settings
path='./'
fig_dir = "./figures/"
prj_name = "axisym_theis"
prj_file = f"{prj_name}.prj"
pvd_name = "liquid_pcs"
vtu_name = "axisym_theis.vtu"
title = "H process: Theis solution (Pumping well)"
out_dir = os.environ.get('OGS_TESTRUNNER_OUT_DIR', '_out')
if not os.path.exists(out_dir):
os.makedirs(out_dir)
Image(filename = fig_dir + "ogs-jupyter-lab.png", width=150, height=100)

Image(filename = fig_dir + "h-tet-1.png", width=150, height=100)

## H process: Theis solution

Problem description

Theis’ problem examines the transient lowering of the water table induced by a pumping well. The assumptions required by the Theis solution are:

The aquifer

• is homogeneous, isotropic, confined, infinite in radial extent,
• has uniform thickness, horizontal piezometric surface.

The well

• is fully penetrating the entire aquifer thickness,
• has a constant pumping rate,
• well storage effects can be neglected,
• no other wells or long term changes in regional water levels.

Analytical solution

The analytical solution of the drawdown as a function of time and distance is expressed by $s(r,t) = h_0 - h(r,t) = \frac{Q}{4\pi T}W(u), \quad \mathrm{where}\quad u = \frac{r^2S}{4Tt}.$

where

• $$s$$ [$$L$$] is the drawdown or change in hydraulic head,
• $$h_0$$ is the constant initial hydraulic head,
• $$h$$ is the hydrauic head at distance $$r$$ at time $$t$$
• $$Q$$ [$$L^3T^{-1}$$] is the constant pumping (discharge) rate
• $$S$$ [$$-$$] is the aquifer storage coefficient (volume of water released per unit decrease in $$h$$ per unit area)
• $$T$$ [$$L^2T^{-1}$$] is the transmissivity (a measure of how much water is transported horizontally per unit time).

The Well Function, $$W(u)$$ is the exponential integral, $$E_1(u).$$ $$W(u)$$ is defined by an infinite series: $W(u) = - \gamma - \ln u + \sum_{k=1}^\infty \frac{(-1)^{k+1} u^k}{k \cdot k!}$ where

• $$\gamma=0.577215664$$ is the Euler-Mascheroni constant

For practical applications an approximation to the exponential integral is used often: $W(u) \approx -\gamma - \ln u$

This results in an expression for $$s(r,t)$$ known as the Jacob equation: $s(r,t) = -\frac{Q}{4\pi T}\left(\gamma + \ln u \right).$ For more details we refer to Srivastava and Guzman-Guzman (1998).

#source: https://scipython.com/blog/linear-and-non-linear-fitting-of-the-theis-equation/

def calc_u(r, S, T, t):
"""Calculate and return the dimensionless time parameter, u."""

return r**2 * S / 4 / T / t

def theis_drawdown(t, S, T, Q, r):
"""Calculate and return the drawdown s(r,t) for parameters S, T.

This version uses the Theis equation, s(r,t) = Q * W(u) / (4.pi.T),
where W(u) is the Well function for u = Sr^2 / (4Tt).
S is the aquifer storage coefficient,
T is the transmissivity (m2/day),
r is the distance from the well (m), and
Q is the pumping rate (m3/day).

"""

u = calc_u(r, S, T, t)
s_theis = Q/4/np.pi/T * exp1(u)
return s_theis

Q = 2000        # Pumping rate from well (m3/day)
r = 10          # Distance from well (m)

# Time grid, days.
t = np.array([1, 2, 4, 8, 12, 16, 20, 30, 40, 50, 60, 70, 80, 90, 100])

# Calculate some synthetic data to fit.
S, T = 0.0003, 1000
s = theis_drawdown(t, S, T, Q, r)

# Plot the data
titlestring = "Theis: Analytical solution"
plt.title(titlestring)
plt.plot(t, s, label='r = '+str(r)+' m')
plt.xlabel(r'$t\;/\;\mathrm{days}$')
plt.ylabel(r'$s\;/\;\mathrm{m}$')
plt.legend()
plt.grid()
plt.show()

# Recalculation from days in sec
Q = 0.016      # Pumping rate from well (m3/s)
t = 864000     # Time in s.

# Distance from well (m)
##r = np.array([0.5, 1, 2, 4, 8, 12, 16, 20, 25, 30, 35, 40])
r = np.arange(1,41,1)
##print(r)

# Calculate some synthetic data to fit.
S = 0.001
T = 9.2903e-4
u = calc_u(r, S, T, t)
s = theis_drawdown(t, S, T, Q, r)

# Plot the data
titlestring = "Theis: Analytical solution"
plt.title(titlestring)
plt.plot(r, s, label='t = '+str(t)+' days')
plt.xlabel(r'$r\;/\mathrm{m}$')
plt.ylabel(r'$hydraulic head\;/\;\mathrm{m}$')
plt.legend()
plt.grid()
plt.show()

Numerical solution

mesh = pv.read(vtu_name)
print("inspecting vtu-file")
mesh
inspecting vtu-file

UnstructuredGridInformation
N Cells354
N Points476
X Bounds3.048e-01, 3.048e+02
Y Bounds0.000e+00, 1.000e+00
Z Bounds0.000e+00, 0.000e+00
N Arrays1
NameFieldTypeN CompMinMax
OGS5_pressurePointsfloat6410.000e+001.247e+01
print("inspecting mesh and initial conditions")
#file
reader.Update()  # Needed because of GetScalarRange
pressure = data.GetPointData().GetArray("OGS5_pressure")
#points
points = data.GetPoints()
npts = points.GetNumberOfPoints()
x = vtk_to_numpy(points.GetData())
triang = tri.Triangulation(x[:,0], x[:,1])
#plt.triplot(triang, 'go-', lw=1.0)
plt.triplot(triang,lw=0.2)
plt.tricontour(triang, pressure, 16)
inspecting mesh and initial conditions

<matplotlib.tri.tricontour.TriContourSet at 0x7feec8541fd0>


Running OGS

#run ogs
t0 = time.time()
print("run ogs")
print(f"ogs {prj_file} > log.txt")
! ogs {prj_file} -o {out_dir} > {out_dir}/log.txt
tf = time.time()
print("computation time: ", round(tf - t0, 2), " s.")
run ogs
ogs axisym_theis.prj > log.txt
computation time:  1.85  s.


Spatial Profiles

import vtuIO
import numpy as np
import matplotlib.pyplot as plt
pvdfile = vtuIO.PVDIO(f"{out_dir}/{pvd_name}.pvd", dim=2)
xaxis =  [(i,0,0) for i in np.linspace(start=1.0, stop=40, num=40)]

r_x = np.array(xaxis)[:,0]
time = [8.64,86.4,1728,24192,172800,604800,864000]

pressure_xaxis_t = pvdfile.read_set_data(t, 'OGS5_pressure', data_type="point", pointsetarray=xaxis)

plt.plot(r_x, pressure_xaxis_t, 'x', label='OGS5, t = 1728 sec')

for t in time:
pressure_xaxis_t = pvdfile.read_set_data(t, 'pressure', data_type="point", pointsetarray=xaxis)
plt.plot(r_x, pressure_xaxis_t, label='t = '+str(t)+' sec')
titlestring = "Theis: Numerical solution"
plt.title(titlestring)
plt.xlabel(r'$r\;/\mathrm{m}$')
plt.ylabel(r'$hydraulic head\;/\;\mathrm{m}$')
plt.legend()
plt.grid()
#plt.savefig("theis.png")
plt.show()
WARNING: Default interpolation backend changed to VTK. This might result in
slight changes of interpolated values if defaults are/were used.


time = [864000]
pressure_xaxis_t = pvdfile.read_set_data(t, 'pressure', data_type="point", pointsetarray=xaxis)
#plot configuration
##plt.rcParams['figure.figsize'] = (16, 6)
##plt.rcParams['font.size'] = 12
##fig1, (ax1, ax2) = plt.subplots(1, 2)

fig, ax=plt.subplots(ncols=2, figsize=(12,4))
titlestring = "Theis: Comparison analytical and numerical solution"
ax[0].set_title(titlestring)
ax[0].set_xlim(0,40)
ax[0].plot(r_x, pressure_xaxis_t, 'x', label='numerical solution (ogs6)')
ax[0].plot(r, s, label='analytical solution')
ax[0].set_xlabel(r'$radius\;/\;\mathrm{m}$')
ax[0].set_ylabel(r'$hydraulic head\;/\;\mathrm{m}$')
ax[0].grid()
ax[0].legend()

##diff = np.setdiff1d(s,pressure_xaxis_t,assume_unique=False)
##print(diff)
titlestring = "Difference between analytical and numerical solutions"
caption = "Differences are due to different boundary conditions"
ax[1].set_title(titlestring)
ax[1].set_xlim(0,40)
ax[1].plot(r, s-pressure_xaxis_t, label='')
ax[1].set_xlabel(r'$radius\;/\;\mathrm{m}$')
ax[1].set_ylabel(r'$diff\;/\;\mathrm{m}$')
ax[1].grid()
ax[1].text(5,0.7,caption,ha='left')

##plt.savefig("theis-ana+num.png")
plt.show()

import time
print(time.ctime())
Thu Jan 26 18:56:55 2023