Go to the end to download the full example code
The extended Theis solution for truncated power laws
We provide an extended theis solution, that incorporates the effectes of a heterogeneous conductivity field following a truncated power law. In addition, it incorporates the assumptions of the general radial flow model and provides an arbitrary flow dimension.
In the following this extended solution is compared to the standard theis solution for well flow. You can nicely see, that the extended solution represents a transition between the theis solutions for the well- and farfield-conductivity.
Reference: (not yet published)
import numpy as np from matplotlib import pyplot as plt from anaflow import ext_theis_tpl, theis
We use three time steps: 10s, 10min, 10h
time_labels = ["10 s", "10 min", "10 h"] time = [10, 600, 36000] # 10s, 10min, 10h
Radius from the pumping well should be in [0, 4].
rad = np.geomspace(0.05, 4)
Parameters of heterogeneity, storage and pumping rate.
var = 0.5 # variance of the log-conductivity len_scale = 20.0 # upper bound for the length scale hurst = 0.5 # hurst coefficient KG = 1e-4 # the geometric mean of the conductivity KH = KG * np.exp(-var / 2) # the harmonic mean of the conductivity S = 1e-4 # storage rate = -1e-4 # pumping rate
Now let’s compare the extended Theis TPL solution to the classical solutions for the near and far field values of transmissivity.
head_KG = theis(time, rad, S, KG, rate) head_KH = theis(time, rad, S, KH, rate) head_ef = ext_theis_tpl( time=time, rad=rad, storage=S, cond_gmean=KG, len_scale=len_scale, hurst=hurst, var=var, rate=rate, ) time_ticks =  for i, step in enumerate(time): label_TG = "Theis($K_G$)" if i == 0 else None label_TH = "Theis($K_H$)" if i == 0 else None label_ef = "extended Theis TPL 2D" if i == 0 else None plt.plot( rad, head_KG[i], label=label_TG, color="C" + str(i), linestyle="--" ) plt.plot( rad, head_KH[i], label=label_TH, color="C" + str(i), linestyle=":" ) plt.plot(rad, head_ef[i], label=label_ef, color="C" + str(i)) time_ticks.append(head_ef[i][-1]) plt.xscale("log") plt.xlabel("r in [m]") plt.ylabel("h in [m]") plt.legend() ylim = plt.gca().get_ylim() plt.gca().set_xlim([rad, rad[-1]]) ax2 = plt.gca().twinx() ax2.set_yticks(time_ticks) ax2.set_yticklabels(time_labels) ax2.set_ylim(ylim) plt.tight_layout() plt.show()
Total running time of the script: ( 0 minutes 0.416 seconds)