# Potential flow


## Potencial functions

**Uniform flow**:

$w = U\exp{(-i\alpha)}z$

**Line source**:

$w = \dfrac{\dot{Q}}{2\pi}\ln z$

## Velocity fields

$\dfrac{dw}{dz} = u + iv$

**Uniform flow:**

$\dfrac{dw}{dz} = U\left(\cos{\alpha}-i\sin{\alpha}\right)$   

**Line source:**

$\dfrac{dw}{dz} = \dfrac{\dot{Q}}{2\pi}\dfrac{1}{z}$

## Combine potentials 

Consider a potential flow with a sink at $(0,0)$ and a source at $(r,0)$, under the influence of a horizontal uniform flow going from left to right ($\alpha = \pi$). The potential that describes that system is given by:

$w = w_{\rm uniform} + w_{\rm source} + w_{\rm sink}$

Hence, 

$w = U\exp{(-i\alpha)}z + \dfrac{\dot{Q_{\rm in}}}{2\pi}\ln (z-r) - \dfrac{\dot{Q_{\rm out}}}{2\pi}\ln z$

$w = -Uz + \dfrac{\dot{Q_{\rm in}}}{2\pi}\ln (z-r) - \dfrac{\dot{Q_{\rm out}}}{2\pi}\ln z$

The uniform flow component will represent the regional flow, which can be described using Darcy's law as $U=-K\dot{I}$, where $K$ is the hydraulic conductivity (m/s) and $\dot{I}$ is the water table gradient (m/m). The source and sink wells represent injection and extraction points, characterized by flow rates per unit depth $\dot{Q_{\rm in}}$ and $\dot{Q_{\rm out}}$ (m²/s). For an aquifer of depth $H$, the volumetric flow rate $Q$ is $\dot{Q}=Q/H$. Now, the extraction rate can be expressed as a proportion to the injection rate as $Q_{\rm out} = fQ_{\rm in}$. With all this replacements, the potential function can be rewritten as:

$w = K\dot{I}z + \dfrac{Q_{\rm in}}{2\pi H} \left( \ln(z-r) - f \ln z \right)$

Notice that in our case, the gradient $\dot{I}$ must be negative so the regional flow goes indeed from right to left. To avoid confusion, let's rewrite it in terms of the abs value of the gradient $I$:

$w = -KIz + \dfrac{Q_{\rm in}}{2\pi H} \left( \ln(z-r) - f \ln z \right)$

## Velocity field 

$\dfrac{dw}{dz} = \dfrac{d}{dz}\left(-KIz + \dfrac{Q_{\rm in}}{2\pi H} \left( \ln(z-r) - f \ln z \right)\right)$

$\dfrac{dw}{dz} = -KI + \dfrac{Q_{\rm in}}{2\pi H} \dfrac{d}{dz}\left( \ln(z-r) - f \ln z \right)$

$\dfrac{dw}{dz} = -KI + \dfrac{Q_{\rm in}}{2\pi H} \left(\dfrac{1}{z-r} - \dfrac{f}{z} \right)$

## Velocity at the midpoint

We'll pick the velocity in the midpoint between the extraction and injection points as the characteristic velocity of the system. This means that we want to evaluate the velocity field at $z = r/2 + 0i$

\begin{equation}
\begin{array}{rl}
    \dfrac{dw}{dz}\bigg\rvert_{z=r/2+0i} =& -KI + \dfrac{Q_{\rm in}}{2\pi H} \left(\dfrac{1}{r/2-r} - \dfrac{f}{r/2} \right)\\
    =& -KI - \dfrac{Q_{\rm in}}{\pi H} \left(\dfrac{1+f}{r}\right)
\end{array}
\end{equation}

Notice that this velocity can be split into the contribution of the regional flow and the contribution from the wells:

\begin{equation}
\begin{array}{rl}
    w'_{\rm regional} =& -KI \\
    w'_{\rm wells} =& -\dfrac{1}{\pi}\dfrac{Q_{\rm in} (1+f)}{r H}
\end{array}
\end{equation}


### Non-dimensional flow number ($\mathcal{F}_L$)

A comparison between the component in the velocity due regional flow and the component due the wells as a ratio between them two. 

\begin{equation}
\begin{array}{rl}
\mathcal{F}_L =& \dfrac{\text{Velocity due regional flow}}{\text{Velocity due well couple}} \\
=& \dfrac{w'_{\rm regional}}{w'_{\rm wells}} \\
=& \dfrac{-KI}{-\dfrac{Q_{\rm in}}{2\pi H} \left(\dfrac{1+f}{r}\right)} \\
=& \dfrac{\pi K I H r}{Q_{in}(1+f)}
\end{array}
\end{equation}

If $\mathcal{F}_L \to 0$, the flow induced by the wells rule the system, whereas if $\mathcal{F}_L \to \infty$, the regional flow rule the system.

## Minimum velocity along a streamline

The minimum velocity is found at:

\begin{equation}
\begin{array}{rl}
    \dfrac{d}{dx}\dfrac{dw}{dz}\bigg\rvert_{z=x+0i} =& 0 \\
    \dfrac{d}{dx}\left(-KI + \dfrac{Q_{\rm in}}{2\pi H} \left(\dfrac{1}{x-r} - \dfrac{f}{x} \right)\right) =& 0\\
    \dfrac{Q_{\rm in}}{2 \pi H} \dfrac{d}{dx}\left(\dfrac{1}{x-r} - \dfrac{f}{x} \right) =& 0\\
    \dfrac{Q_{\rm in}}{2 \pi H} \left(\dfrac{f}{x^2} - \dfrac{1}{\left(x-r\right)^2}\right) =& 0\\
    \\
    \dfrac{f}{x^2} =& \dfrac{1}{\left(x-r\right)^2}\\
    x^2 =& f \left( x-r \right) ^2\\
    x =& \dfrac{\sqrt{f}}{\sqrt{f}\pm1}r
    \\
    x < r \Rightarrow x =& \dfrac{\sqrt{f}}{\sqrt{f}+1}r
\end{array}
\end{equation}

Defining $m = \dfrac{\sqrt{f}}{\sqrt{f}+1}$, evaluating at $z = mr + 0i$ gives the minimum velocity.


\begin{equation}
\begin{array}{rl}
    \dfrac{dw}{dz}\bigg\rvert_{z=mr+0i} =& -KI + \dfrac{Q_{\rm in}}{2\pi H} \left(\dfrac{1}{mr-r} - \dfrac{f}{mr} \right)\\
    =& -KI + \dfrac{Q_{\rm in}}{2 \pi H} \left(\dfrac{m-fm+f}{m(m-1)}\right)\\
    =& -KI - \dfrac{Q_{\rm in}}{2 \pi H} \left(\dfrac{\left(\sqrt{f}+f\right)\left(\sqrt{f}+1\right)}{\sqrt{f}}\right)\\
    =& -KI - \dfrac{Q_{\rm in}}{2 \pi H} \left( \left( 1+ \sqrt{f}\right)\left( 1 + \dfrac{\sqrt{f}}{f} \right) \right)\\
    =& -KI - \dfrac{ Q_{\rm in}\left( 1 + \tfrac{\sqrt{f}}{f} \right) } {\pi H r}
\end{array}\\
\end{equation}

This velocity can also be split into the contribution of the regional flow and the contribution from the wells:

\begin{equation}
\begin{array}{rl}
    w'_{\rm regional} =& -KI \\
    w'_{\rm wells} =& -\dfrac{Q_{\rm in} }{\pi r H}\left(1+\tfrac{\sqrt{f}}{f}\right)
\end{array}
\end{equation}

## Mean velocity along a streamline

Consider the streamline that connects the sink and source points. It cannot be measured at those points because those are poles in the system, hence, some separation $\delta$ from those points has to be considered:  the well radius should be the straighforward election so we'll adopt it and assume that both injection and extraction wells have the same diameter.

Now, with $z = x+iy$, the streamline exists at $y=0$ and $\delta<x<r-\delta$, hence, the mean velocity can be derived as follows:

\begin{equation}
\begin{array}{rl}
\bar{u} =& \dfrac{1}{r-2\delta} \int_{\delta}^{r-\delta}\dfrac{dw}{dz}dx\\
=& \dfrac{1}{r-2\delta} \int_{\delta}^{r-\delta} \left(-KI + \dfrac{Q_{\rm in}}{2\pi H} \left(\dfrac{1}{z-r} - \dfrac{f}{z} \right)\right)dx\\
=& -KI + \dfrac{Q_{\rm in}}{2\pi (r-2\delta) H}\int_{\delta}^{r-\delta}\left(\dfrac{1}{x-r} - \dfrac{f}{x} \right)dx\\
=& -KI + \dfrac{Q_{\rm in}}{2\pi H}\left(\ln({x-r}) - f\ln{x} \right)\bigg\rvert_{x=\delta}^{x=r-\delta}\\
=& -KI + \dfrac{Q_{\rm in}}{2\pi H}\left[ (1+f)(\ln(\delta) -\ln(r-\delta)  \right]\\
=& -KI - \dfrac{Q_{\rm in}(1+f)}{2\pi H}\ln\left(\dfrac{r-\delta}{\delta}\right)\\
\end{array}
\end{equation}

Notice that this velocity can be split into the contribution of the regional flow and the contribution from the wells:

\begin{equation}
\begin{array}{rl}
    w'_{\rm regional} =& -KI \\
    w'_{\rm wells} =& -\dfrac{Q_{\rm in}(1+f)}{2\pi H}\ln\left(\dfrac{r-\delta}{\delta}\right)
\end{array}
\end{equation}

### Non-dimensional flow number ($\mathcal{F}_L$)

A comparison between the component in the velocity due regional flow and the component due the wells as a ratio between them two. 

\begin{equation}
\begin{array}{rl}
\mathcal{F}_L =& \dfrac{\text{Velocity due regional flow}}{\text{Velocity due well couple}} \\
=& \dfrac{w'_{\rm regional}}{w'_{\rm wells}} \\
=& \dfrac{-KI}{-\dfrac{Q_{\rm in}(1+f)}{2\pi H}\ln\left(\dfrac{r-\delta}{\delta}\right)} \\
=& \dfrac{2 \pi K I H (r-2\delta)}{Q_{in}(1+f)\ln\left(\frac{r-\delta}{\delta}\right)}
\end{array}
\end{equation}

If $\mathcal{F}_L \to 0$, the flow induced by the wells rule the system, whereas if $\mathcal{F}_L \to \infty$, the regional flow rule the system.

## Mass flow through an infinite line

The previous choice of characteristic velocity seems arbitrary in the sense it overestimates the effect of the potential induced by the wells on the system. This is because the streamline along the x-axis has the highest velocity magnitude (the sharpest potential gradient), if we were to pick a point elsewhere, the veloicty form the wells would be less. But in the same sense, for $f=1$, this point represents the lowest velocity along that streamline. 

To avoid this apparent bias, let's consider the totality of the flow that crosses an infinite vertical line at $z=r/2+iy$. This won't give us a useful insight on the regional flow as that mass flow will diverge, but it will give us a point of comparison to the wells contribution to the total mass balance in the system. 

So, the velocity field along that vertical line is 

\begin{equation}
\begin{array}{rl}
\dfrac{dw}{dz}\bigg\rvert_{z=\tfrac{r}{2}+iy} =& -KI - \dfrac{Q_{\rm in}}{2\pi H} \left(\dfrac{1}{z-r} - \dfrac{f}{z} \right)\bigg\rvert_{z=\tfrac{r}{2}+iy}\\
=& - KI + \dfrac{Q_{\rm in}}{2\pi H} \left(\dfrac{1}{\tfrac{r}{2}+iy-r} - \dfrac{f}{\tfrac{r}{2}+iy} \right)\\
=& -KI - \dfrac{Q_{\rm in}}{2\pi H} \left(\dfrac{\tfrac{r}{2}(1+f) + iy(1-f)}{\tfrac{r^2}{4}-y^2}\right)
\end{array}
\end{equation}

The mass flow through that vertical line will require only the horizontal component of the velocity field $u$ at that point, meaning that we are only intereste in the real component of the function above:

$u = -KI - \dfrac{Q_{\rm in}}{2\pi H} \left(\dfrac{\tfrac{r}{2}(1+f)}{\tfrac{r^2}{4}-y^2}\right) = -KI - \dfrac{Q_{\rm in}r(1+f)}{4\pi H} \left(\dfrac{1}{\tfrac{r^2}{4}-y^2}\right) $

The mass flow will be described as:

\begin{equation}
\begin{array}{rl}
    \int_{-\infty}^{\infty} u \, dy =& \int_{-\infty}^{\infty}\left[ -KI - \dfrac{Q_{\rm in}r(1+f)}{4\pi H} \left(\dfrac{1}{\tfrac{r^2}{4}-y^2}\right)\right] \, dy \\
    =& -\underbrace{\int_{-\infty}^{\infty}KI \, dy}_{\text{Diverges}} - \underbrace{\dfrac{Q_{\rm in}r(1+f)}{4\pi H} \int_{-\infty}^{\infty}\left(\dfrac{1}{\tfrac{r^2}{4}-y^2}\right) \, dy}_{\text{Converges}}\\
    =& -\int_{-\infty}^{\infty}KI \, dy - \dfrac{Q_{\rm in}r(1+f)}{4\pi H} \left[ \dfrac{2}{r} \tan^{-1}{\left(\dfrac{2y}{r}\right)}\right]\bigg\rvert_{y \to -\infty}^{y \to \infty}\\
    =& -\int_{-\infty}^{\infty}KI \, dy - \dfrac{Q_{\rm in}r(1+f)}{4\pi H} \left[ \dfrac{2\pi}{r} \right] \\
    =& -\int_{-\infty}^{\infty}KI \, dy - \dfrac{Q_{\rm in}(1+f)}{2 H}
\end{array}
\end{equation}

This end term, we can divide it by the lenght of the line $l_r$ in order to obtain a velocity, hence, the contribution from the wells are:

\begin{equation}
\begin{array}{rl}
    w'_{\rm regional} =& -\dfrac{1}{l_r}\int_{-\infty}^{\infty}KI \, dy \\
    w'_{\rm wells} =& -\dfrac{1}{2l_r}\dfrac{Q_{\rm in}(1+f)}{rH}
\end{array}
\end{equation}



### Non-dimensional flow number ($\mathcal{F}_L$)

In this case, this comparison don't make much sense as one of the integrals diverged. 

\begin{equation}
\begin{array}{rl}
\mathcal{F}_L =& \dfrac{\text{Velocity due regional flow}}{\text{Velocity due well couple}} \\
=& \dfrac{w'_{\rm regional}}{w'_{\rm wells}} \\
=& \dfrac{-\dfrac{1}{l_r}\int_{-\infty}^{\infty}KI \, dy }{-\dfrac{1}{2l_r}\dfrac{Q_{\rm in}(1+f)}{rH}} \\
=& \dfrac{2rH \int_{-\infty}^{\infty}KI \, dy}{Q_{in}(1+f)}
\end{array}
\end{equation}

We could infer that in this approach, the regional flow will always dwarf the potential induced by the pair of wells. But this is not useful for our purpose :(

## Mass flow through an finite line

In order to obtain a convergent integral for the regional component, let's pick a vertical line that spans from $r/2(1-i)$ to $r/2(1+i)$, that way, the mass flow will be described as:

\begin{equation}
\begin{array}{rl}
    \int_{-r/2}^{r/2} u \, dy =& \int_{-r/2}^{r/2}\left[ -KI - \dfrac{Q_{\rm in}r(1+f)}{4\pi H} \left(\dfrac{1}{\tfrac{r^2}{4}-y^2}\right)\right] \, dy \\
    =& \int_{-r/2}^{r/2}KI \, dy - \dfrac{Q_{\rm in}r(1+f)}{4\pi H} \int_{-r/2}^{r/2}\left(\dfrac{1}{\tfrac{r^2}{4}-y^2}\right) \, dy\\
    =& -KIy\bigg\rvert_{y = -r/2}^{y = r/2} - \dfrac{Q_{\rm in}r(1+f)}{4\pi H} \left[ \dfrac{2}{r} \tan^{-1}{\left(\dfrac{2y}{r}\right)}\right]\bigg\rvert_{y = -r/2}^{y = r/2}\\
    =& -KIr - \dfrac{Q_{\rm in}r(1+f)}{4\pi H} \left[ \dfrac{\pi}{r} \right] \\
    =& -KIr - \dfrac{Q_{\rm in}(1+f)}{4 H}
\end{array}
\end{equation}

Jusy as before, this end term can be divided by $r$ in order to obtain a velocity, hence, both contributions are:

\begin{equation}
\begin{array}{rl}
    w'_{\rm regional} =& -KI \\
    w'_{\rm wells} =& -\dfrac{1}{4}\dfrac{Q_{\rm in}(1+f)}{rH}
\end{array}
\end{equation}


### Non-dimensional flow number (2) ($\mathcal{F}_L$)

From this analysis we could also draw a comparison between the two contributors of the system:

\begin{equation}
\begin{array}{rl}
\mathcal{F}_L =& \dfrac{\text{Velocity due regional flow}}{\text{Velocity due well couple}}\\
=& \dfrac{w'_{\rm regional}}{w'_{wells}}\\
=& \dfrac{-KI}{-\dfrac{Q_{\rm in}}{4 H} \left(\dfrac{1+f}{r}\right)} \\
=& \dfrac{4 K I H r}{Q_{in}(1+f)}
\end{array}
\end{equation}

There is not much of a difference in this "more conservative" approach, and it holds that if $\mathcal{F}_L \to 0$, the flow induced by the wells rule the system, whereas if $\mathcal{F}_L \to \infty$, the regional flow rule the system. We well conserve the first non-dimensional flow number definition, i.e., 

\begin{equation}
\mathcal{F}_L = \dfrac{4 K I H r}{Q_{in}(1+f)}
\end{equation}

## Characteristic velocity ($q_c$)

We found out before an expression for the velocity at the midpoint between the well couple. This will be used as our characteristic velocity, which can be rewritten in terms of the non-dimensional flow number:

\begin{equation}
\begin{array}{rl}
    q_c =& \dfrac{1}{r}\int_{-r/2}^{r/2}\dfrac{dw}{dz}dy\\
    =& -KI - \dfrac{Q_{\rm in}(1+f)}{4 r H}\\
    =& -KI \left( 1 + \dfrac{1}{\mathcal{F}_L} \right)\\
    =& - KI \left( \dfrac{\mathcal{F}_L}{\mathcal{F}_L + 1} \right)
\end{array}
\end{equation}

The actual pore-water velocity $u_c$ is found by dividing by the porosity $\theta$

\begin{equation}
\begin{array}{rl}
    u_c =& \dfrac{q_c}{\theta}\\
    =& - \dfrac{KI}{\theta} \left( \dfrac{\mathcal{F}_L}{\mathcal{F}_L + 1} \right)
\end{array}
\end{equation}

## Flow time scale ($\tau$)
Picking the half-distance between the two wells as our characteristic lenght $r/2$, a characteristic time $\tau$ can be derived from $u_c$

\begin{equation}
\begin{array}{rl}
    \tau =& \dfrac{r}{|u_c|}\\
    =& \dfrac{r\theta}{KI} \left( \dfrac{\mathcal{F}_L + 1}{\mathcal{F}_L} \right)\\
\end{array}
\end{equation}

## Decay rate and timescale comparison

Considering a tracer $C$ that is injected in the source point. This tracer follows a first-order decay reaction with a rate $\lambda$. 

$\dfrac{dC}{dt} = -\lambda C$

With an initial condition $C(0) = C_0$, the concentration of that tracer as a function of space and time is,

$C(t) = C_{0}\exp{-\lambda t}$

$C_\tau = C_0 \exp{-\lambda \tau} =  C_0 \exp{\left(-\lambda\dfrac{r}{|u_c|}\right)}$

Hence, for our characteristic time ($t=\tau$), the relative concentration of that tracer will be:

$C_\tau = C_0\exp{ \left( - \lambda \dfrac{r \theta (\mathcal{F}_L + 1)}{KI \mathcal{F}_L} \right) }$

## Dilution effect

Assuming that perfect mixing has between the injection point and the distance between the wells, we could calculate a relative concentration that takes into account this mass balance. Notice that this analysis cannot be done directly from our previous assumptions because there is no mixing in potential flow. However, we can keep the characteristic velocity.

\begin{equation}
\begin{array}{rl}
    C_0 u_c \Delta y \Delta z \theta=& C_{\rm in} Q_{\rm in} + C_{\rm reg} Q_{\rm reg}\\
    C_0 u_c \Delta y \Delta z \theta=& C_{\rm in} Q_{\rm in}\\
    C_0 =& \dfrac{Q_{\rm in}}{u_c \Delta y \Delta z \theta} C_{in} \\
    =& \dfrac{Q_{\rm in}}{u_c \Delta y \Delta z \theta} C_{\rm in}\\
    =& \dfrac{Q_{\rm in}}{- \dfrac{KI}{\theta} \left( \dfrac{\mathcal{F}_L}{\mathcal{F}_L + 1} \right) \Delta y \Delta z \theta} C_{\rm in}\\
    =& \dfrac{Q_{\rm in}(\mathcal{F}_L + 1)}{KI \mathcal{F}_L \Delta y \Delta z} C_{\rm in}
\end{array}
\end{equation}

$\Delta y$ and $\Delta z$ are just lenghts around the arbitrary control volume. For our 2D case, $\Delta z = H$ and $\Delta y = $ size of the element in PFLOTRAN

$C_0 =
\dfrac{Q_{\rm in}(\mathcal{F}_L + 1)}{KI \mathcal{F}_L \Delta y H} C_{\rm in}$

### Which process is more dominant?

The effects of each process is more enhanced dependind of the system. For example, if $I$ is big, our characteristic time will be small, then, decay won't be as important as dilution could be in the aquifer. We are interested in finding which set of parameters maximize $C_\tau$, or minimize the log-reductions of the tracer concentration.

\begin{equation}
\begin{array}{rl}
C_\tau =& C_0 \exp{\left(-\lambda\dfrac{r}{|u_c|}\right)}\\
C_\tau =& \dfrac{C_{\rm in}Q_{\rm in}}{u_c \Delta y H \theta} \exp{\left(-\lambda\dfrac{r}{|u_c|}\right)}
\end{array}
\end{equation}

# Potential flow visualization

## Well flow is dominant  

\begin{equation}
\begin{array}{rl}
\dfrac{dw}{dz} =& -KI - \dfrac{Q_{\rm in}}{2\pi H} \left(\dfrac{f}{z} - \dfrac{1}{z-r}\right)\\
\\
\mathcal{F}_L =& \dfrac{K I}{\dfrac{Q_{\rm in}}{2\pi H} \left(\dfrac{f}{z} - \dfrac{1}{z-r}\right)}
\end{array}
\end{equation}


```python
%reset -f
import numpy as np
import matplotlib.pyplot as plt
import matplotlib

def ZFix(x,y):
    XX,YY = np.meshgrid(X,Y)
    return XX + (YY*1j)
    
def regionalFlow():
    return -K*I

def wellFlow(z):
    w = -(QIN/(2*PI*H))*(F/z - 1/(z-R))
    wAbs = np.abs(w)    
    return (w, w.real, -w.imag, wAbs)

def totalFlow(z):
    reg  = regionalFlow()
    well = wellFlow(z)
    tot  = reg + well[0]
    totAbs = np.abs(tot)
    return (tot, tot.real,-tot.imag,totAbs)

def flowNumber(z):
    reg  = regionalFlow() 
    well = wellFlow(z)[3]
    FL = np.abs(reg)/well
    return FL
```


```python
''' GLOBAL CONSTANTS '''
PI = 3.141592
K = 1.0E-2
#I = 6.6E-14
I = 1.0E-10
R = 40.
F = 1.5
QIN = 1.0/86400.
H = 1
#H = 20.

POINT_OUT = {"x":0,
             "y":0,
             "c":"r",
             "s":100,
             "linewidths":1,
             "edgecolors":"w",
             "label":"Extraction"}

POINT_QIN = {"x":R,
             "y":0,
             "c":"b",
             "s":80,
             "linewidths":1,
             "edgecolors":"w",
             "label":"Injection"}

X = np.linspace(-50,150,151)
Y = np.linspace(-100,100,33)
Z = ZFix(X,Y)
_,U,V,C = totalFlow(Z)
FL = flowNumber(Z)

fig, axs = plt.subplots(2,1,figsize=(12,10))
ax = axs[0]
st = ax.streamplot(X,Y,U,V,color="k",density=[0.8, 0.8])
fi = ax.pcolormesh(X,Y,C,shading='auto',vmin=0,vmax=5.0E-7,cmap="cool")
pi = ax.scatter(**POINT_QIN,zorder=3)
po = ax.scatter(**POINT_OUT,zorder=3)

cbar = ax.figure.colorbar(fi,ax=ax,orientation="vertical")
cbar.set_label(r'$\bf{|U|}$')

ax.legend(loc='lower right')

ax = axs[1]
fi = ax.pcolormesh(X,Y,np.log10(FL),shading='auto')
pi = ax.scatter(**POINT_QIN)
po = ax.scatter(**POINT_OUT)

ax.xaxis.set_tick_params(which="both",labelbottom=False,bottom=False)
ax.yaxis.set_tick_params(which="both",labelleft=False,left=False)
cbar = ax.figure.colorbar(fi,ax=ax,orientation="vertical")
ax.legend(loc='lower right')

plt.show()
```


```python
X = np.linspace(0.01,R-0.01,151)
Y = np.array([0.])
Z = ZFix(X,Y)
_,U,V,C = totalFlow(Z)
FL = flowNumber(Z)
fig, axs = plt.subplots(2,1,figsize=(8,6),sharex=True,\
                        gridspec_kw={"height_ratios":[3,1],"hspace":0})
ax = axs[0]
li = ax.plot(X,C[0],lw=5)
ax.set(xlim=[0,R],ylim=[0,1.0E-6],ylabel=r"$\bf{|U|}$ [m/s]")

ax = axs[1]
li = ax.plot(X,np.log10(FL[0]),lw=3,c="gray")
ax.set(ylim=[-7.5,-4.5],xlabel=r"$\bf{X}$ [m]",ylabel=r"$\bf{\log(\mathcal{F}_L)}$")

plt.show()
```

## Regional flow is dominant


```python
%reset -f
import numpy as np
import matplotlib.pyplot as plt
import matplotlib

def ZFix(x,y):
    XX,YY = np.meshgrid(X,Y)
    return XX + (YY*1j)
    
def regionalFlow():
    return -K*I

def wellFlow(z):
    w = -(QIN/(2*PI*H))*(F/z - 1/(z-R))
    wAbs = np.abs(w)    
    return (w, w.real, -w.imag, wAbs)

def totalFlow(z):
    reg  = regionalFlow()
    well = wellFlow(z)
    tot  = reg + well[0]
    totAbs = np.abs(tot)
    return (tot, tot.real,-tot.imag,totAbs)

def flowNumber(z):
    reg  = regionalFlow() 
    well = wellFlow(z)[3]
    FL = np.abs(reg)/well
    return FL
```


```python
''' GLOBAL CONSTANTS '''
PI = 3.141592
K = 1.0E-2
I = 6.6E-4
R = 40.
F = 10
QIN = 1.0/86400.
H = 20.

POINT_OUT = {"x":0,
             "y":0,
             "c":"r",
             "s":100,
             "linewidths":1,
             "edgecolors":"w",
             "label":"Extraction"}

POINT_QIN = {"x":R,
             "y":0,
             "c":"b",
             "s":80,
             "linewidths":1,
             "edgecolors":"w",
             "label":"Injection"}

X = np.linspace(-50,150,151)
Y = np.linspace(-100,100,33)
Z = ZFix(X,Y)
_,U,V,C = totalFlow(Z)
FL = flowNumber(Z)

fig, axs = plt.subplots(2,1,figsize=(12,10))
ax = axs[0]
st = ax.streamplot(X,Y,U,V,color="k",density=[0.8, 0.8])
fi = ax.pcolormesh(X,Y,C,shading='auto',vmin=0,vmax=2.0E-5,cmap="cool")
pi = ax.scatter(**POINT_QIN,zorder=3)
po = ax.scatter(**POINT_OUT,zorder=3)

cbar = ax.figure.colorbar(fi,ax=ax,orientation="vertical")
cbar.set_label(r'$\bf{|U|}$')

ax.legend(loc='lower right')

ax = axs[1]
fi = ax.pcolormesh(X,Y,np.log10(FL),shading='auto')
pi = ax.scatter(**POINT_QIN)
po = ax.scatter(**POINT_OUT)

ax.xaxis.set_tick_params(which="both",labelbottom=False,bottom=False)
ax.yaxis.set_tick_params(which="both",labelleft=False,left=False)
cbar = ax.figure.colorbar(fi,ax=ax,orientation="vertical")
ax.legend(loc='lower right')

plt.show()
```


```python
X = np.linspace(0.01,R-0.01,151)
Y = np.array([0.])
Z = ZFix(X,Y)
_,U,V,C = totalFlow(Z)
FL = flowNumber(Z)
fig, axs = plt.subplots(2,1,figsize=(8,6),sharex=True,\
                        gridspec_kw={"height_ratios":[3,1],"hspace":0})
ax = axs[0]
li = ax.plot(X,C[0],lw=5)
ax.set(xlim=[0,R],ylim=[0,1.0E-5],ylabel=r"$\bf{|U|}$ [m/s]")

ax = axs[1]
li = ax.plot(X,np.log10(FL[0]),lw=3,c="gray")
ax.set(xlabel=r"$\bf{X}$ [m]",ylabel=r"$\bf{\log(\mathcal{F}_L)}$")
ax.yaxis.set_tick_params(which="both",labelright=True,right=True,labelleft=False,left=False)

plt.show()
```

___
# Implementation to find worst cases


```python
%reset -f
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
from os import system
from matplotlib.gridspec import GridSpec 

''' GLOBAL CONSTANTS '''
PI = 3.141592
THETA = 0.35
```

\begin{equation}
\begin{array}{rl}
    \mathcal{F}_L =& \dfrac{4 K I H r}{Q_{in}(1+f)}\\
    u_{c} =& \dfrac{-KI\mathcal{F}_L}{\theta \left(\mathcal{F}_L+1\right)}\\
    \tau =& -\dfrac{r}{|u_{c}|}\\
    C_{\tau,{\rm decay}}=& C_0 \exp{\left(-\lambda \tau \right)}\\
    C_{\tau,{\rm dilut}} =& C_{in} \left( \dfrac{Q_{in}}{u_c \Delta y \Delta z} \right)\\
    C_{\tau,{\rm both}} =& \dfrac{C_{\rm in}Q_{\rm in}}{u_c \Delta y H \theta} \exp{\left(-\lambda\dfrac{r}{|u_c|}\right)}
\end{array}
\end{equation}


```python
def flowNumber():
    return (4.0*K*I*H*r) / (Qin*(1+f))

def uChar():
    return -(K*I*flowNumber())/(THETA*(flowNumber() + 1))

def tChar():
    return -r/uChar()

def cDecay():
    return C0 * np.exp(-decayRate * tChar())

def cDilut():
    return (C0 * Qin) / (-uChar() * delY * delZ * THETA)

def cBoth():
    return (C0 * Qin) / (-uChar() * delY * delZ * THETA) * np.exp(-decayRate * tChar())

def findSweet():
    deltaConc = np.abs(cBoth() - np.max(cBoth()))
    return np.argmin(deltaConc)
```

**Some constant Parameters**


```python
K = 10**-2
Qin = 0.24/86400
f = 10
H = 20
r = 40
I = 0.001
C0 = 1.0
decayRate = 3.5353E-06
delY,delZ = 1.35,H
```


```python
print("Nondim Flow  = {:.2E}".format(flowNumber()))
print("Charac. Vel  = {:.2E} m/s".format(uChar()))
print("Charac. time = {:.2E} s".format(tChar()))
```

    Nondim Flow  = 1.05E+03
    Charac. Vel  = -2.85E-05 m/s
    Charac. time = 1.40E+06 s



```python
print("Rel concenc. due decay = {:.2E}".format(cDecay()))
print("Rel conc. due dilution = {:.2E}".format(cDilut()))
print("Rel conc. due both eff = {:.2E}".format(cBoth()))
```

    Rel concenc. due decay = 7.05E-03
    Rel conc. due dilution = 1.03E-02
    Rel conc. due both eff = 7.26E-05



```python
I = 10**np.linspace(-5,0,num=100)

c1 = cDecay()
c2 = cDilut()
c3 = cBoth()
i = findSweet()

#worstC = np.mean([c1[i],c2[i]])
worstC = c3[i]
worstI = I[i]
```


```python
print(I[i])
```

    0.0047508101621027985



```python
# create objects 
fig, axs = plt.subplots(2,2,sharex=True, sharey=False,\
                        figsize=(12,8),gridspec_kw={"height_ratios":[1,3],"hspace":0.04,"wspace":0.35}) 

bbox = dict(boxstyle='round', facecolor='mintcream', alpha=0.90)

arrowprops = dict(
    arrowstyle="->",
    connectionstyle="angle,angleA=90,angleB=10,rad=5")

annotation = \
    r"$\bf{-\log(C/C_0)} = $" + "{:.1f}".format(-np.log10(worstC)) + \
    "\n@" + r" $\bf{I} = $" + "{:.1E}".format(worstI)

information = \
    r"$\bf{K}$" + " = {:.1E} m/s".format(K) + "\n"\
    r"$\bf{H}$" + " = {:.1f} m".format(H) + "\n"\
    r"$\bf{r}$" + " = {:.1f} m".format(r) + "\n"\
    r"$\bf{Q_{in}}$" + " = {:.2f} m³/d".format(Qin*86400) + "\n"\
    r"$\bf{f}$" + " = {:.1f}".format(f) + "\n"\
    r"$\bf{\lambda}$" + " = {:.2E} 1/s".format(decayRate)

# Ax1 - Relative concentration
ax = axs[1,0]
ax.plot(I,c1,label="Due decay",lw=3,ls="dashed",alpha=0.8)
ax.plot(I,c2,label="Due dilution",lw=3,ls="dashed",alpha=0.8)
#ax.plot(I,np.minimum(c1,c2),label="Overall effect",lw=3,c='k',alpha=0.9)
ax.plot(I,c3,label="Overall effect",lw=3,c='k',alpha=0.9)

ax.set_xscale("log")
ax.set_yscale("log")

ax.set_ylim(1.0E-10,1)
ax.set_xlim(1.0E-4,1.0E-1)

ax.set_xlabel("Water table gradient\n$I$ [m/m]")
ax.set_ylabel("Relative Concentration\n$C/C_0$ [-]")
ax.legend(loc="lower left",shadow=True)

ax.annotate(annotation,(worstI,worstC),
    xytext=(0.05,0.85), textcoords='axes fraction',
    bbox=bbox, arrowprops=arrowprops)

ax.text(0.65,0.05,information,bbox=bbox,transform=ax.transAxes)
ax.axvline(x=I[i], lw = 1, ls = "dashed", c = "red")

####################################
# Ax2 - log-removals
ax = axs[1,1]
ax.plot(I,-np.log10(c1),label="Due decay",lw=3,ls="dashed",alpha=0.8)
ax.plot(I,-np.log10(c2),label="Due dilution",lw=3,ls="dashed",alpha=0.8)
#ax.plot(I,-np.log10(np.minimum(c1,c2)),label="Overall effect",lw=3,c='k',alpha=0.9)
ax.plot(I,-np.log10(c3),label="Overall effect",lw=3,c='k',alpha=0.9)

ax.set_xscale("log")
ax.set_ylim(0,10)
ax.set_xlim(1.0E-4,1.0E-1)
ax.set_xlabel("Water table gradient\n$I$ [m/m]")
ax.set_ylabel("log-reductions\n$-\log(C/C_0)$ [-]")
ax.legend(loc="lower right",shadow=True)

ax.annotate(annotation,(worstI,-np.log10(worstC)),
    xytext=(0.05,0.10), textcoords='axes fraction',
    bbox=bbox, arrowprops=arrowprops)

ax.text(0.65,0.70,information,bbox=bbox,transform=ax.transAxes)
ax.axvline(x=I[i], lw = 1, ls = "dashed", c = "red")

####################################
#Ax3 - Flow number
ax = axs[0,1]
ax.plot(I,flowNumber(),label="flowNumber",lw=3,c="gray")
ax.axhline(y=1.0)
ax.set_xscale("log")
ax.set_yscale("log")
ax.xaxis.set_tick_params(which="both",labeltop='on',top=True,bottom=False)

ax.set_ylabel("Nondim. flow number\n$\mathcal{F}_L$ [-]")
ax.axvline(x=I[i], lw = 1, ls = "dashed", c = "red")

#Ax4 - Flow number
ax = axs[0,0]
ax.plot(I,flowNumber(),label="flowNumber",lw=3,c="gray")
ax.axhline(y=1.0)
ax.set_xscale("log")
ax.set_yscale("log")
ax.xaxis.set_tick_params(which="both",labeltop='on',top=True,bottom=False)

ax.set_ylabel("Nondim. flow number\n$\mathcal{F}_L$ [-]")
ax.axvline(x=I[i], lw = 1, ls = "dashed", c = "red")

plt.show()
```

____
# Find the worst case 
## >> Geometric parameters $H$ and $r$


```python
from drawStuff import *
```


```python
K = 10**-2
Qin = 0.24/86400
f = 10
C0 = 1.0
decayRate = 3.5353E-06

Harray = np.array([2.,5.,10.,20.,50.])
rarray = np.array([5.,10.,40.,100.])
Iarray = 10**np.linspace(-5,0,num=100)

Ci = np.zeros([len(rarray),len(Harray)])
Ii = np.zeros([len(rarray),len(Harray)])
FLi = np.zeros([len(rarray),len(Harray)])

for hi,H in enumerate(Harray):
    for ri,r in enumerate(rarray):
        i = findSweet()

        worstC = -np.log10(cBoth()[i])
        worstGradient = Iarray[i]
        worstFlowNumber = flowNumber()[i]
        
        Ci[ri,hi]  = worstC
        Ii[ri,hi]  = worstGradient
        FLi[ri,hi] = worstFlowNumber
```


```python
myLabels={"Title": { 0: r"$\bf{-\log (C_{\tau}/C_0)}$",
                     1: r"$\bf{I}$ (%)",
                     2: r"$\log(\mathcal{F}_L)$"},
          "Y": "Aquifer thickness\n$\\bf{H}$ (m)",
          "X": "Setback distance\n$\\bf{r}$ (m)"}

threeHeatplots(data={"I":Ii.T,"C":Ci.T,"FL":FLi.T},\
               xlabel=Harray,ylabel=rarray,myLabels=myLabels)
```

## >>Well parameters


```python
K = 10**-2
H = 20
r = 40
C0 = 1.0
decayRate = 3.5353E-06

Qin_array = np.array([0.24,1.0,10.0,100.])/86400.
f_array = np.array([1,10.,100.,1000.,10000.])
Iarray = 10**np.linspace(-5,0,num=100)

Ci = np.zeros([len(Qin_array),len(f_array)])
Ii = np.zeros([len(Qin_array),len(f_array)])
FLi = np.zeros([len(Qin_array),len(f_array)])

for fi,f in enumerate(f_array):
    for qi,Qin in enumerate(Qin_array):
        i = findSweet()
        worstC = -np.log10(cBoth()[i])
        worstGradient = Iarray[i]
        worstFlowNumber = flowNumber()[i]
        
        Ci[qi,fi] = worstC
        Ii[qi,fi] = worstGradient
        FLi[qi,fi] = worstFlowNumber
```


```python
myLabels={"Title": { 0: r"$\bf{-\log (C_{\tau}/C_0)}$",
                     1: r"$\bf{I}$ (%)",
                     2: r"$\log(\mathcal{F}_L)$"},
          "Y": "Extraction to injection ratio\n$\\bf{f}$ (-)",
          "X": "Injection flow rate\n$\\bf{Q_{in}}$ (m³/d)"}

threeHeatplots(data={"I":Ii.T,"C":Ci.T,"FL":FLi.T},\
               xlabel=f_array,ylabel=np.round(Qin_array*86400,decimals=2),myLabels=myLabels)
```

## Hydraulic conductivity


```python
Qin = 0.24/86400
f = 10
C0 = 1.0
decayRate = 3.5353E-06

Karray = 10.**np.array([-1.,-2.,-3.,-4.,-5.])
rarray = np.array([5.,10.,40.,100.])
Iarray = 10**np.linspace(-5,0,num=100)

Ci = np.zeros([len(rarray),len(Karray)])
Ii = np.zeros([len(rarray),len(Karray)])
FLi = np.zeros([len(rarray),len(Karray)])

for ki,K in enumerate(Karray):
    for ri,r in enumerate(rarray):
        i = findSweet()

        worstC = -np.log10(cBoth()[i])
        worstGradient = Iarray[i]
        worstFlowNumber = flowNumber()[i]
        
        Ci[ri,ki] = worstC
        Ii[ri,ki] = worstGradient
        FLi[ri,ki] = worstFlowNumber
```


```python
myLabels={"Title": { 0: r"$\bf{-\log (C_{\tau}/C_0)}$",
                     1: r"$\bf{I}$ (%)",
                     2: r"$\log(\mathcal{F}_L)$"},
          "Y": "Hydraulic conductivity\n$\\bf{K}$ (m/s)",
          "X": "Setback distance\n$\\bf{r}$ (m)"}

threeHeatplots(data={"I":Ii.T,"C":Ci.T,"FL":FLi.T},\
               xlabel=Karray,ylabel=rarray,myLabels=myLabels)
```

>> EXPERIMENTAL


```python

```

## r and Qin


```python
K = 10**-2
H = 20
f = 10
C0 = 1.0
decayRate = 3.5353E-06

Qin_array = np.array([0.24,1.0,10.0,100.])/86400.
rarray = np.array([5.,10.,40.,100.])
Iarray = 10**np.linspace(-5,0,num=100)

Ci = np.zeros([len(rarray),len(Qin_array)])
Ii = np.zeros([len(rarray),len(Qin_array)])
FLi = np.zeros([len(rarray),len(Qin_array)])

for qi,Qin in enumerate(Qin_array):
    for ri,r in enumerate(rarray):
        i = findSweet()

        worstC = -np.log10(cBoth()[i])
        worstGradient = Iarray[i]
        worstFlowNumber = flowNumber()[i]
        
        Ci[ri,qi] = worstC
        Ii[ri,qi] = worstGradient
        FLi[ri,qi] = worstFlowNumber
```


```python
myLabels={"Title": { 0: r"$\bf{-\log (C_{\tau}/C_0)}$",
                     1: r"$\bf{I}$ (%)",
                     2: r"$\log(\mathcal{F}_L)$"},
          "Y": "Injection flow rate\n$\\bf{Q_{in}}$ (m³/d)",
          "X": "Setback distance\n$\\bf{r}$ (m)"}

threeHeatplots(data={"I":Ii.T,"C":Ci.T,"FL":FLi.T},\
               ylabel=rarray,xlabel=np.round(Qin_array*86400,decimals=2),myLabels=myLabels)
```


```python

```

96ebd40a8e71779a84c6b1d7037137930aadfd03
ipynb
Jupyter Notebook
notebooks/Concepts/.ipynb_checkpoints/Potential flow-checkpoint.ipynb
edsaac/bioparticle
67e191329ef191fc539b290069524b42fbaf7e21
MIT
notebooks/Concepts/.ipynb_checkpoints/Potential flow-checkpoint.ipynb
edsaac/bioparticle
67e191329ef191fc539b290069524b42fbaf7e21
MIT
2020-09-25T23:31:21.000Z
2020-09-25T23:31:21.000Z
notebooks/Concepts/.ipynb_checkpoints/Potential flow-checkpoint.ipynb
edsaac/VirusTransport_RxSandbox
67e191329ef191fc539b290069524b42fbaf7e21
MIT
2021-09-30T05:00:58.000Z
2021-09-30T05:00:58.000Z
Qwen/Qwen-72B
1. YES
2. YES