1 Introduction

Anomalous or fractional heat and mass diffusion defines deviations from the traditional diffusion model, which describes the movement of molecules as a random walk with a constant diffusion coefficient. In anomalous diffusion, the mean squared displacement of molecules grows with time slower or faster than linearly. With a slower growth, the regime is subdiffusive, meaning molecules spread out more slowly than expected. When the growth is faster than linear, the process is superdiffusive, meaning molecules spread out more rapidly than in standard diffusion. A thorough review on this topic can be found in Henry et al. [1]. In recent work, Barletta [2] has proposed a model for convective motion of a fluid in a saturated porous medium where the diffusion coefficient is time dependent but may be of subdiffusion or superdiffusion type. He showed, for the linearized theory, that regardless of the size of the Rayleigh number, the solution to the perturbation system may grow in the short term, but eventually the perturbation will decay for superdiffusion and will grow indefinitely for subdiffusion. The object of this paper is to analyse the completely nonlinear problem, and we show that exponential decay holds in this case for superdiffusion. We establish this result in a saturated porous material but also allow for variable gravity effects, cf. Straughan [3], or for penetrative convection driven by an internal heat source, cf. Straughan [4,  p. 343]. This extension is important because mathematically the coefficients in the perturbation equations often depend on the vertical coordinate z. For such spatial dependence, it is usually not possible to obtain an analytical solution by a normal mode procedure, cf. Barletta [5, 6], and one must resort to numerical techniques. We additionally establish nonlinear decay of the solution when the porous medium is bidisperse, i.e. the medium has macropores, but the solid skeleton possesses fissures or cracks, necessitating the inclusion of micropores, and hence a bidisperse or double porosity structure. Furthermore, we prove a similar asymptotic decay result when thermal and salt effects are present, i.e. in the case of thermosolutal or double diffusive convection

The extension to spatially dependent coefficients is essential for real life applications. Also, consideration of double diffusive flow in a bidisperse porous medium is a topic of immense importance in current everyday life. For example, this phenomenon is proving important in renewable energy research, especially in solar pond technology, see Dineshkumar and Raja [7], Wang et al. [8]. A particularly important application of double diffusion in a bidisperse porous medium is to magma flow in a volcano, see e.g. Vieira et al. [9], Allocca et al. [10], Toy et al. [11], Singh et al. [12], Bagdassarov and Fradkov [13], DeCampos et al. [14]. Given the recent seismic and volcanic activity in the Campi Flegrei region near Pozzuoli, an understanding of this scenario, and the potential effects of anomalous diffusion, is of vital importance.

We now present the equations for convection models and deliver a fully nonlinear analysis of asymptotic solution behaviour.

2 Thermal convection in a Darcy porous material

The subject of thermal convection in a Darcy porous material is investigated by Barletta [2] who allows for anomalous thermal diffusion by introducing a statistically motivated time dependent diffusion coefficient leading to a diffusion of form \(\varphi D_rrt^{r-1}\Delta C\), where C may be a concentration or a temperature field. The coefficients \(\varphi ,D_r,r\) are positive constants being porosity, constant diffusion coefficient, and an exponent. The variable t denotes time and \(\Delta \) is the three-dimensional Laplace operator

$$\begin{aligned} \Delta = \frac{\partial ^2}{\partial x^2} +\frac{\partial ^2}{\partial y^2} +\frac{\partial ^2}{\partial z^2}. \end{aligned}$$

The basic solution of Barletta [2] is the classical one of Chandrasekhar [15], where C (or T in the case of temperature) is linear in the vertical coordinate z and the velocity is 0.

We commence with the non-dimensional perturbation equations of Barletta [2, eqs. 36–38], which employ a Boussinesq approximation, cf. Barletta [16], and have form (in our notation)

$$\begin{aligned} \begin{aligned}&u_i+\pi _{,i}- Ra\theta k_i=0,\\&u_{i,i}=0,\\&\theta _{,t}+u_i\theta _{,i} = w + k t^{r-1}\Delta \theta , \end{aligned} \end{aligned}$$
(1)

where \(u_i,\pi ,\theta \) are the velocity perturbation, pressure perturbation, and temperature perturbation, \(w=u_3\), \(\textbf{k}=(0,0,1)\), and \(\textrm{Ra}\) is the Rayleigh number. It is convenient to rescale these equations and employ the parameter \(R=\sqrt{Ra}\) and we write (1) in the form

$$\begin{aligned} \begin{aligned}&u_i=-\pi _{,i}+ R\theta k_i,\\&u_{i,i}=0,\\&\theta _{,t}+u_i\theta _{,i}=Rw+kt^{r-1}\Delta \theta . \end{aligned} \end{aligned}$$
(2)

Throughout we investigate the superdiffusion case where \(r>1\).

Barletta [2] essentially linearizes (2) and employs a normal mode analysis, cf. Barletta [5, 6], to derive very interesting novel behaviour. In the interests of encompassing greater physical behaviour, we allow for the effect of variable gravity, cf. Straughan [3], and for penetrative convection driven by an internal heat source, cf. Straughan [4, p. 343], and the perturbation Eq. (2) are replaced by

$$\begin{aligned} \begin{aligned}&u_i=-\pi _{,i}+ Rg(z)\theta k_i,\\&u_{i,i}=0,\\&\theta _{,t}+u_i\theta _{,i}=Rh(z)w+kt^{r-1}\Delta \theta , \end{aligned} \end{aligned}$$
(3)

where g, h are bounded functions of z, k is a positive constant and \(r>1\). The domain for Eq. (3) is \((x,y)\in \mathbb {R}^2\), \(\{z\in (0,1)\}\) and \(t>0\). The boundary conditions are

$$\begin{aligned} w=0,\quad \theta =0,\quad z=0,1, \end{aligned}$$
(4)

together with \(u_i,\pi ,\theta \) being periodic in x, y. The periodicity ensures cellular structure of the convection cells as explained in detail by Chandrasekhar [15, pp. 43–52].

Suppose

$$\begin{aligned} \vert g(z)\vert \le c_2,\qquad \vert g+h\vert \le c_1, \qquad \forall z\in [0,1], \end{aligned}$$
(5)

for constants \(c_1\) and \(c_2\). Then let V be a periodic cell for the perturbation solution to (3). Let \((\cdot , \cdot )\) and \(\Vert \cdot \Vert \) be the inner product and norm on \(L^2(V)\).

Multiply (3)\(_1\) by \(u_i\) and integrate over V, and multiply (3)\(_3\) by \(\theta \) and integrate over V. After integration by parts and use of the boundary conditions, one may show

$$\begin{aligned} 0=-\Vert \textbf{u}\Vert ^2+R(gw,\theta ), \end{aligned}$$
(6)

and

$$\begin{aligned} \frac{\textrm{d}}{\textrm{d}t}\frac{1}{2}\Vert \theta \Vert ^2=R(hw,\theta )-kt^{r-1}\Vert \nabla \theta \Vert ^2. \end{aligned}$$
(7)

We next add (6) and (7) and then use the arithmetic–geometric mean inequality on the result to obtain

$$\begin{aligned} \frac{\textrm{d}}{\textrm{d}t}\frac{1}{2}\Vert \theta \Vert ^2\le \frac{R^2c_1^2}{4}\Vert \theta \Vert ^2-kt^{r-1}\Vert \nabla \theta \Vert ^2. \end{aligned}$$
(8)

The function \(\theta \) satisfies Poincaré’s inequality \(\Vert \nabla \theta \Vert ^2\ge \pi ^2\Vert \theta \Vert ^2\) and then from (8) we derive

$$\begin{aligned} \frac{\textrm{d}}{\textrm{d}t}\Vert \theta \Vert ^2+\Bigl (2kt^{r-1}\pi ^2-\frac{R^2c_1^2}{2}\Bigr )\Vert \theta \Vert ^2\le 0. \end{aligned}$$
(9)

Now employ an integrating factor and integrate (9) to see that

$$\begin{aligned} \Vert \theta (t)\Vert ^2\le \Vert \theta (0)\Vert ^2\,\exp \Bigl (\frac{R^2c_1^2t}{2}-\frac{2k\pi ^2}{r}t^r\Bigr ). \end{aligned}$$
(10)

It follows from (10) that as \(t\rightarrow \infty \), \(\Vert \theta (t)\Vert \) tends to zero very rapidly no matter how large \(\Vert \theta (0)\Vert \) is.

From (6), one shows

$$\begin{aligned} \Vert \textbf{u}\Vert ^2\le R^2c_1^2\Vert \theta \Vert ^2, \end{aligned}$$

and so decay of \(\Vert \theta (t)\Vert \) also implies decay of \(\Vert \textbf{u}(t)\Vert \).

Remark

The exponential decay is obtained regardless of the size of the Rayleigh number \(\textrm{Ra}=R^2\). The result (10) applies also to the Barletta model where \(g=h=1\). Observe that we here employ no linearization and show decay for a solution to the fully nonlinear equations.

3 Convection in a bidisperse porous material

For a single temperature, T, in the macro- and micropores, equations for thermal convection in a bidisperse Darcy porous material are given by Gentile and Straughan [17]. If we adopt the anomalous diffusion term of Barletta [2], then the non-dimensional perturbation equations for thermal convection in a bidisperse porous medium may be shown to be, cf. Gentile and Straughan [17], Straughan [18],

$$\begin{aligned} \begin{aligned}&u_i^f+\xi (u_i^f-u_i^p)=-\pi ^f_{,i}+R\theta k_i,\\&u^f_{i,i}=0,\\&K_ru_i^p+\xi (u_i^p-u_i^f)=-\pi ^p_{,i}+R\theta k_i,\\&u^p_{i,i}=0,\\&\theta _{,t}+(u^f_i+u^p_i)\theta _{,i}=R(w^f+w^p)+kt^{r-1}\Delta \theta . \end{aligned} \end{aligned}$$
(11)

Here, \(u^f_i,u^p_i,\pi ^f,\pi ^p\) and \(\theta \) are nonlinear perturbations to the velocity in the macropores, velocity in the micropores, pressure in the macropores, pressure in the micropores, and temperature, respectively, with \(w^f=u^f_3,w^p=u^p_3.\) The parameters \(\xi \) and \(K_r\) are an interaction coefficient and the relative permeability \(K_r=K^f/K^p\), where \(K^f\) and \(K^p\) are the permeabilities in the macro- and micropores.

Equation (11) hold on \(\{(x,y)\in \mathbb {R}^2\},\) \(\{z\in (0,1)\}, t>0\), and the boundary conditions are

$$\begin{aligned} w^f=0, \quad w^p=0, \quad \theta =0,\qquad \text {on}\quad z=0,1, \end{aligned}$$
(12)

together with \(u_i^f,u_i^p,\pi ^f,\pi ^p,\theta \) being periodic in x, y.

To derive an asymptotic behaviour result, we multiply (11)\(_1\) by \(u_i^f\), (11)\(_3\) by \(u_i^p\) and (11)\(_5\) by \(\theta \) and integrate each over V using integration by parts and the boundary conditions. After addition of the equations for \(u_i^f\) and \(u_i^p\), this leads to

$$\begin{aligned} 0=-\Vert \textbf{u}^f\Vert ^2-K_r\Vert \textbf{u}^p\Vert ^2-\xi \Vert \textbf{u}^f-\textbf{u}^p\Vert ^2 +R(\theta ,w^f+w^p), \end{aligned}$$
(13)

and

$$\begin{aligned} \frac{\textrm{d}}{\textrm{d}t}\frac{1}{2}\Vert \theta \Vert ^2= R(\theta ,w^f+w^p) -kt^{r-1}\Vert \nabla \theta \Vert ^2. \end{aligned}$$
(14)

We now add Eqs. (13) and (14) and we use the arithmetic–geometric mean inequality on the \((\theta ,w^f+w^p)\) terms to arrive at

$$\begin{aligned} \frac{\textrm{d}}{\textrm{d}t}\frac{1}{2}\Vert \theta \Vert ^2\le R^2 \left( 1+\frac{1}{K_r}\right) \Vert \theta \Vert ^2 -kt^{r-1}\Vert \nabla \theta \Vert ^2. \end{aligned}$$
(15)

Now employ Poincaré’s inequality on the last term and use an integrating factor as before to obtain

$$\begin{aligned} \Vert \theta (t)\Vert ^2\le \Vert \theta (0)\Vert ^2\left[ 2R^2\left( 1+\frac{1}{K_r}\right) t-\frac{2k\pi ^2}{r}t^r \right] . \end{aligned}$$
(16)

This establishes that since \(r>1\), \(\Vert \theta (t)\Vert \) decays rapidly as \(t\rightarrow \infty \), for any \(\textrm{Ra}\) and \(\Vert \theta (0)\Vert \).

From (13), one shows

$$\begin{aligned} \Vert \textbf{u}^f\Vert ^2+K_r\Vert \textbf{u}^p\Vert ^2\le R^2\Bigl (1+\frac{1}{K_r}\Bigr )\Vert \theta \Vert ^2. \end{aligned}$$
(17)

Hence, from (16) decay of \(\Vert \theta (t)\Vert \) guarantees decay of both \(\Vert \textbf{u}^f(t)\Vert \) and \(\Vert \textbf{u}^p(t)\Vert \).

4 Double diffusive porous convection

In the case of double diffusive convection in a Darcy porous material when the layer is heated below and salted above or below, the basic solution is as in e.g. Straughan [4, pp. 238–241], and the perturbation equations are given as, Straughan [4, eqs. 14.20–14.22]. If we started at the outset with the anomalous diffusion of Barletta [2] for the temperature of form \(k_1t^{r-1}\Delta \theta \) and for salt as \(k_2t^{s-1}\Delta \phi \), for constants \(k_1,k_2\) and \(r>1, s>1\), where \(\phi \) is now the salt perturbation, then the non-dimensional perturbation equations are

$$\begin{aligned} \begin{aligned}&0=-u_i+R\theta k_i-R_s\phi k_i,\\&u_{i,i}=0,\\&\theta _{,t}+u_i\theta _{,i}=Rw+k_1t^{r-1}\Delta \theta ,\\&\textrm{Le}(\phi _{,t}+u_i\phi _{,i})=\pm R_sw+k_2t^{s-1}\Delta \phi , \end{aligned} \end{aligned}$$
(18)

where \(\textrm{Ra}=R^2\) is the Rayleigh number, \(\mathcal {C}=R_s^2\) is the salt Rayleigh number and \(\textrm{Le}\) is a constant known as the Lewis number.

The boundary conditions are

$$\begin{aligned} w=0,\quad \theta =0,\quad \phi =0,\qquad z=0,1, \end{aligned}$$
(19)

together with periodicity in x, y.

We allow for different exponents of anomalous diffusion and without loss of generality we here assume \(s>r>1\). We consider the case of the minus sign in (18)\(_4\) which corresponds to salting from above. The analysis in the other case of the plus sign is easier and we omit details.

To derive an asymptotic result, we note that for \(t\in (0,1), t^{r-1}>t^{s-1}\). From (18)\(_1\), we may obtain

$$\begin{aligned} 0=-\Vert \textbf{u}\Vert ^2+R(\theta ,w)-R_s(\phi ,w). \end{aligned}$$
(20)

From (18)\(_3\) and (18)\(_4\), one finds

$$\begin{aligned} \begin{aligned} \frac{\textrm{d}}{\textrm{d}t}\Bigl ( \frac{1}{2}\Vert \theta \Vert ^2 +\frac{\textrm{Le}}{2}\Vert \phi \Vert ^2 \Bigr ) =&R(w,\theta )-k_1t^{r-1}\Vert \nabla \theta \Vert ^2\\&-R_s(w,\phi )-k_2t^{s-1}\Vert \nabla \phi \Vert ^2. \end{aligned} \end{aligned}$$
(21)

We add (20) and (21) and use the arithmetic–geometric mean inequality and Poincaré’s inequality to obtain

$$\begin{aligned} \frac{\textrm{d}}{\textrm{d}t}\Bigl ( \frac{1}{2}\Vert \theta \Vert ^2 +\frac{\textrm{Le}}{2}\Vert \phi \Vert ^2 \Bigr ) \le&(2R^2-k_1t^{r-1}\pi ^2)\Vert \theta \Vert ^2 +(2R^2_s-k_2t^{s-1}\pi ^2)\Vert \phi \Vert ^2\nonumber \\ \le&(2R^2-k_1t^{s-1}\pi ^2)\Vert \theta \Vert ^2 +(2R^2_s-k_2t^{s-1}\pi ^2)\Vert \phi \Vert ^2. \end{aligned}$$
(22)

Now put \(A=\max \{4R^2,4R_s^2/\textrm{Le}\},\) \(B=\min \{2k_1,2k_2/\textrm{Le}\}\), and from (22) we may obtain

$$\begin{aligned} \frac{\textrm{d}F}{\textrm{d}t}+(B\pi ^2t^{s-1}-A)F\le 0, \end{aligned}$$
(23)

for \(t\in (0,1]\), where \(F(t)=\Vert \theta \Vert ^2+\textrm{Le}\Vert \phi \Vert ^2\). Integrate (23) with an integrating factor to obtain

$$\begin{aligned} F(1)\exp \Bigl (\frac{B\pi ^2}{s}-A\Bigr )\le F(0). \end{aligned}$$

We now apply a similar argument to the above on the time interval \((1,\infty )\), observing then that \(t^{s-1}>t^{r-1}>1.\) In this case, we obtain instead of (23),

$$\begin{aligned} \frac{\textrm{d}F}{\textrm{d}t}+(B\pi ^2t^{r-1}-A)F\le 0, \end{aligned}$$
(24)

This is integrated with an integrating factor to arrive at

$$\begin{aligned} F(t)\exp \Bigl (\frac{B\pi ^2}{r}t^r-At\Bigr )\le&F(1)\exp \Bigl (\frac{B\pi ^2}{r}-A\Bigr )\\ \le&F(0)\exp \Bigl [\frac{B\pi ^2(s-r)}{rs}\Bigr ]=G(0), \end{aligned}$$

where G(0) is as defined. Thus,

$$\begin{aligned} \Vert \theta (t)\Vert ^2+\textrm{Le}\Vert \phi (t)\Vert ^2\le G(0)\exp \Bigl (At-\frac{B\pi ^2}{r}t^r\Bigr ). \end{aligned}$$
(25)

Inequality (25) shows that as \(t\rightarrow \infty \), \(\Vert \theta (t)\Vert \) and \(\Vert \phi (t)\Vert \) both decay rapidly regardless of the size of \(\textrm{Ra}\), \(\mathcal {C}\), or the initial data.

By using the arithmetic–geometric mean inequality in (20), one shows

$$\begin{aligned} \Vert \textbf{u}\Vert ^2\le 2R^2\Vert \theta \Vert ^2+2R_s^2\Vert \phi \Vert ^2. \end{aligned}$$
(26)

Then from (25) and (26), \(\Vert \textbf{u}\Vert \) likewise decays as \(t\rightarrow \infty \).

5 Conclusions

We have extended the interesting result of Barletta [2] for the asymptotic behaviour of the solution to convection in a Darcy porous material for a superdiffusion model to the fully nonlinear case and we have shown that the perturbation velocity and temperature will always decay to zero, at least in \(L^2\) norm. We have shown that this result may be extended to other convection in porous media scenarios. In particular, we allow for effects such as variable gravity or convection with internal heat source. We also established an asymptotic decay result in the important problems of bidisperse convection and double diffusive convection.