Contents

AnaFlow Quickstart

AnaFlow provides several analytical and semi-analytical solutions for the groundwater-flow equation.

Installation

The package can be installed via pip. On Windows you can install WinPython to get Python and pip running.

pip install anaflow

Provided Functions

The following functions are provided directly

• thiem Thiem solution for steady state pumping

• theis Theis solution for transient pumping

• ext_thiem_2d extended Thiem solution in 2D from Zech 2013

• ext_theis_2d extended Theis solution in 2D from Mueller 2015

• ext_thiem_3d extended Thiem solution in 3D from Zech 2013

• ext_theis_3d extended Theis solution in 3D from Mueller 2015

• neuman2004 transient solution from Neuman 2004

• neuman2004_steady steady solution from Neuman 2004

• grf “General Radial Flow” Model from Barker 1988

• ext_grf the transient extended GRF model

• ext_grf_steady the steady extended GRF model

• ext_thiem_tpl extended Thiem solution for truncated power laws

• ext_theis_tpl extended Theis solution for truncated power laws

• ext_thiem_tpl_3d extended Thiem solution in 3D for truncated power laws

• ext_theis_tpl_3d extended Theis solution in 3D for truncated power laws

Laplace Transformation

We provide routines to calculate the laplace-transformation as well as the inverse laplace-transformation of a given function

• get_lap Get the laplace transformation of a function

• get_lap_inv Get the inverse laplace transformation of a function

Requirements

• NumPy >= 1.14.5

• SciPy >= 1.1.0

• pentapy

License

MIT

AnaFlow Tutorial

In the following you will find several Tutorials on how to use AnaFlow to explore its whole beauty and power.

Gallery

The Theis solution

The Theis solution

The extended Theis solution in 2D

The extended Theis solution in 2D

The extended Theis solution in 3D

The extended Theis solution in 3D

The extended Theis solution for truncated power laws

The extended Theis solution for truncated power laws

The transient heterogeneous Neuman solution from 2004

The transient heterogeneous Neuman solution from 2004

extended Thiem 2D vs. steady solution for coarse graining transmissivity

extended Thiem 2D vs. steady solution for coarse graining transmissivity

extended Thiem 3D vs. steady solution for coarse graining conductivity

extended Thiem 3D vs. steady solution for coarse graining conductivity

extended Thiem 2D vs. steady solution for apparent transmissivity from Neuman

extended Thiem 2D vs. steady solution for apparent transmissivity from Neuman

extended Theis 2D vs. transient solution for apparent transmissivity from Neuman

extended Theis 2D vs. transient solution for apparent transmissivity from Neuman

Convergence of the extended Theis solutions for truncated power laws

Convergence of the extended Theis solutions for truncated power laws

Convergence of the general radial flow model (GRF)

Convergence of the general radial flow model (GRF)

Quasi steady convergence

Quasi steady convergence

Self defined radial conductivity or transmissivity

Self defined radial conductivity or transmissivity

Interval pumping

Interval pumping

Accruing pumping rate

Accruing pumping rate

The Theis solution

In the following the well known Theis function is called an plotted for three different time-steps.

Reference: Theis 1935

import numpy as np
from matplotlib import pyplot as plt

from anaflow import theis

time = [10, 100, 1000]
rad = np.geomspace(0.1, 10)

head = theis(time=time, rad=rad, storage=1e-4, transmissivity=1e-4, rate=-1e-4)

for i, step in enumerate(time):
    plt.plot(rad, head[i], label="Theis(t={})".format(step))

plt.legend()
plt.tight_layout()
plt.show()

The extended Theis solution in 2D

We provide an extended theis solution, that incorporates the effectes of a heterogeneous transmissivity field on a pumping test.

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 geometric- and harmonic-mean transmissivity.

Reference: Zech et. al. 2016

import numpy as np
from matplotlib import pyplot as plt

from anaflow import ext_theis_2d, 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-transmissivity
len_scale = 10.0  # correlation length of the log-transmissivity
TG = 1e-4  # the geometric mean of the transmissivity
TH = TG * np.exp(-var / 2.0)  # the harmonic mean of the transmissivity

S = 1e-4  # storativity
rate = -1e-4  # pumping rate

Now let’s compare the extended Theis solution to the classical solutions for the near and far field values of transmissivity.

head_TG = theis(time, rad, S, TG, rate)
head_TH = theis(time, rad, S, TH, rate)
head_ef = ext_theis_2d(time, rad, S, TG, var, len_scale, rate)
time_ticks = []
for i, step in enumerate(time):
    label_TG = "Theis($T_G$)" if i == 0 else None
    label_TH = "Theis($T_H$)" if i == 0 else None
    label_ef = "extended Theis" if i == 0 else None
    plt.plot(
        rad, head_TG[i], label=label_TG, color="C" + str(i), linestyle="--"
    )
    plt.plot(
        rad, head_TH[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.xlabel("r in [m]")
plt.ylabel("h in [m]")
plt.legend()
ylim = plt.gca().get_ylim()
plt.gca().set_xlim([0, 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()

The extended Theis solution in 3D

We provide an extended theis solution, that incorporates the effectes of a heterogeneous conductivity field on a pumping test. It also includes an anisotropy ratio of the horizontal and vertical length scales.

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 effective and harmonic-mean conductivity.

Reference: Müller 2015

import numpy as np
from matplotlib import pyplot as plt

from anaflow import ext_theis_3d, theis
from anaflow.tools.special import aniso

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, extend and pumping rate.

var = 0.5  # variance of the log-conductivity
len_scale = 10.0  # correlation length of the log-conductivity
anis = 0.75  # anisotropy ratio of the log-conductivity
KG = 1e-4  # the geometric mean of the conductivity
Kefu = KG * np.exp(
    var * (0.5 - aniso(anis))
)  # the effective conductivity for uniform flow
KH = KG * np.exp(-var / 2.0)  # the harmonic mean of the conductivity

S = 1e-4  # storage
L = 1.0  # vertical extend of the aquifer
rate = -1e-4  # pumping rate

Now let’s compare the extended Theis solution to the classical solutions for the near and far field values of transmissivity.

head_Kefu = theis(time, rad, S, Kefu * L, rate)
head_KH = theis(time, rad, S, KH * L, rate)
head_ef = ext_theis_3d(time, rad, S, KG, var, len_scale, anis, L, rate)
time_ticks = []
for i, step in enumerate(time):
    label_TG = "Theis($K_{efu}$)" if i == 0 else None
    label_TH = "Theis($K_H$)" if i == 0 else None
    label_ef = "extended Theis 3D" if i == 0 else None
    plt.plot(
        rad, head_Kefu[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.xlabel("r in [m]")
plt.ylabel("h in [m]")
plt.legend()
ylim = plt.gca().get_ylim()
plt.gca().set_xlim([0, 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()

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[0], 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()

The transient heterogeneous Neuman solution from 2004

We provide the transient pumping solution for the apparent transmissivity from Neuman 2004. This solution is build on the apparent transmissivity from Neuman 2004, which represents a mean drawdown in an ensemble of pumping tests in heterogeneous transmissivity fields following an exponential covariance.

In the following this 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: Neuman 2004

import numpy as np
from matplotlib import pyplot as plt

from anaflow import neuman2004, 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-transmissivity
len_scale = 10.0  # correlation length of the log-transmissivity
TG = 1e-4  # the geometric mean of the transmissivity
TH = TG * np.exp(-var / 2.0)  # the harmonic mean of the transmissivity

S = 1e-4  # storativity
rate = -1e-4  # pumping rate

Now let’s compare the apparent Neuman solution to the classical solutions for the near and far field values of transmissivity.

head_TG = theis(time, rad, S, TG, rate)
head_TH = theis(time, rad, S, TH, rate)
head_ef = neuman2004(
    time=time,
    rad=rad,
    trans_gmean=TG,
    var=var,
    len_scale=len_scale,
    storage=S,
    rate=rate,
)
time_ticks = []
for i, step in enumerate(time):
    label_TG = "Theis($T_G$)" if i == 0 else None
    label_TH = "Theis($T_H$)" if i == 0 else None
    label_ef = "transient Neuman [2004]" if i == 0 else None
    plt.plot(
        rad, head_TG[i], label=label_TG, color="C" + str(i), linestyle="--"
    )
    plt.plot(
        rad, head_TH[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[0], 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()

extended Thiem 2D vs. steady solution for coarse graining transmissivity

The extended Thiem 2D solutions is the analytical solution of the groundwater flow equation for the coarse graining transmissivity for pumping tests. Therefore the results should coincide.

References:

• Schneider & Attinger 2008

• Zech & Attinger 2016

import numpy as np
from matplotlib import pyplot as plt

from anaflow import ext_grf_steady, ext_thiem_2d
from anaflow.tools.coarse_graining import T_CG

rad = np.geomspace(0.05, 4)  # radius from the pumping well in [0, 4]
r_ref = 2.0  # reference radius

var = 0.5  # variance of the log-transmissivity
len_scale = 10.0  # correlation length of the log-transmissivity
TG = 1e-4  # the geometric mean of the transmissivity

rate = -1e-4  # pumping rate

head1 = ext_thiem_2d(rad, r_ref, TG, var, len_scale, rate)
head2 = ext_grf_steady(
    rad, r_ref, T_CG, rate=rate, trans_gmean=TG, var=var, len_scale=len_scale
)

plt.plot(rad, head1, label="Ext Thiem 2D")
plt.plot(rad, head2, label="grf(T_CG)", linestyle="--")

plt.xlabel("r in [m]")
plt.ylabel("h in [m]")
plt.legend()
plt.tight_layout()
plt.show()

extended Thiem 3D vs. steady solution for coarse graining conductivity

The extended Thiem 3D solutions is the analytical solution of the groundwater flow equation for the coarse graining conductivity for pumping tests. Therefore the results should coincide.

Reference: Zech et. al. 2012

import numpy as np
from matplotlib import pyplot as plt

from anaflow import ext_grf_steady, ext_thiem_3d
from anaflow.tools.coarse_graining import K_CG

rad = np.geomspace(0.05, 4)  # radius from the pumping well in [0, 4]
r_ref = 2.0  # reference radius

var = 0.5  # variance of the log-transmissivity
len_scale = 10.0  # correlation length of the log-transmissivity
KG = 1e-4  # the geometric mean of the transmissivity
anis = 0.7  # aniso ratio

rate = -1e-4  # pumping rate

head1 = ext_thiem_3d(rad, r_ref, KG, var, len_scale, anis, 1, rate)
head2 = ext_grf_steady(
    rad,
    r_ref,
    K_CG,
    rate=rate,
    cond_gmean=KG,
    var=var,
    len_scale=len_scale,
    anis=anis,
)

plt.plot(rad, head1, label="Ext Thiem 3D")
plt.plot(rad, head2, label="grf(K_CG)", linestyle="--")

plt.xlabel("r in [m]")
plt.ylabel("h in [m]")
plt.legend()
plt.tight_layout()
plt.show()

extended Thiem 2D vs. steady solution for apparent transmissivity from Neuman

Both, the extended Thiem and the Neuman solution, represent an effective steady drawdown in a heterogeneous aquifer. In both cases the heterogeneity is represented by two point statistics, characterized by mean, variance and length scale of the log transmissivity field. Therefore these approaches should lead to similar results.

References:

• Neuman 2004

• Zech & Attinger 2016

import numpy as np
from matplotlib import pyplot as plt

from anaflow import ext_thiem_2d, neuman2004_steady

rad = np.geomspace(0.05, 4)  # radius from the pumping well in [0, 4]
r_ref = 30.0  # reference radius

var = 0.5  # variance of the log-transmissivity
len_scale = 10.0  # correlation length of the log-transmissivity
TG = 1e-4  # the geometric mean of the transmissivity

rate = -1e-4  # pumping rate

head1 = ext_thiem_2d(rad, r_ref, TG, var, len_scale, rate)
head2 = neuman2004_steady(rad, r_ref, TG, var, len_scale, rate)

plt.plot(rad, head1, label="extended Thiem 2D")
plt.plot(rad, head2, label="Steady Neuman 2004", linestyle="--")

plt.xlabel("r in [m]")
plt.ylabel("h in [m]")
plt.legend()
plt.tight_layout()
plt.show()

extended Theis 2D vs. transient solution for apparent transmissivity from Neuman

Both, the extended Theis and the Neuman solution, represent an effective transient drawdown in a heterogeneous aquifer. In both cases the heterogeneity is represented by two point statistics, characterized by mean, variance and length scale of the log transmissivity field. Therefore these approaches should lead to similar results.

References:

• Neuman 2004

• Zech et. al. 2016

import numpy as np
from matplotlib import pyplot as plt

from anaflow import ext_theis_2d, neuman2004

time_labels = ["10 s", "10 min", "10 h"]
time = [10, 600, 36000]  # 10s, 10min, 10h

rad = np.geomspace(0.05, 4)  # radius from the pumping well in [0, 4]

TG = 1e-4  # the geometric mean of the transmissivity
var = 0.5  # correlation length of the log-transmissivity
len_scale = 10.0  # variance of the log-transmissivity

S = 1e-4  # storativity
rate = -1e-4  # pumping rate

head1 = ext_theis_2d(time, rad, S, TG, var, len_scale, rate)
head2 = neuman2004(time, rad, S, TG, var, len_scale, rate)
time_ticks = []
for i, step in enumerate(time):
    label1 = "extended Theis 2D" if i == 0 else None
    label2 = "Transient Neuman 2004" if i == 0 else None
    plt.plot(rad, head1[i], label=label1, color="C" + str(i))
    plt.plot(rad, head2[i], label=label2, color="C" + str(i), linestyle="--")
    time_ticks.append(head1[i][-1])

plt.title(
    "$T_G={}$, $\sigma^2={}$, $\ell={}$, $S={}$".format(TG, var, len_scale, S)
)
plt.xlabel("r in [m]")
plt.ylabel("h in [m]")
plt.legend()
ylim = plt.gca().get_ylim()
plt.gca().set_xlim([0, 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()

Convergence of the extended Theis solutions for truncated power laws

Here we set an outer boundary to the transient solution, so this condition coincides with the references head of the steady solution.

Reference: (not yet published)

import numpy as np
from matplotlib import pyplot as plt

from anaflow import ext_theis_tpl, ext_thiem_tpl

time = 1e4  # time point for steady state
rad = np.geomspace(0.1, 10)  # radius from the pumping well in [0, 4]
r_ref = 10.0  # reference radius

KG = 1e-4  # the geometric mean of the transmissivity
len_scale = 5.0  # correlation length of the log-transmissivity
hurst = 0.5  # hurst coefficient
var = 0.5  # variance of the log-transmissivity
dim = 1.5  # using a fractional dimension

S = 1e-4  # storativity
rate = -1e-4  # pumping rate

head1 = ext_thiem_tpl(
    rad, r_ref, KG, len_scale, hurst, var, dim=dim, rate=rate
)
head2 = ext_theis_tpl(
    time, rad, S, KG, len_scale, hurst, var, dim=dim, rate=rate, r_bound=r_ref
)

plt.plot(rad, head1, label="Ext Thiem TPL")
plt.plot(rad, head2, label="Ext Theis TPL (t={})".format(time), linestyle="--")

plt.xlabel("r in [m]")
plt.ylabel("h in [m]")
plt.legend()
plt.tight_layout()
plt.show()

Convergence of the general radial flow model (GRF)

The GRF model introduces an arbitrary flow dimension and was presented to analyze groundwater flow in rock formations. In the following we compare the bounded transient solution for late times, the unbounded quasi steady solution and the steady state.

Reference: Barker 1988

import numpy as np
from matplotlib import pyplot as plt

from anaflow import ext_grf, ext_grf_steady, grf

time = 1e4  # time point for steady state
rad = np.geomspace(0.1, 10)  # radius from the pumping well in [0, 4]
r_ref = 10.0  # reference radius

K = 1e-4  # the geometric mean of the transmissivity
dim = 1.5  # using a fractional dimension

rate = -1e-4  # pumping rate

head1 = ext_grf_steady(rad, r_ref, K, dim=dim, rate=rate)
head2 = ext_grf(time, rad, [1e-4], [K], [0, r_ref], dim=dim, rate=rate)
head3 = grf(time, rad, 1e-4, K, dim=dim, rate=rate)
head3 -= head3[-1]  # quasi-steady

plt.plot(rad, head1, label="Ext GRF steady")
plt.plot(rad, head2, label="Ext GRF (t={})".format(time), linestyle="--")
plt.plot(
    rad, head3, label="GRF quasi-steady (t={})".format(time), linestyle=":"
)

plt.xlabel("r in [m]")
plt.ylabel("h in [m]")
plt.legend()
plt.tight_layout()
plt.show()

Quasi steady convergence

The quasi steady is reached, when the radial shape of the drawdown in not changing anymore.

import numpy as np
from matplotlib import pyplot as plt

from anaflow import theis, thiem

time = [10, 100, 1000]
rad = np.geomspace(0.1, 10)
r_ref = 10.0

head_ref = theis(
    time,
    np.full_like(rad, r_ref),
    storage=1e-3,
    transmissivity=1e-4,
    rate=-1e-4,
)
head1 = (
    theis(time, rad, storage=1e-3, transmissivity=1e-4, rate=-1e-4) - head_ref
)
head2 = theis(
    time, rad, storage=1e-3, transmissivity=1e-4, rate=-1e-4, r_bound=r_ref
)
head3 = thiem(rad, r_ref, transmissivity=1e-4, rate=-1e-4)

for i, step in enumerate(time):
    label_1 = "Theis quasi steady" if i == 0 else None
    label_2 = "Theis bounded" if i == 0 else None
    plt.plot(rad, head1[i], label=label_1, color="C" + str(i), linestyle="--")
    plt.plot(rad, head2[i], label=label_2, color="C" + str(i))

plt.plot(rad, head3, label="Thiem", color="k", linestyle=":")

plt.xlabel("r in [m]")
plt.ylabel("h in [m]")
plt.legend()
plt.tight_layout()
plt.show()

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 matplotlib.gridspec as gridspec
import numpy as np
from matplotlib import pyplot as plt

from anaflow import ext_grf, ext_grf_steady
from anaflow.tools import annular_hmean, specialrange_cut, step_f


def cond(rad, K_far, K_well, len_scale):
    """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]
rad = np.geomspace(0.1, 6)

K_well = 1e-5
K_far = 1e-4
len_scale = 5.0
dim = 1.5

rate = -1e-4
S = 1e-4

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)
# calculate transient and steady heads
head1 = ext_grf(time, rad, S_part, K_part, R_part, dim=dim, rate=rate)
head2 = ext_grf_steady(
    rad,
    r_bound,
    cond,
    dim=dim,
    rate=-1e-4,
    K_far=K_far,
    K_well=K_well,
    len_scale=len_scale,
)

# 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
    ax2.plot(rad, head1[i], label=label, color="C" + str(i))
    time_ticks.append(head1[i][-1])

ax2.plot(rad, head2, label="Steady", color="k", linestyle=":")

rad_lin = np.linspace(rad[0], rad[-1], 1000)
ax1.plot(rad_lin, step_f(rad_lin, R_part, K_part), label="step Conductivity")
ax1.plot(
    rad_lin, cond(rad_lin, K_far, K_well, len_scale), label="Conductivity"
)
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()
ax2.set_xlim([0, rad[-1]])
ax3 = ax2.twinx()
ax3.set_yticks(time_ticks)
ax3.set_yticklabels(time_labels)
ax3.set_ylim(ax2.get_ylim())

plt.tight_layout()
plt.show()

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)
rad = [1, 5]

# Q(t) = Q * characteristic([0, a])
lap_kwargs = {"cond": 3, "cond_kw": {"a": 100}}

head = theis(
    time=time,
    rad=rad,
    storage=1e-4,
    transmissivity=1e-4,
    rate=-1e-4,
    lap_kwargs=lap_kwargs,
)

for i, step in enumerate(rad):
    plt.plot(time, head[:, i], label="Theis(r={})".format(step))

plt.title("The Stehfest algorithm is not suitable for this!")
plt.legend()
plt.tight_layout()
plt.show()

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 matplotlib.gridspec as gridspec
import numpy as np
from matplotlib import pyplot as plt
from scipy.special import erf

from anaflow import theis

time = np.geomspace(1, 600)
rad = [1, 5]

# Q(t) = Q * erf(t / a)
a = 120
lap_kwargs = {"cond": 4, "cond_kw": {"a": a}}

head1 = theis(
    time=time,
    rad=rad,
    storage=1e-4,
    transmissivity=1e-4,
    rate=-1e-4,
    lap_kwargs=lap_kwargs,
)
head2 = theis(
    time=time,
    rad=rad,
    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)

for i, step in enumerate(rad):
    ax2.plot(
        time,
        head1[:, i],
        color="C" + str(i),
        label="accruing Theis(r={})".format(step),
    )
    ax2.plot(
        time,
        head2[:, i],
        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()

AnaFlow API

Purpose

Anaflow provides several analytical and semi-analytical solutions for the groundwater-flow-equation.

Subpackages

flow	Anaflow subpackage providing flow-solutions for the groundwater flow equation.
tools	Anaflow subpackage providing miscellaneous tools.

Solutions

Homogeneous

Solutions for homogeneous aquifers

thiem(rad, r_ref, transmissivity[, rate, h_ref])	The Thiem solution.
theis(time, rad, storage, transmissivity[, ...])	The Theis solution.
grf(time, rad, storage, conductivity[, dim, ...])	The general radial flow (GRF) model for a pumping test.

Heterogeneous

Solutions for heterogeneous aquifers

ext_thiem_2d(rad, r_ref, trans_gmean, var, ...)	The extended Thiem solution in 2D.
ext_thiem_3d(rad, r_ref, cond_gmean, var, ...)	The extended Thiem solution in 3D.
ext_thiem_tpl(rad, r_ref, cond_gmean, ...[, ...])	The extended Thiem solution for truncated power-law fields.
ext_thiem_tpl_3d(rad, r_ref, cond_gmean, ...)	The extended Theis solution for truncated power-law fields in 3D.
ext_theis_2d(time, rad, storage, ...[, ...])	The extended Theis solution in 2D.
ext_theis_3d(time, rad, storage, cond_gmean, ...)	The extended Theis solution in 3D.
ext_theis_tpl(time, rad, storage, ...[, ...])	The extended Theis solution for truncated power-law fields.
ext_thiem_tpl_3d(rad, r_ref, cond_gmean, ...)	The extended Theis solution for truncated power-law fields in 3D.
neuman2004(time, rad, storage, trans_gmean, ...)	The transient solution for the apparent transmissivity from [Neuman2004].
neuman2004_steady(rad, r_ref, trans_gmean, ...)	The steady solution for the apparent transmissivity from [Neuman2004].

Extended GRF

The extended general radial flow model.

ext_grf(time, rad, S_part, K_part, R_part[, ...])	The extended "General radial flow" model for transient flow.
ext_grf_steady(rad, r_ref, conductivity[, ...])	The extended "General radial flow" model for steady flow.

Laplace

Helping functions related to the laplace-transformation

get_lap(func[, arg_dict])	Callable Laplace transform.
get_lap_inv(func[, method, method_dict, ...])	Callable Laplace inversion.

Tools

Helping functions.

step_f(rad, r_part, f_part)	Callalbe step function.
specialrange(val_min, val_max, steps[, typ])	Calculation of special point ranges.
specialrange_cut(val_min, val_max, steps[, ...])	Calculation of special point ranges.

Changelog

All notable changes to AnaFlow will be documented in this file.

[1.1.0] - 2023-04

See #11

Enhancements

• move to src/ base package structure

• drop py36 support

• added archive support

• simplify documentation

1.0.1 - 2020-04-02

Bugfixes

• ModuleNotFoundError not present in py35

• np.asscalar deprecated, use array.item()

• CHANGELOG.md links updated

1.0.0 - 2020-03-22

Enhancements

• new TPL Solution

• new tools sub-module

• using pentapy to solve LES in laplace space

• solution for aparent transmissivity from neuman 2004

• added extended GRF model

• convenient functions for (inverse-)laplace transformation

Bugfixes

• lat_ext was ignored by ext_theis_3d

Changes

• py2.7 support dropped

0.4.0 - 2019-03-07

Enhancements

• the output for transient tests now preserves the shapes of time and rad (better for plotting in 3D)

• the grf model is now the default way of calculating pumping tests in laplace space

• the grf_laplace routine was optimized to estimate the radius of the cone of depression

• the grf_laplace uses now the pentapy solver, so we get rid of the umf_pack dependency

• grf_model and grf_disk are now part of the standard routines

Changes

• the input for transient tests changed from “rad, time” to “time, rad” as order of input (in line with output format)

Bugfixes

• multiple minor bugfixes

0.3.0 - 2019-02-28

Enhancements

• GRF model added

• new documetation

• added examples

• code restructured

Changes

Bugfixes

0.2.4 - 2018-04-26

Enhancements

• Released on PyPI

0.1.0 - 2018-01-05

First release of AnaFlow. Containing:

• thiem - Thiem solution for steady state pumping

• theis - Theis solution for transient pumping

• ext_thiem2D - extended Thiem solution in 2D

• ext_theis2D - extended Theis solution in 2D

• ext_thiem3D - extended Thiem solution in 3D

• ext_theis3D - extended Theis solution in 3D

• diskmodel - Solution for a diskmodel

• lap_transgwflow_cyl - Solution for a diskmodel in laplace inversion