Test Theis Storage

import numpy as np
from scipy.integrate import quad, solve_ivp
from scipy.special import exp1

Q = -100
npor = 0.3
k = 10
H = 10
T = k * H
Ss = 1e-4
S = Ss * H

def headtheis(r, t, T=T, S=S, Q=Q):
    u = r**2 * S / (4 * T * t)
    h = -Q / (4 * np.pi * T) * exp1(u)
    return h


def voltheis(r, t, T, S, Q):
    u = r**2 * S / (4 * T * t)
    h = -Q / (4 * np.pi * T) * exp1(u)
    vol = h * 2 * np.pi * r
    return vol * S

# demonstrate headtheis works correctly
quad(voltheis, a=1e-3, b=10000, args=(10, T, S, Q))  # gives 1000 m^3

(999.9999999944237, 1.4198832332112943e-06)

def volume(r, t, T=T, S=S, Q=Q):
    return quad(voltheis, a=1e-3, b=r, args=(10, T, S, Q))[0] + npor * np.pi * (r**2)

def vxytheis(t, xy):
    x, y = xy
    r = np.sqrt(x**2 + y**2)
    u = S * r**2 / (4 * T * t)
    Qr = -Q / (2 * np.pi) / r * np.exp(-u)
    vr = Qr / (H * npor)
    vx = vr * x / r
    vy = vr * y / r
    return np.array([vx, vy])

def vxytheisnew(t, xy):
    x, y = xy
    r = np.sqrt(x**2 + y**2)
    u = S * r**2 / (4 * T * t)
    Qr = -Q / (2 * np.pi) / r * np.exp(-u)
    vr = Qr / (H * (npor + S * headtheis(x, t)))
    vx = vr * x / r
    vy = vr * y / r
    return np.array([vx, vy])

# trace pathline for 10 days
path0 = solve_ivp(vxytheis, (1e-5, 10), y0=[1e-5, 0], method="DOP853")
R0 = path0.y[0, -1]
R0, np.pi * R0**2 * npor * H

(np.float64(10.300504380350029), np.float64(999.9726219155026))

path0 = solve_ivp(vxytheisnew, (1e-5, 10), y0=[1e-5, 0], method="DOP853")
R1 = path0.y[0, -1]
R1, np.pi * R1**2 * npor * H

(np.float64(10.28692090332354), np.float64(997.3369938953236))

headtheis(5, 10) * S

np.float64(0.0009076383332409223)

npor * H

3.0

R = path0.y[0, -1]
print("R, volume", R, np.pi * R**2 * npor * H)

R, volume 10.28692090332354 997.3369938953236

from scipy.integrate import quad

quad(headtheis, 1e-5, R, args=(100, T, S, Q))[0] + np.pi * R**2 * npor * H

np.float64(1009.0147425061984)

t0 = 0
t1 = 10


def vxytheis2(t, xy):
    x, y = xy
    r = np.sqrt(x**2 + y**2)
    u1 = S * r**2 / (4 * T * (t - t0))
    u2 = S * r**2 / (4 * T * (t - t1))
    Qr = -Q / (2 * np.pi) / r * np.exp(-u1)
    if t >= t1:
        Qr += 2 * Q / (2 * np.pi) / r * np.exp(-u2)
    vr = Qr / (H * npor)
    vx = vr * x / r
    vy = vr * y / r
    return np.array([vx, vy])

path1 = solve_ivp(vxytheis2, (10 + 1e-6, 20), y0=[R, 0], method="DOP853")
path1.y[0, -1]

np.float64(0.23642465667946888)

path1

  message: The solver successfully reached the end of the integration interval.
  success: True
   status: 0
        t: [ 1.000e+01  1.027e+01  1.060e+01  1.268e+01  1.700e+01
             1.878e+01  1.951e+01  1.981e+01  1.994e+01  1.999e+01
             2.000e+01]
        y: [[ 1.029e+01  1.016e+01 ...  4.482e-01  2.364e-01]
            [ 0.000e+00  0.000e+00 ...  0.000e+00  0.000e+00]]
      sol: None
 t_events: None
 y_events: None
     nfev: 194
     njev: 0
      nlu: 0

vxytheis2(12 + 1e-6, (10, 0))

array([-0.53039491, -0.        ])