Skip to main content
Log in

Numerical modeling of a memory-based diffusivity equation and determination of its fractional order value

  • Original Paper
  • Published:
Computational Geosciences Aims and scope Submit manuscript

Abstract

Conventional diffusion equations for fluid flow through porous media do not consider the effects of the history of rock, fluid, and flow. This limitation can be overcome by the incorporation of “memory” in the model, using fractional-order derivatives. Inclusion of fractional-order derivatives in the diffusion equation, however, adds complexity to both the equation and its numerical approximation. Of particular importance is the choice of temporal mesh, which can dramatically affect the convergence of the scheme. In this article, we consider a memory-based radial diffusivity equation, discretized on either uniform or graded meshes. Numerical solutions obtained from these models are compared against analytical solutions, and it is found that the simulation using properly chosen graded meshes gives substantially smaller errors compared to that using uniform meshes. Experimental data from one-dimensional flow through a porous layer with constant pressure gradient are collected from the literature and used to fit the fractional order in the diffusivity equation considered here. A reasonable value of the fractional order is found to be 0.05; this is further validated by performing numerical simulations to match these experiments, demonstrating substantial improvement over the classical Darcy’s model.

This is a preview of subscription content, log in via an institution to check access.

Access this article

Price excludes VAT (USA)
Tax calculation will be finalised during checkout.

Instant access to the full article PDF.

Institutional subscriptions

Similar content being viewed by others

References

  1. Caputo, M: Diffusion of fluids in porous media with memory. Geothermics 28(1), 113–130 (1999)

    Article  Google Scholar 

  2. Caputo, M: Models of flux in porous media with memory. Water Resour. Res. 36(3), 693–705 (2000)

    Article  Google Scholar 

  3. Christensen, R: Theory of viscoelasticity: an introduction. Elsevier, Amsterdam (2012)

  4. Hossain, M E, Mousavizadegan, S H, Ketata, C, Islam, M R: A novel memory based stress strain model for reservoir characterization. Nat Sci Sustainable Technol Res Prog 1, 1–29 (2008)

    Google Scholar 

  5. Zavala-Sanchez, V, Dentz, M, Sanchez-Vila, X: Characterization of mixing and spreading in a bounded stratified medium. Adv. Water Resour. 32(5), 635–648 (2009)

    Article  Google Scholar 

  6. Zhang, HM: Driver memory, traffic viscosity and a viscous vehicular traffic flow model. Transp. Res. B Methodol. 37(1), 27–41 (2003)

    Article  Google Scholar 

  7. Caputo, M, Plastino, W: Diffusion in porous layers with memory. Geophys. J. Int. 158(1), 385–396 (2004)

    Article  Google Scholar 

  8. Iaffaldano, G, Caputo, M, Martino, S: Experimental and theoretical memory diffusion of water in sand. Hydrol. Earth Syst. Sci. Discuss. 2(4), 1329–1357 (2005)

    Google Scholar 

  9. Di Giuseppe, E, Moroni, M, Caputo, M: Flux in porous media with memory: models and experiments. Transp. Porous Media 83(3), 479–500 (2010)

    Article  Google Scholar 

  10. Hossain, M E, Mousavizadegan, S H, Islam, M R, et al.: A new porous media diffusivity equation with the inclusion of rock and fluid memories. E-Library, SPE– 114287 (2008)

  11. Abu-Saman, A M, Assaf, A M: Stability and convergence of Crank-Nicholson method for fractional advection dispersion equation. Advances in Applied Mathematical Analysis 2(2), 117–125 (2007)

    Google Scholar 

  12. Chen, C-M, Liu, F, Turner, I, Anh, V: A Fourier method for the fractional diffusion equation describing sub-diffusion. J. Comput. Phys. 227(2), 886–897 (2007)

    Article  Google Scholar 

  13. Chen, S, Liu, F, Zhuang, P, Anh, V: Finite difference approximations for the fractional Fokker–Planck equation. Appl. Math. Model. 33(1), 256–273 (2009)

    Article  Google Scholar 

  14. Chen, C-M, Liu, F, Anh, V, Turner, I: Numerical schemes with high spatial accuracy for a variable-order anomalous subdiffusion equation. SIAM J. Sci. Comput. 32(4), 1740–1760 (2010)

    Article  Google Scholar 

  15. Cui, M: Compact finite difference method for the fractional diffusion equation. J. Comput. Phys. 228(20), 7792–7804 (2009)

    Article  Google Scholar 

  16. Du, R, Cao, WR, Sun, ZZ: A compact difference scheme for the fractional diffusion-wave equation. Appl. Math. Model. 34(10), 2998–3007 (2010)

    Article  Google Scholar 

  17. Gao, G-H, Sun, Z-Z: A compact finite difference scheme for the fractional sub-diffusion equations. J. Comput. Phys. 230(3), 586–595 (2011)

    Article  Google Scholar 

  18. Langlands, TAM, Henry, B I: The accuracy and stability of an implicit solution method for the fractional diffusion equation. J. Comput. Phys. 205(2), 719–736 (2005)

    Article  Google Scholar 

  19. Liu, F, Zhuang, P, Anh, V, Turner, I: A fractional-order implicit difference approximation for the space-time fractional diffusion equation. ANZIAM Journal 47, 48–68 (2006)

    Article  Google Scholar 

  20. Liu, Q, Gu, Y, Zhuang, P, Liu, F, Nie, YF: An implicit RBF meshless approach for time fractional diffusion equations. Comput. Mech. 48(1), 1–12 (2011)

    Article  Google Scholar 

  21. Lynch, V E, Carreras, B A, del Castillo-Negrete, D, Ferreira-Mejias, KM, Hicks, HR: Numerical methods for the solution of partial differential equations of fractional order. J. Comput. Phys. 192(2), 406–421 (2003)

    Article  Google Scholar 

  22. Meerschaert, M M, Tadjeran, C: Finite difference approximations for fractional advection–dispersion flow equations. J. Comput. Appl. Math. 172(1), 65–77 (2004)

    Article  Google Scholar 

  23. Murillo, J Q, Yuste, S B: On three explicit difference schemes for fractional diffusion and diffusion-wave equations. Phys. Scr. 2009(T136), 014025 (2009)

    Article  Google Scholar 

  24. Murillo, J Q, Yuste, S B: An explicit difference method for solving fractional diffusion and diffusion-wave equations in the Caputo form. J. Comput. Nonlinear Dyn. 6(2), 021014 (2011)

    Article  Google Scholar 

  25. Sun, Z-Z, Wu, X: A fully discrete difference scheme for a diffusion-wave system. Appl. Numer. Math. 56(2), 193–209 (2006)

    Article  Google Scholar 

  26. Tadjeran, C, Meerschaert, M M, Scheffler, H-P: A second-order accurate numerical approximation for the fractional diffusion equation. J. Comput. Phys. 213(1), 205–213 (2006)

    Article  Google Scholar 

  27. Wang, K, Wang, H: A fast characteristic finite difference method for fractional advection–diffusion equations. Adv. Water Resour. 34(7), 810–816 (2011)

    Article  Google Scholar 

  28. Yuste, S B, Acedo, L: An explicit finite difference method and a new Von Neumann-type stability analysis for fractional diffusion equations. SIAM J. Numer. Anal. 42(5), 1862–1874 (2005)

    Article  Google Scholar 

  29. Zhuang, P, Liu, F: Implicit difference approximation for the time fractional diffusion equation. J. Appl. Math. Comput. 22(3), 87–99 (2006)

    Article  Google Scholar 

  30. Zhuang, P, Liu, F, Anh, V, Turner, I: New solution and analytical techniques of the implicit numerical method for the anomalous subdiffusion equation. SIAM J. Numer. Anal. 46(2), 1079–1095 (2008)

    Article  Google Scholar 

  31. Deng, W: Finite element method for the space and time fractional Fokker–Planck equation. SIAM J. Numer. Anal. 47(1), 204–226 (2008)

    Article  Google Scholar 

  32. Roop, J P: Computational aspects of FEM approximation of fractional advection dispersion equations on bounded domains in R2. J. Comput. Appl. Math. 193(1), 243–268 (2006)

    Article  Google Scholar 

  33. Chen, C-M, Liu, F, Burrage, K: Finite difference methods and a Fourier analysis for the fractional reaction–subdiffusion equation. Appl. Math. Comput. 198(2), 754–769 (2008)

    Google Scholar 

  34. Stynes, M, O’Riordan, E, Gracia, J L: Error analysis of a finite difference method on graded meshes for a time-fractional diffusion equation. SIAM J. Numer. Anal. 55(2), 1057–1079 (2017)

    Article  Google Scholar 

  35. Oldham, K, Spanier, J: The fractional calculus theory and applications of differentiation and integration to arbitrary order. Elsevier, Amsterdam, vol. 111 (1974)

  36. Sochi, T: Non-Newtonian flow in porous media. Polymer 51(22), 5007–5023 (2010)

    Article  Google Scholar 

  37. Caputo, M: Diffusion with space memory modelled with distributed order space fractional differential equations. Annals of Geophysics, 46(2) (2003)

  38. Gaspar, F J, Rodrigo, C: Multigrid waveform relaxation for the time-fractional heat equation. SIAM J. Sci. Comput. 39(4), A1201–A1224 (2017)

    Article  Google Scholar 

  39. Elias, B P, Hajash, A: Changes in quartz solubility and porosity due to effective stress: An experimental investigation of pressure solution. Geology 20(5), 451–454 (1992)

    Article  Google Scholar 

  40. Bear, J: Dynamics of fluids in porous media. Courier Corporation, North Chelmsford (2013)

  41. Velarde, J, Blasingame, TA, McCain, W.D. Jr, et al: Correlation of black oil properties at pressures below bubble point pressure-a new approach. Annual Technical Meeting, Petroleum Society of Canada (1997)

Download references

Acknowledgments

The authors would like to thank the Natural Sciences and Engineering Research Council of Canada (NSERC); Research and Development Corporation of Newfoundland and Labrador (RDC), funding no. 210992; and Statoil Canada Ltd., funding no. 211162 for providing financial support to accomplish this research. The work of Scott MacLachlan was partially supported by an NSERC Discovery Grant.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Tareq Uz Zaman.

Additional information

Publisher’s note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Appendices

Appendix A

To find an analytical solution, we consider C1 = C2 = 1 in Eq. 3 giving

$$ \frac{\partial}{\partial x} \left[\frac{\partial^{\alpha}}{\partial t^{\alpha}} \left( \frac{\partial p}{\partial x}\right)\right] = \frac{\partial p}{\partial t}, $$
(A1)

with boundary conditions p(0) = p(1) = 0. We write the solution in series form as

$$ p(x,t) = \sum\limits_{k=1}^{\infty} T_{k}(t) \sin(k\pi x), $$

noting that

$$ \begin{array}{@{}rcl@{}} \frac{\partial p}{\partial t} (x,t) &=& \sum\limits_{k=1}^{\infty} T_{k}^{\prime} (t) \sin(k\pi x) = \\ &-&\sum\limits_{k=1}^{\infty} k^{2} \pi^{2} \frac{\partial^{\alpha} T_{k}(t)}{\partial t^{\alpha}} \sin(k\pi x) = \frac{\partial}{\partial x} \left[\frac{\partial^{\alpha}}{\partial t^{\alpha}} \left( \frac{\partial p}{\partial x}\right)\right]. \end{array} $$

To be a solution, we require that

$$ T_{k}^{\prime} (t) = -k^{2} \pi^{2} \frac{\partial^{\alpha} T_{k}(t)}{\partial t^{\alpha}}, $$

and that

$$ T_{k} (0) = {\upbeta}_{k}, $$

where βk comes from the sine series expansion of the initial data, \(p(x,0) = {\sum }_{k=1}^{\infty } {\upbeta }_{k} \sin \limits (k\pi x)\). Taking Laplace transforms, we have

$$ \mathscr{L} \left[ T_{k}^{\prime} (t) \right] = s \widehat{T_{k}} (s) - \beta_{k}. $$

For the Riemann-Liouville definition of the fractional-order derivative,

$$ \mathscr{L} \left[\frac{\partial^{\alpha} T_{k}(t)}{\partial t^{\alpha}}\right] = s^{\alpha} \widehat{T_{k}} (s) - \left[D^{\alpha -1}T_{k}(t)\right]_{t=0}. $$
(A2)

Therefore,

$$ \begin{array}{@{}rcl@{}} s \widehat{T_{k}}(s)-\beta_{k} = -k^{2} \pi^{2}\left( s^{\alpha} \widehat{T_{k}}(s)-c_{k}\right) \qquad \\ \text{where} \quad c_{k} = \left[D^{\alpha - 1}T_{k} (t)\right]_{t=0} \end{array} $$
$$ \text{or,} \quad \widehat{T_{k}}(s) = \frac{\beta_{k} + k^{2} \pi^{2} c_{k}}{s+k^{2} \pi^{2} s^{\alpha}}=\frac{(\beta_{k}+k^{2} \pi^{2} c_{k})s^{-\alpha}}{s^{1-\alpha}+k^{2}\pi^{2}}. $$

This gives

$$ T_{k}(t) = \left( \beta_{k}+k^{2} \pi^{2} c_{k}\right)E_{1-\alpha}\left( -k^{2} \pi^{2} t^{1-\alpha}\right), $$

where E1−α(v) is the Mittag-Leffler function, and is defined for (1 − α)> 0 as

$$ E_{1-\alpha}(v) = \sum\limits_{k=0}^{\infty} \frac{v^{k}}{{\Gamma} \left( \left( 1-\alpha\right)k+1\right)}, $$
$$ \text{since} \quad \mathscr{L} \left[E_{1-\alpha}\left( -k^{2} \pi^{2} t^{1-\alpha}\right)\right] = \frac{s^{-\alpha}}{s^{1-\alpha}+k^{2} \pi^{2}}. $$

Now,

$$ \begin{array}{@{}rcl@{}} D^{\alpha - 1} \left[({\upbeta}_{k}+k^{2}\pi^{2}c_{k})E_{1-\alpha}(-k^{2}\pi^{2}t^{1-\alpha})\right] = \\ \frac{\beta_{k}+k^{2}\pi^{2}c_{k}}{-k^{2} \pi^{2}} \left( E_{1-\alpha}\left( -k^{2}\pi^{2}t^{1-\alpha}\right)-1\right), \end{array} $$

which forces \(D^{\alpha - 1} \left [T_{k}(t)\right ]_{t=0}=0\) as E1−α(0) = 1, giving ck = 0, \(T_{k}(t) = {\upbeta }_{k} E_{1-\alpha }\left (-k^{2} \pi ^{2} t^{1-\alpha }\right )\), and

$$ p(x,t) = \sum\limits_{k=1}^{\infty} {\upbeta}_{k} E_{1-\alpha} \left( -k^{2} \pi^{2} t^{1-\alpha}\right) \sin (k \pi x). $$
(A3)

Equation A3 is the general analytical solution of Eq. A1 for any initial condition. Now, taking the initial condition considered above,

$$ p_{0}(x) = x(1-x), $$

we have (from standard orthogonality)

$$ \beta_{k} = \frac{{{\int}_{0}^{1}} p_{0}(x) \sin(k \pi x) dx}{{{\int}_{0}^{1}} \sin^{2} (k \pi x) dx}. $$

Calculation gives

$$ \beta_{k} = \frac{4}{k^{3} \pi^{3}}\left[1-(-1)^{k}\right]. $$
(A4)

Note also that were we to replace the Riemann-Liouville definition of the fractional order derivative by the Caputo definition, Eq. A2 would be replaced by

$$ \mathscr{L} \left[\frac{\partial^{\alpha} T_{k}(t)}{\partial t^{\alpha}}\right] = s^{\alpha} \widehat{T_{k}} (s) - s^{\alpha-1}{\upbeta}_{k}, $$

leading to

$$ p(x,t) = p(x,0), $$

which is a constant-in-time solution. Here, we see that this constant-in-time solution must arise when the Caputo derivative is used in the fractional diffusion equation containing a first derivative in time and the fractional derivative mixed within the spatial derivatives.

Cross-sections of the analytical solution given in Eq. A3 for βk as in Eq. A4 are plotted in Fig. 9 for different values of α. These plots show that the inclusion of a temporal derivative of order α result in solutions that dampen more quickly towards constant solutions, and that higher values of α lead to faster initial decay. This is emphasized in Fig. 10, where we plot the temporal behavior of the dominant term in Eq. A3, showing the fast initial decay for α > 0. This behavior of the fractional temporal derivative makes the solution difficult to approximate on uniform meshes. Hence, it is helpful to grade the temporal mesh so that the time steps are smaller in the beginning.

Fig. 9
figure 9

Analytical solution (Eq. A3) of Eq. A1 for different values of α

Fig. 10
figure 10

Change in the values of E1−α(−π2t1−α) with t for different α values. For α = 0, \(E_{1-\alpha }(-\pi ^{2}t^{1-\alpha }) = e^{-\pi ^{2}t}\)

Appendix B

We compute density values at different pressures below bubblepoint pressure using the correlations developed in [41]. From the given density and pressure values, we find the linear equation of best fit as

$$ \rho = 0.7298-0.0004p, $$
(B1)

where ρ is in g/cc, and p is in atm. This relation is derived by considering “typical” values for the physical variables, of reservoir temperature = 185 0F, stock tank oil gravity = 43.7 0API, solution gas-oil ratio at bubble point pressure = 941 scf/stb, gas specific gravity measured at separator = 0.735, and gas specific gravity = 0.865. Rescaling x, t, p, ρ as \(\hat {x}=\frac {x}{x_{max}}\), \(\hat {t}=\frac {t}{t_{max}}\), \(\hat {p}=\frac {p}{p_{i}}\), \(\hat {\rho }=\frac {\rho }{0.7298}\), respectively and taking pi = 250 atm(3650 psi), Eq. B1 can be written as

$$ \hat{\rho} = 1-0.137 \hat{p}. $$
(B2)

We use Eq. B2 on a nondimensional unit domain in Section 2.7.

Rights and permissions

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Zaman, T.U., MacLachlan, S. & Hossain, M.E. Numerical modeling of a memory-based diffusivity equation and determination of its fractional order value. Comput Geosci 25, 655–669 (2021). https://doi.org/10.1007/s10596-020-09986-x

Download citation

  • Received:

  • Accepted:

  • Published:

  • Issue Date:

  • DOI: https://doi.org/10.1007/s10596-020-09986-x

Keywords

Navigation