The transient Theis solution

 The partial differential equation that describes hydraulic head in a confined aquifer over time is given as follows \[ \frac{\partial h}{\partial t} = \frac{T}{S}\nabla^2 h \]

 where \(T\) is the transmissivity and \(S\) is the storativity of the aquifer. The laplacian can be written in cylindrical coordinates

\[ \begin{align} \frac{\partial h}{\partial t}=\frac{T}{S}\nabla^2 h =& \frac{T}{S}\left( \frac{\partial^2 h}{\partial r^2}+\frac{1}{r}\frac{\partial h}{\partial r}+\frac{1}{r^2}\frac{\partial^2 h}{\partial \theta^2}+\frac{\partial^2 h}{\partial z^2} \right) \\=& \frac{T}{S} \left( \frac{\partial^2 h}{\partial r^2}+\frac{1}{r}\frac{\partial h}{\partial r} \right)=\frac{T}{S}\frac{\partial}{\partial r} \left( r\frac{\partial h}{\partial r} \right)\end{align}\]

where the flow is assumed to travel only in the radial direction.

It can also be assumed that at the initial time or at an infinite distance from the pumping well the hydraulic head is a constant \(h_0\)
\[
h(r,0) = h_0 \qquad \lim_{r\rightarrow \infty}h(r,t)=h_0
\]

Since the flow rate from the well is \(Q\) and there is only flow in radial direction in the reservoir mass conservation gives
\[ \begin{align}
&Q=KA\frac{\partial h}{\partial r}=Kb\cdot2\pi r\frac{\partial h}{\partial r}\\
\Rightarrow& \lim_{r\rightarrow 0} \left( r \frac{\partial h}{\partial r}\right) =\frac{Q}{2 \pi T}
\end{align} \]

It might be more convenient to descibe those equations in terms of drawdown, rather than hydraulic head, where they can be related by
\[
s=h_0-h \Rightarrow h=h_0-s
\]
Inserting this change of varibles into the partial differential equation that describes the hydraulic head gives
\[
\frac{\partial s}{\partial t} = \frac{T}{S}\left( \frac{\partial^2 s}{\partial r^2}+\frac{1}{r}\frac{\partial s}{\partial r} \right)
\]
with the boundary conditions
\[
s(r,0) = 0 \qquad \lim_{r\rightarrow \infty}s(r,t)=0
\]
and
\[
\lim_{r\rightarrow 0} \left( r\frac{\partial s}{\partial r} \right) =-\frac{Q}{2\pi T}
\]

We then apply another change of variables as suggested by Theis in order to obtain a similarity solution
\[
u=\frac{r^2 S}{4Tt}\quad\text{and}\quad s=\frac{Q}{4\pi T}f(u)
\]
which by applying the chain rule gives
\[ \begin{align}
&\frac{\partial s}{\partial u}\frac{\partial u}{\partial t}= \frac{T}{S} \left( \frac{\partial^2 s}{\partial u^2}\left( \frac{\partial u }{\partial r} \right)^2+\frac{\partial s}{\partial u}\frac{\partial^2 u}{\partial r^2}+\frac{1}{r}\frac{\partial s}{\partial u}\frac{\partial u}{\partial r} \right) \\
\Rightarrow&-\frac{Q}{4\pi Tt}u\frac{\partial f}{\partial u}= \frac{T}{S}\cdot\frac{Q}{4\pi T} \left( \frac{4u^2}{r^2}\frac{\partial^2 f}{\partial u^2}+2\frac{u}{r^2}\frac{\partial f}{\partial u}+2\frac{u}{r^2}\frac{\partial f}{\partial u} \right)\\
\Rightarrow& \frac{\partial f}{\partial u}=-\frac{1}{u}\left( u\frac{\partial^2f}{\partial u^2}+\frac{\partial f}{\partial u} \right)\\
\Rightarrow&\frac{\partial^2 f}{\partial u^2}+\left(1+\frac{1}{u}\right) \frac{\partial f}{\partial u}=0
\end{align} \]
we then multiply through the equation by \(u e^u\) and get the following
\[ \begin{align}
&ue^u\frac{\partial^2 f}{\partial u^2}+(ue^u+e^u)\frac{\partial f}{\partial u}=0\\
\Rightarrow& \frac{\partial}{\partial u} \left( ue^u\frac{\partial f}{\partial u} \right)=0
\end{align} \]

The change of variables is then also applied to the boundary conditions such that
\[
\lim_{r\rightarrow\infty}s(r,t) = \lim_{t\rightarrow 0} s(r,t) = \lim_{u\rightarrow\infty}f(u)=0
\]
and
\[ \begin{align}
&r\frac{\partial s}{\partial r}=\frac{Qr}{4\pi T}\frac{\partial u}{\partial r}\frac{\partial f}{\partial u}=\frac{Qr}{4\pi T}\cdot\frac{2u}{r}\frac{\partial f}{\partial u}=-\frac{Q}{2\pi T}\\
\Rightarrow& \lim_{u\rightarrow 0} \left( u\frac{\partial f}{\partial u} \right) =-1
\end{align} \]

The equation now simplifies to an ordinary differential equation where the boundary conditions above can be applied as follows
\[ \begin{align}
&\frac{d}{d u}\left( ue^u \frac{d f}{d u} \right) = 0\\
\Rightarrow& ue^u\frac{df}{du}=C\\
\Rightarrow& \lim_{u\rightarrow 0}\left( e^u u\frac{df}{du} \right)=e^0\cdot(-1) = -1 \\ \Rightarrow& C=-1
\end{align} \]
which gives
\[
ue^u \frac{df}{du}=-1 \Rightarrow \frac{df}{du}=-\frac{e^{-u}}{u}
\]

The fundamental theorem of calculus states that
\[
\frac{d}{du}\int_a^u g(t) dt = g(u)
\]
where \(a\) is some point on the interval where \(g\) is continious. This gives
\[ \begin{align}
&\frac{df}{du}=\frac{d}{du}\int_{\infty}^u \frac{-e^{-t}}{t}dt=-\frac{e^{-u}}{u}\\
\Rightarrow& f(u) = \int_u^\infty\frac{e^{-t}}{t}dt
\end{align} \]
This function can not be expressed in terms of a finite sum of elementary functions and is known outside of hydrogeology literature as the exponential integral, or
\[
f(u)=W(u)=-Ei(-u)=E_1(u)
\]

The exponential integral, or the well function, can be witten as an infinite Taylor series
\[
E_1 (u) = W(u) = -\gamma-\ln u+\sum_{k=1}^{\infty}\frac{(-1)^{k+1}u^k}{k\cdot k!}
\]
where \(\gamma\approx0.57721\) is the Euler-Mascheroni constant.