Beside the implementations of solutions form literatur, AnaFlow provides some advanced features to pimp the solutions for the groundwater flow equation.

## 1. Self defined radial conductivity or transmissivity¶

All heterogeneous solutions of AnaFlow are derived by calculating an equivalent step function of a radial symmetric transmissivity resp. conductivity function.

The following code shows how to apply this workflow to a self defined transmissivity function. The function in use represents a linear transition from a local to a far field value of transmissivity within a given range.

The step function is calculated as the harmonic mean within given bounds, since the groundwater flow under a pumping condition is perpendicular to the different annular regions of transmissivity.

Reference: (not yet published)

import numpy as np
from matplotlib import pyplot as plt
import matplotlib.gridspec as gridspec
from anaflow.tools import specialrange_cut, annular_hmean, step_f

"""Conductivity with linear increase from K_well to K_far."""
return np.minimum(np.abs(rad) / len_scale, 1.0) * (K_far - K_well) + K_well

time_labels = ["10 s", "100 s", "1000 s"]
time = [10, 100, 1000]
S = 1e-4
K_well = 1e-5
K_far = 1e-4
len_scale = 5.0
rate = -1e-4
dim = 1.5

cut_off = len_scale
parts = 30
r_well = 0.0
r_bound = 50.0

# calculate a disk-distribution of "trans" by calculating harmonic means
R_part = specialrange_cut(r_well, r_bound, parts, cut_off)
K_part = annular_hmean(cond, R_part, ann_dim=dim, K_far=K_far, K_well=K_well, len_scale=len_scale)
S_part = np.full_like(K_part, S)

# plotting
gs = gridspec.GridSpec(2, 1, height_ratios=[1, 3])
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1], sharex=ax1)
time_ticks=[]
for i, step in enumerate(time):
label = "Transient" if i == 0 else None

ax1.set_yticks([K_well, K_far])
ax1.set_ylabel(r"$K$ in $[\frac{m}{s}]$")
plt.setp(ax1.get_xticklabels(), visible=False)
ax1.legend()
ax2.set_xlabel("r in [m]")
ax2.set_ylabel("h in [m]")
ax2.legend()
ax3 = ax2.twinx()
ax3.set_yticks(time_ticks)
ax3.set_yticklabels(time_labels)
ax3.set_ylim(ax2.get_ylim())

plt.tight_layout()
plt.show()


## 2. Accruing pumping rate¶

AnaFlow provides different representations for the pumping condition. One is an accruing pumping rate represented by the error function. This could be interpreted as that the water pump needs a certain time to reach its constant rate state.

import numpy as np
from scipy.special import erf
from matplotlib import pyplot as plt
import matplotlib.gridspec as gridspec
from anaflow import theis

time = np.geomspace(1, 600)

# Q(t) = Q * erf(t / a)
a = 120
lap_kwargs = {"cond": 4, "cond_kw": {"a": a}}
time=time,
storage=1e-4,
transmissivity=1e-4,
rate=-1e-4,
lap_kwargs=lap_kwargs,
)
time=time,
storage=1e-4,
transmissivity=1e-4,
rate=-1e-4,
)
gs = gridspec.GridSpec(2, 1, height_ratios=[1, 3])
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1], sharex=ax1)

ax2.plot(
time,
color="C" + str(i),
label="accruing Theis(r={})".format(step),
)
ax2.plot(
time,
color="C" + str(i),
label="constant Theis(r={})".format(step),
linestyle="--"
)
ax1.plot(time, 1e-4 * erf(time / a), label="accruing Q")
ax2.set_xlabel("t in [s]")
ax2.set_ylabel("h in [m]")
ax1.set_ylabel(r"|Q| in [$\frac{m^3}{s}$]")
ax1.legend()
ax2.legend()
plt.tight_layout()
plt.show()


## 3. Interval pumping¶

Another case often discussed in literatur is interval pumping, where the pumping is just done in a certain time frame.

Unfortunatly the Stehfest algorithm is not suitable for this kind of solution, which is demonstrated in the following script.

import numpy as np
from matplotlib import pyplot as plt
from anaflow import theis

time = np.linspace(10, 200)

# Q(t) = Q * characteristic([0, T])
lap_kwargs = {"cond": 3, "cond_kw": {"a": 100}}
time=time,
storage=1e-4,
transmissivity=1e-4,
rate=-1e-4,
lap_kwargs=lap_kwargs,
)