# Modeling Transient Flows in Heterogeneous Layered Porous Media Using the Space-Time Trefftz Method
## 2 The Methodology
### The Space-Time Trefftz Method  
The space-time Trefftz method utilizing Trefftz functions originates from the Trefftz method [12; 21; 24]. For a simply connected domain, one usually locates the source point in the space-time domain; the number of source points is only one for the proposed space-time Trefftz method. The approximation of Equation (1) can be expressed as follows:  
\[h(r,\ \theta,\ t)\approx\sum\limits_{p=1}^{m}\sum\limits_{q=1}^{s}A_{pq}\mathbf{ T}_{pq}(r,\ \theta,\ t), \tag{5}\]  
where \(A_{pq}\) is a unknown coefficient, \(\mathbf{T}_{pq}(r,\ \theta,\ t)\) is the basis functions [24], and \(m\) and \(s\) denote the order of the basis functions. Considering the positive basis functions, we have the following:  
\[\mathbf{T}_{pq}(r,\ \theta,\ t)=\Big{\{}r^{q}\cos(q\theta),r^{q}\sin(q\theta),e^{ -p^{2}Dt}J_{0}\Big{(}pr)e^{-p^{2}Dt}J_{q}(pr)\cos(q\theta),e^{-p^{2}Dt}J_{q}( pr)\sin(q\theta)\Big{\}}, \tag{6}\]  
where \(D\) is the hydraulic diffusivity expressed as \(D=k/S_{S}\), \(J_{0}\) is the Bessel function of the first kind, and \(J_{q}\) is the Bessel function of the first kind of the \(q\) order. Considering transient flows in heterogeneous layered porous media, we may write the following:  
\[p=g\sqrt{\frac{D_{1}}{D_{i}}}, \tag{7}\]  
where \(g\) ranges from \(1\) to \(m\); \(D_{1}\) is the hydraulic diffusivity of first subdomain; and \(i\) denotes the \(i\)th subdomain, which means that the proposed method is applicable to transient flows in porous media including many layers. Inserting Equation (7) into Equation (6) obtains the following:  
\[\mathbf{T}_{pqi}(r,\ \theta,\ t)=\Bigg{\{}r^{q}\cos(q\theta),r^{q}\sin(q \theta),e^{-g^{2}D_{1}t}J_{0}(g\sqrt{\frac{D_{1}}{D_{i}}}r)e^{-g^{2}D_{1}t}I_{ q}(g\sqrt{\frac{D_{1}}{D_{i}}}r)\cos(q\theta),e^{-g^{2}D_{1}t}I_{q}(g\sqrt{ \frac{D_{1}}{D_{i}}}r)\sin(q\theta)\Bigg{\}}. \tag{8}\]  
Applying the Dirichlet boundary condition obtains the following:  
\[\begin{array}{l}h_{i}(r,\ \theta,\ t)\approx A_{1}+\sum\limits_{q=1}^{s}A_{ 2q}r^{q}\cos(q\theta)+\sum\limits_{q=1}^{s}A_{3q}r^{q}\sin(q\theta)+\sum \limits_{g=1}^{m}A_{4g}e^{-g^{2}\hat{I}}J_{0}(g\sqrt{\frac{D_{1}}{D_{i}}}r)+ \\ \sum\limits_{g=1}^{m}\sum\limits_{q=1}^{s}I_{q}(g\sqrt{\frac{D_{1}}{D_{i}}}\hat {r})[A_{5 gq}\cos(q\theta)+A_{6 gq}\sin(q\theta)]e^{-g^{2}\hat{I}},\end{array} \tag{9}\]  
where \(A_{1},A_{2q},\cdots,A_{6 gq}\) denote arbitrary constants that have to be evaluated, \(\hat{r}\) is the dimensionless radius described as \(\hat{r}=r/R_{0}\), \(R_{0}\) is the characteristic length, and \(\hat{t}\) is the dimensionless time written as \(\hat{t}=tD_{1}/R_{0}^{2}\). Applying the Neumann boundary condition obtains the following:  
\[\frac{\partial h_{i}(r,\ \theta,\ t)}{\partial n}\approx\sum\limits_{q=1}^{s}A_{ 2q}T_{2}{}^{\prime}+\sum\limits_{q=1}^{s}A_{3q}T_{3}{}^{\prime}+\sum\limits_{g =1}^{m}A_{4g}T_{4}{}^{\prime}+\sum\limits_{g=1}^{m}\sum\limits_{q=1}^{s}A_{5 gq }T_{5}{}^{\prime}+\sum\limits_{g=1}^{m}\sum\limits_{q=1}^{s}A_{6 gq}T_{6}{}^{ \prime}. \tag{10}\]  
where \(n\) is the outward normal direction. The notation in Equation (10) is listed in Table 1.