## Abstract

A bidispersive porous material is one which has usual pores but additionally contains a system of micro pores due to cracks or fissures in the solid skeleton. We present general equations for thermal convection in a bidispersive porous medium when the permeabilities, interaction coefficient and thermal conductivity are anisotropic but symmetric tensors. In this case, we show exchange of stabilities holds and fluid movement will commence via stationary convection, and additionally we show the global nonlinear stability threshold is the same as the linear instability one. Attention is then focused on the case where the interaction coefficient and thermal conductivity are isotropic, and the permeability is isotropic in the horizontal directions, although the permeability in the vertical direction is different. The nonlinear stability threshold is calculated in this case and numerical results are presented and discussed in detail.

## 1. Introduction

Double porosity, or bidispersive, materials are occupying much research attention. A double porosity material is one where there are the usual pores, known as macro pores, but there are also cracks or fissures in the solid skeleton and these give rise to a micro porosity. Theories for non-isothermal fluid flow in a bidispersive porous material were developed by Nield & Kuznetsov [1–3], see also Nield [4], and these theories which allow for independent velocity, pressure and temperature fields in the macro and micro phases are described in detail in the books by Straughan [5,6]. A simpler theory which retains independent velocity and pressure fields but restricts attention to a single temperature field, *T*(**x**,*t*), where **x** is a spatial variable and *t* denotes time, is developed from the Nield & Kuznetsov [1] model in Falsaperla *et al.* [7] and further studies in this theory are analysed in Gentile & Straughan [8,9].

Fluid flow, and in particular, thermal convection in a single porosity fluid saturated porous medium has been the subject of much recent research. Particular interest has been in cases where the permeability is anisotropic (cf. [10–19]). A case of special practical interest is where the permeability is isotropic in the horizontal, *x*,*y*, directions but has a different value in the vertical, *z*, direction. This is known as horizontally isotropic permeability and thermal convection in this theory is analysed in, e.g. Capone *et al.* [10–12], and in Hill & Morad [15]. In general, when the permeability in one direction is different from that in the orthogonal directions, the material is known as transversely isotropic. We adopt the notation horizontally isotropic permeability to distinguish with the general case of transverse isotropy where the axis is not necessarily in the vertical direction (cf. [19,20]).

The major reason for the interest in horizontally isotropic permeability is that many real life porous materials possess this property, such as soils or rocks, and as Karmakar & Raja Sekhar [16] point out, analysis of such porous media are important in hydrocarbon recovery.

Let us denote by *K*_{V} the value of the vertical permeability and by *K*_{H} the value of the horizontal permeability in a horizontally isotropic porous medium. We define the ratio ℓ^{2} by ℓ^{2}=*K*_{H}/*K*_{V}. Fazelalani [21] analyses various real rocks and while *K*_{V}/*K*_{H} is often less than one there are cases where *K*_{V}/*K*_{H}=2.4484 and *K*_{V}/*K*_{H}=82.624. Ayan *et al.* [22] study many routine core samples and find *K*_{V}/*K*_{H} often has small values, for example, for limestone they report values of 0.0362, 0.03679 and 0.270. They also write that, …‘the importance of permeability anisotropy to sound reservoir management is not in dispute …so that reservoir models may be refined, leading to better field development strategies, such as enhanced oil recovery programs and infill well placement’. Widarsono *et al.* [23] analyse rock samples from Indonesia and they report a range of *K*_{V}/*K*_{H} values for sandstone and for carbonate rocks, the values they find being typically in the range 0<*K*_{V}/*K*_{H}<1.2 for sandstone, whereas 0<*K*_{V}/*K*_{H}<4.2 for carbonate rocks.

Application areas for double porosity media, and especially those with anisotropic permeability, are many and include carbon sequestration, [15,24]; landslides, [25–27]; hydraulic fracturing for natural gas (fracking), [28]; drinking water recovery from an aquifer, [29,30]; oil recovery from an underground reservoir, [16,31]; and nuclear waste treatment, [32]. These and many other applications are discussed in the books by Straughan [5,6].

The purpose of the present work is to present what we believe is the first analysis of thermal convection in a bidispersive porous medium allowing for an anisotropic permeability. We begin with a general case where the permeability tensors for the macro and micro phases are symmetric along with the interaction coefficients between the macro and micro velocities. We show in the general case that the onset of convection is by stationary convection and oscillatory convection will not occur. We then show that an analysis of the theory of linear instability will also guarantee nonlinear stability. This is a powerful result which shows the linear theory is correctly capturing the physics of the onset of thermal convection. To establish the coincidence of the linear instability and nonlinear global stability threshold, we employ an energy method. Such techniques are very much in vogue in the literature (e.g. [15,33–40]); see also the books by Straughan [5,41]. We conclude with a detailed numerical analysis of thermal convection in a horizontally isotropic bidispersive porous medium where we study the effects of the term ℓ^{2}=*K*_{H}/*K*_{V}, the constant interaction parameter which arises due to fluid interactions between the macro and micro pores, and the ratio of permeabilities in the macro and micro phases.

## 2. Equations of motion

We let the porous medium be contained in the horizontal layer 0<*z*<*d* and the temperatures on the boundaries are kept at fixed constant values, *T*=*T*_{L} at *z*=0, *T*=*T*_{U} at *z*=*d*, where *T*_{L}>*T*_{U}, where all unscaled temperatures are in degrees Celsius.

Let the porosity of the macropores be *ϕ*, i.e. *ϕ* is the ratio of the volume of the macropores to the total volume of the fluid saturated porous material. The microporosity is denoted by *ϵ*, i.e. *ϵ* is the volume occupied by the micropores to the volume of the porous body which remains once the macropores are removed. Hence, the fraction of volume occupied by the micropores is *ϵ*(1−*ϕ*) and the fraction of volume occupied by the solid skeleton is (1−*ϵ*)(1−*ϕ*).

We let *p*’ denotes fluid in the micropores. The temperature in the porous medium is *T*. If we allow for general anisotropic permeabilities and interaction coefficients, then the governing equations for non-isothermal flow in an anisotropic bidispersive porous medium, employing a Boussinesq approximation, may be derived from Gentile & Straughan [8] as
*p*^{f} and *p*^{p} are the pressures in the macro and micro pores, *ζ*_{ij} are interaction coefficients, *ρ*_{F},*α*,*g* are a reference density, coefficient of thermal expansion and gravity, and **k**=(0,0,1). Throughout we employ standard indicial notation in conjunction with the Einstein summation convention, and subscript, *i* denotes ∂/∂*x*_{i}. The quantities (*ρc*)_{m} and *c* denotes the specific heat at constant pressure, and *s* denotes the solid skeleton. The tensors *M*^{f}_{ij} and *K*^{f}_{ij} and *μ* is the dynamic viscosity of the saturating fluid. We shall require that *ζ*_{ij} is a symmetric positive semi-definite tensor. The coefficients in *ζ*_{ij} are here assumed constant.

Equations (2.1) possess the steady conduction solution
*β*=(*T*_{L}−*T*_{U})/*d* is the temperature gradient.

We derive the equations for a perturbation to (2.3), *d*, time scale *m*_{11}=*M*^{f}_{11}. We rescale *M*^{f}_{ij} by taking out a factor *m*_{11} and likewise use the same factor on *Ra*=*R*^{2} by
_{ij}=*ζ*_{ij}/*m*_{11} and denote **u**^{f}=(*u*^{f},*v*^{f},*w*^{f}), **u**^{p}=(*u*^{p},*v*^{p},*w*^{p}). Then the non-dimensional perturbation equations may be written in the form
*κ*_{ij} is the non-dimensional form of *x*,*y*) directions, and we denote the periodic cell by *V* . The boundary conditions are

## 3. Exchange of stabilities

We linearize equations (2.4) and write the solutions as (cf. [42])
*σ* is real and the marginal states are characterized by *σ*=0, i.e. *Re*(*σ*)=0,*Im*(*σ*)=0. We say that if *Im*(*σ*)≠0 implies *Re*(*σ*)<0. This was used by E. A. Speigel in penetrative convection [43], and further use is described in [44] and [41], pp. 84,85.

Let * denote the complex conjugate of a quantity and let 〈⋅〉 denote integration over a period cell *V* . Multiply (3.1)_{1} by _{2} by _{3} by *θ**, and integrate each over *V* . After use of the boundary conditions, we may show
*σ*=*σ*_{r}+i*σ*_{1} and take the imaginary part of (3.2) to obtain
*θθ**〉≠0 and so *σ*_{1}=0 and

## 4. Nonlinear stability

Let ∥⋅∥ denote the norm on *L*^{2}(*V*). To investigate nonlinear stability, we multiply equation (2.4)_{1} by _{2} by _{3} by *θ* and integrate each in turn over the domain *V* . By integrating by parts, using the boundary conditions, and using the fact that

Define *R*_{E} by
*H* is the space of admissible solutions, i.e. *L*^{2}(*V*), are divergence free, *θ* is in *H*^{1}(*V*), all satisfy the boundary conditions (2.5) together with periodicity in (*x*,*y*).

From (4.3), we derive
*R*<*R*_{E} put *a*=1−*R*/*R*_{E} and use Poincaré’s inequality to deduce from (4.6)
*κ*_{0}>0 is the constant in the positive-definiteness of *κ*_{ij}. As *a*>0 an integration of (4.7) shows
*a*>0 guarantees rapid decay of ∥*θ*(*t*)∥. Let *M*_{1} and *M*_{2} be constants such that
_{ij} is positive semi-definite, we may use the arithmetic–geometric mean inequality to find
*a*>0 also guarantees decay of ∥**u**^{f}(*t*)∥ and ∥**u**^{p}(*t*)∥. Thus, the global nonlinear stability threshold is *R*_{E}.

To determine *R*_{E}, one calculates the Euler–Lagrange equations from (4.5), and for Lagrange multipliers λ^{f} and λ^{p} one finds these are

One may observe that (4.10) are identical to equations (3.1) when *σ*=0 and since we have shown exchange of stabilities holds we conclude that the nonlinear stability boundary *Ra*_{L}=*R*^{2}.

## 5. Horizontally isotropic equations

We have shown the linear instability boundary for the conduction solution (2.3) is the same as the global nonlinear stability one in the case of general symmetry of *κ*_{ij}. We now specialize to the case where λ_{ij}=*λδ*_{ij} and *δ*_{ij} being the Kronecker delta, with λ>0, *κ*_{m}>0 constants, and we analyse the case where *M*^{f}_{ij} and

We thus suppose
*K*^{f}_{V} are the horizontal and vertical permeabilities associated with the macro porosity. We shall suppose *M*^{f}_{ij} and **M**^{f}=diag(*a*_{1},*a*_{1},*a*_{3}), **M**^{p}=diag(*b*_{1},*b*_{1},*b*_{3}), and not impose the restriction involving *ω*. However, this introduces a further parameter into the analysis and as this is the first time we have seen such work we believe it is acceptable to present the simpler theory. Let *a*_{1}=*μ*/*K*^{f}_{H} and *a*_{3}=*μ*/*K*^{f}_{V} and recall ^{2} is a measure of the horizontal to vertical permeability in both the macro and micro porosity systems. As we suppose *M*^{f}_{ij}, it follows that:

The conduction solution is again (2.3) and the governing equations are still (2.1), *mutatis mutandis*. One may write out the perturbation equations and in this case we use the non-dimensional variables of *d* for length, *T*^{♯} for time, pressure, velocity and temperature, respectively, where
*ζ*/*a*_{1}, where the interaction coefficient in equations (2.1) is *ζ*_{ij}=*ζδ*_{ij}, and define the Rayleigh number by
*D*_{ij}=diag(1,1,ℓ^{2}). These equations hold on the domain *x*,*y*).

## 6. Stability threshold

To find the linear instability (and global nonlinear stability) threshold value critical Rayleigh number from (5.2), we discard the nonlinear terms and seek a time dependence like *e*^{σt} for *π*^{p}. We then discard the resulting *σθ* term since exchange of stabilities holds. We then take the double curl of (5.2)_{1} and (5.2)_{2} and there remains the following system of equations:
^{2}/∂*x*^{2}+∂^{2}/∂*y*^{2} is the horizontal Laplacian.

Let *f* be a planform for the solution (cf. [42], pp. 43–52), so that Δ**f*=−*a*^{2}*f*, where *a* is the wavenumber. Then we write *W*^{f} an amplitude, with a similar representation for *w*^{p} and *θ* (cf. [42]). After some manipulation, we find
*Λ*_{n}=*n*^{2}*π*^{2}+*a*^{2} and *Λ*_{ℓn}=*n*^{2}*π*^{2}+ℓ^{2}*a*^{2}. One may show ∂*R*^{2}/∂*n*^{2}≥0 and so since we must minimize *R*^{2} one selects *n*=1. Then
*Λ*=*π*^{2}+*a*^{2} and *Λ*_{ℓ}=*π*^{2}+ℓ^{2}*a*^{2}.

The critical value of *Ra*=*R*^{2} is found by fixing λ,*ω* and ℓ^{2} and minimizing *R*^{2} as given by (6.2) in *a*^{2}. This we do numerically and results are reported in the next section.

## 7. Numerical results and conclusion

We now report numerical results for the horizontally isotropic bidispersive convection problem described in §§5 and 6. We have computed many results and those presented represent a selection chosen to describe the type of behaviour found.

The numerical results presented are displayed in a series of tables. Tables 1–3 show critical Rayleigh and wavenumbers, *Ra*,*a*^{2} for various values of *ω*,λ and ℓ^{2}. In table 4, we again show *Ra* and *a*^{2} but now we keep λ fixed and vary *ω* for a selection of ℓ^{2} values. In table 5, we show detail of *a*^{2} for λ fixed with *ω* varying for some ℓ^{2} values. Table 6 shows the variation in *Ra* and *a*^{2} as the interaction parameter λ is varied for fixed *ω* and for two values of ℓ^{2}. Table 7 displays *Ra* and *a*^{2} for fixed *ω* and λ when ℓ^{2} is very small or relatively large. Finally in tables 8 and 9, we fix λ and show the variation of *Ra* and *a*^{2} when ℓ^{2} is varied for a selection of *ω* values. Recall from §5 that *ω* is a measure between the permeability in the macro and micro states. Also ^{2} measures the permeability between the horizontal and vertical directions in the porous layer.

For many of the rocks discussed in the Introduction, ℓ^{2} is relatively large and certainly ℓ^{2}>1. In tables 1–3, ℓ^{2} takes values 2, 3 and 10, and in all cases we see that increasing *ω* from 0.5 to 1.5 results in a relatively strong increase in the critical Rayleigh number, *Ra*. This means that as *ω* increases *Ra* increases and it becomes more difficult for convection to occur. Thus, increasing the permeability ratio from the macro to micro phases results in convection occurring less easily. The wavenumber shows little variation as *ω* increases although there is a minimum or maximum achieved and this is discussed further below.

In table 4, we see specifically how *Ra* increases as *ω* increases for fixed ℓ^{2} values of 0.6, 5 and 10. The Rayleigh number increases relatively strongly in all cases as *ω* increases. For each fixed value of ℓ^{2}, we see there is little variation in *a*^{2} as *ω* increases. However, when ℓ^{2}<1, *a*^{2} decreases from *ω*=0.5 to a minimum when *ω*=1 and thereafter increases. When ℓ^{2}=1, *a*^{2} stays the same. When ℓ^{2}>1, *a*^{2} increases from *ω*=0.5 to a maximum when *ω*=1 and thereafter decreases. The aspect ratio of a convection cell, *L*, is the width/depth ratio and *L*∝1/*a*. Thus, when ℓ^{2}<1 the cells increase in width as *ω* increases toward 1 and then decrease in width afterward. For ℓ^{2}>1 the effect is exactly the opposite. Table 5 presents details of the wavenumber variation as *ω* increases with λ=0.1 and ℓ^{2} taking the values 0.9, 1.0 and 1.1. Again, we observe a maximum in *a*^{2} at *ω*=1.0 when ℓ^{2}=0.9 and a minimum in *a*^{2} at *ω*=1.0 when ℓ^{2}=1.1. When ℓ^{2}=1.0 the wavenumber always stays the same.

Table 6 shows that as λ increases from 0.1 to 10, *Ra* increases, but relatively slowly, although the *Ra* values depend strongly on the value of ℓ^{2}. The wavenumber also displays little variation over the same range of λ, and so one may conclude that the interaction effect displayed via the λ term is having less of an effect on convection and convection cell shape than variation in *ω* or ℓ^{2}. We have computed the λ variation for several other values of ℓ^{2} and the effect observed in table 6 persists. Namely, for ℓ^{2}<1, *a*^{2} displays a minimum value as λ increases, with λ≠1 at the minimum, and when ℓ^{2}>1, *a*^{2} displays a maximum value, again with λ≠1. When ℓ^{2}=1, *a*^{2} is always found to have value 9.86960 (to 5 decimal places), regardless of the value of λ.

In table 7, we fix λ and *ω* and show critical *Ra* and *a*^{2} values for ℓ^{2}=0.02,0.03,0.04 and ℓ^{2}=200,300,400. The same trend as already reported is again found with *Ra* being relatively small whereas *a*^{2} is relatively large for small ℓ^{2}, but *Ra* is relatively large and *a*^{2} relatively small when ℓ^{2} is large. This shows that the horizontal to vertical permeability variation plays a major role in quantitative assessment of bidispersive thermal convection. Thus, in interpreting any experimental results it is very important to have accurate values for the horizontal and vertical permeabilities.

Tables 8 and 9 show that as ℓ^{2} increases *Ra* increases relatively strongly and *a*^{2} decreases relatively strongly. These values again show how important the horizontal to vertical permeability ratio is upon thermal convection.

## Data accessibility

This work has not supplementary data.

## Competing interests

There are no competing interests.

## Funding

No funding has been received for this article.

## Acknowledgements

I am indebted to a referee for incisive remarks which clarified several technical points.

- Received January 7, 2018.
- Accepted April 13, 2018.

- © 2018 The Author(s)

Published by the Royal Society. All rights reserved.