Effect of anisotropy and boundary conditions on Darcy and Brinkman porous penetrative convection

We investigate the effects of anisotropic permeability and changing boundary conditions upon the onset of penetrative convection in a porous medium of Darcy type and of Brinkman type. Attention is focussed on the critical eigenfunctions which show how many convection cells will be found in the porous layer. The number of cells is shown to depend critically upon the ratio of vertical to horizontal permeability, upon the Brinkman coefficient, and upon the upper boundary condition for the velocity which may be of Dirichlet type or constant pressure. The critical Rayleigh numbers and wave numbers are determined, and it is shown how an unconditional threshold for nonlinear stability may be derived. Shows how number of convection cells depends upon the temperature of the upper layer and the anisotropy of the permeability Shows how number of convection ceels depends upon the temperature of the upper layer and the Brinkman coefficient Shows how number of convection cells patters depends upon the upper boundary condition on the velocity or the ambient pressure Shows how number of convection cells depends upon the temperature of the upper layer and the anisotropy of the permeability Shows how number of convection ceels depends upon the temperature of the upper layer and the Brinkman coefficient Shows how number of convection cells patters depends upon the upper boundary condition on the velocity or the ambient pressure


Introduction
Penetrative convection is a phenomenon whereby thermal convection may commence in a sub-layer of a horizontal layer of fluid, or in a horizontal layer of fluid saturated porous medium, and the ensuing convective motion will induce motion in other part(s) of the layer. It typically induces counter rotating convection cells. Mathematical models for penetrative 1 3 convection typically involve either a quadratic density in the buoyancy force or an internal heat source, see e.g. Straughan [1, pages 97-102]. We here concentrate on the model where density is a quadratic function of temperature, as introduced by Veronis [2], for the case of a layer of incompressible viscous fluid.
DietrichWicht [3] write,..."Many celestial objects are thought to host interfaces between convective and stable stratified interior regions. The interaction between both, e.g. the transfer of heat, mass, or angular momentum depends on whether and how flows penetrate into the stable layer. Powered from the unstable, convective regions, radial flows can pierce into the stable region depending on their inertia (overshooting). Veronis [2] developed and analysed a model for penetrative convection in an infinite horizontal layer of water where the temperature at the bottom of the layer is kept fixed at 0 • C while the temperature of the upper plane is kept fixed at temperature T U (≥ 4 • C) . Water possesses a density maximum at approximately 4 • C and thus the Veronis [2] situation has water in the 0 • C to 4 • C range in a potentially gravitationally unstable configuration. When convective motion commences it can penetrate into the part of the layer where the temperature is above 4 • C, and if the upper temperature is sufficiently high then a second counter rotating fluid cell may arise in the upper part of the layer, see e.g. the streamlines shown in Musman [17]. In this article we are interested in finding critical Rayleigh numbers and wave numbers for when penetrative convection may occur in a porous medium saturated with water where the geometric configuration is the Veronis [2] one, i.e. the lower boundary temperature held at 0 • C with the upper temperature fixed at T U . Of particular interest to us is to determine the conditions under which one convection cell occurs, and when two or more will occur. Since we are employing a porous medium we have several influences to consider, such as whether the porous medium is isotropic or anisotropic, what conditions are imposed on the fluid at the upper boundary, and what theory is employed to describe the porous medium, e.g. Darcy theory or Brinkman theory. Flow patterns in the porous medium are important in transporting micro particles or contaminants which may subsequently be distributed into the surrounding environment where they may degrade quantities such as air quality. Therefore, an understanding of the flow patterns in saturated porous media due to penetrative motions is essential for a complete knowledge of the physics of the environment.
The goal of this work is to analyse models for penetrative convection in a porous material of Darcy or Brinkman type allowing for the porous structure to be of horizontally isotropic type. We consider a fixed boundary condition for the fluid at the upper boundary of the layer but alternatively we allow a condition of constant pressure. We derive critical Rayleigh numbers and wave numbers for the onset of penetrative 1 3 convection and investigate in detail the conditions on the upper boundary temperature which give rise to a single convection cell or multiple cells. This is achieved by finding the associated eigenfunctions of the instability problem and our results vary strongly depending on the degree of anisotropy, the upper boundary condition on the velocity, or which porous medium theeory is utilized.

Equations for penetrative convection
Suppose the porous medium occupies the horizontal layer ℝ 2 × {z ∈ (0, d)} with gravity g acting in the downward direction, i.e. g i = −gk i , where = (0, 0, 1) . The lower plane z = 0 is held at fixed temperature 0 • C while the upper plane at z = d is held at fixed temperature T U ≥ 4 • C.
For a linear density-temperature relationship the governing equations for an isotropic Darcy porous medium are given by Straughan [30] in equations (4.1), p. 148, and we repeat these here but we allow for a density which is quadratic in temperature T. Thus, the governing equations are where the density is given by with 0 the density of water at 4 • C and where is the coefficient of thermal expansion. In (1) v i , , K, p, g and are the velocity field, the dynamic viscosity of water, the permeability of the porous material, the pressure of the water in the saturated porous medium, gravity, and the thermal diffusivity of the porous medium. Standard indicial notation is used throughout in conjunction with the Einstein summation convention, so that for example, if , then we write the divergence of the velocity field in the forms For an example involving a nonlinearity, we write Many porous materials display distinct anisotropy in their structure and this may be manifest by replacing the scalar permeability K with a tensor K ij . In this case the relevant equations are cf. Straughan [ where p a is a constant ambient surface pressure. The steady temperature field for all cases is where = T U ∕d . The steady velocity field in whose stability we are interested is

Effects on penetrative convection
Our main concern here is to investigate the effect of anisotropy via K ij , the changes due to the Brinkman term with coefficient ̃ , and the effect of the boundary conditions (6) or (7), upon the critical parameters of penetrative convection. Anisotropic effects upon thermal convection in saturated porous media have been the subject of many recent investigations, see e.g. Capone et al. [31][32][33][34][35], Hemanthkumar et al. [36][37][38]. In particular the ramifications of anisotropy in porous media have proved of interest in the application to healthcare materials and in human tissues, see Fang et al. [39], Mirbod et al. [40].
The Brinkman effect upon thermal convection in porous media has also inspired recent research such as Rees [41], Gentile and Straughan [28,42] and Wu and Mirbod [43]. Of particular note is the observation by Wu and Mirbod [43] that and ̃ in (4) 1 will not in general be the same. In [28] it is also shown that the Darcy theory in equations (1) and the Brinkman theory in (4) can lead to very different physical effects and thus the two theories should always be considered separately.
In order to place the mathematical theories of flow in porous media on a firm mathematical footing there have been many recent articles studying structural stability aspects of the equations themselves, see e.g. Li et al. [56],Liu and Xiao [57], Liu et al. [58], Liu et al. [59], Liu et al. [60],Gentile and Straughan [61].
We next investigate instability of the basic solution (8), (9), under variation of the effects of anisotropy, the Brinkman coefficient, boundary conditions, and the upper temperature T U .

Instability
To investigate instability of the conduction solution (8), (9), we introduce perturbations u i , , to v i , p, T as and we then non-dimensionalize with the scales given in George et al. [15] and Straughan [37]. The scalings needed are with K replaced by K H in the anisotropic permeability case. This yields the non-dimensional perturbation equations for (1) as

and R is defined by
When dealing with the Darcy models we present results in terms of a Rayleigh number which reflects the depth of the destabilizing layer in the steady state, cf. George et al. [15].
For the Darcy anisotropic case we suppose the permeability tensor is that appropriate to horizontal isotropy so that K ij ≡ diag{K H , K H , K V } cf. Straughan [37]. Then with It is worth observing that Straughan [37] provides many instances of real geophysical situations where the horizontally isotropic form D ij is valid. For the Brinkman system (4), a form of the non-dimensional perturbation equations which incorporates a horizontally isotropic permeability is For equations (10) the boundary conditions are that (u i , , ) satisfies a plane tiling periodicity in x and y commensurate with the form of periodic cells (typically hexagonal shaped) and on the planes z = 0, 1, We shall also consider the case where the upper boundary is such that the pressure is constant there at the ambient atmospheric pressure p a , and then the boundary conditions (15) are replaced by We only consider equations (16) for the Darcy isotropic case. From equation (10) 1 applied on the boundary z = 1 , = 0 when z = 1 , and this yields u, v zero there, and so from the equation of continuity w z = 0 when z = 1 . This case is referred by Barletta et al. [44] as corresponding to a perfectly permeable upper boundary. For the non -penetrative convection situation boundary conditions (16) are analysed in detail by Barletta et al. [44], by taking a 1 = a 2 = b 1 = 0 , b 2 = ∞ , in their notation. They analyse this class of solutions in section 4.2.4 of their paper and the Rayleigh number against wavenumber curve is given in their Fig. 5 (lower frame) with a 1 = 0.
For the Brinkman perturbation equations (14) the boundary conditions are again periodicity in the (x, y) plane and on the boundaries z = 0, 1, To find the critical Rayleigh number, wavenumber, and associated eigenfunction, we discard the Pr 2 terms and take curlcurl of each of equations (10) 1 , (13) 1 , (14) 1 . We then look for a normal mode solution of form where h(x, y) is a plane tiling function, which satisfies * h = −a 2 h , for a wavenumber a where * is the horizontal Laplacian, cf. [29, page 51]. This results in having to solve the system of equations in the isotropic Darcy case, in the horizontally isotropic Darcy case, and in the horizontally isotropic Brinkman case, where D = d∕dz and z ∈ (0, 1). The boundary conditions for (18) and (19) are or for the constant pressure case and for (20) for two fixed surfaces

Global nonlinear stability
For any of the systems of equations (10), (13), or (14), a standard energy stability analysis will not lead to global nonlinear stability due to the presence of the quadratic term Pr 2 k i . Instead, we may employ an energy function which has a weight, of form where > 2 is a constant at our disposal and V is a period cell for the solution, cf. [ where the dissipation D and the production term I are given by  (27) and solve these for the critical value of R E . This calculation is similar to that done for the linear instability problem in Sect. 7 of this paper. We do not present numerical results for this calculation here since we are primarily interested in the effects of anisotropy and boundary conditions upon penetrative convection. However, the numerical results follow those of Straughan and Walker [16] on a different problem and show that the nonlinear stability threshold is close to the linear instability one for T U values not exceeding 8 • C.Veronis [2] shows that subcritical instabilities are possible in the pure fluid penetrative convection problem and so we do not expect coincidence of the energy stability and linear instability Rayleigh number values, especially for T U larger, where in Sect. 7 we find multi-cellular structures form.

Numerical method
To solve equations (18), (19), (20) with boundary conditions (21), (22), (23) numerically we employ a D 2 Chebyshev tau method, as described in Dongarra et al. [62]. We treat as the eigenvalue and recast each system into a generalized matrix eigenvalue problem of form where for the Darcy cases with where T i (z) are Chebyshev polynomials. For the Brinkman case we introduce a variable by = (D 2 − a 2 )W and then we still have a system like (28) to solve but now 1

3
The boundary conditions (21) and (23) N − 1, N, 2N − 1, 2N, 3N − 1, 3N for the Brinkman problem. For these cases the resulting generalized matrix eigenvalue problem was solved by the QZ algorithm of Moler and Stewart [63]. The eigenfunctions W(z) are calculated from the eigenvalues W i as in (29) or (30). For the constant pressure boundary condition (22) we found significant problems with the production of spurious eigenvalues if we used the method of writing in the boundary conditions as rows of the matrices. To overcome this it was necessary to use the discrete form of the boundary conditions to remove variables W N+1 and W N+2 and in this manner to incorporate the boundary conditions into all rows of the matrices. The The expressions for W N+1 and W N+2 are obtained by elimination from these two equations.
We put = r + 1 , r , 1 ∈ ℝ , and the secant method is employed to locate where r = 0 . We then minimize the value of Ra so found in a 2 to yield the critical values of Ra and a. It is of interest to note that for all the cases we have performed we found 1 = 0 at criticality. While we do not have a rigorous proof that ∈ ℝ , i.e. that the principle of exchange of stabilties holds, it is certainly found numerically in all the results shown here.

Numerical results
Numerical results for critical wavenumbers, a, and Rayleigh numbers, Ra, are presented in tables 1 -6. The W eigenfunctions associated to the critical values of a and Ra are displayed in Figs. 1, 2, 3, 4.
Tables 1, 2, 3 concentrate on isotropic theory and concern, respectively, the Darcy model with W = 0 on the upper surface, the Darcy model with constant pressure on the upper surface, and the Brinkman model for two fixed surfaces. The corresponding eigenfunctions are presented in Figs. 1, 2, 3. For the Brinkman model the results are presented in terms of a Rayleigh number so that the depth of the destabilizing layer is reflected in that case. Table 1 shows that the critical wavenumber increases as T U increases, which is equivalent to the aspect ratio of the convection cell (width/depth) for constant depth decreasing. We do not witness counter cells for T U = 4, 6, but there is one counter cell when T U = 8, but the aspect ratio is less. When T U = 12 we find 3 cells and when T U = 16 there are 4 cells, with the cell width being approximately one third of that when T U = 4 . The pattern of cell formation is repeated in Tables 2, 3, for Darcy theory with a constant pressure boundary condition, and Brinkman theory, respectively. As T U increases the critical Ra values tend to a constant which is different in the Brinkman case to the two Darcy theories. The value of Ra as T U increases is the same for both of the Darcy theories and the eigenfunctions are very similar which suggests that for T U ≥ 12 the upper boundary condition on W is less relevant. The narrowness of the convection cells and greater temperature variation appears to be more influential to cell formation.

3
The transition value from one cell to two cells varies for each model. For the Darcy theory with W = 0 on the upper boundary we find two cells are present when T U = 6.5 , whereas with a constant pressure boundary condition T U = 6.6 . For the Brinkman theory with B = 1 or 10 we find two cells when T U = 7.2 . When B = 0.1 the relevant value is T U = 7.1 and when B = 0.01 , T U = 6.8.
The eigenfunctions for Darcy theory with W = 0 on the upper surface or for a constant pressure boundary condition are very close for T U = 8, 12, 16 , apart from near to the upper surface z = 1 . When T U = 8 the strength of the counter cell is 0.0877 for the constant pressure boundary condition whereas it is 0.0718 for W = 0 at z = 1 . When T U = 12 the second cells have the same strength with the strength of the third cell being 1.6 times greater for the constant pressure case. When T U = 16 the strengths of the second and third cells are the same for W = 0 or constant pressure boundary conditions, but the fourth cell is approximately 1/3 times stronger for the constant pressure boundary condition.
There is a greater variation of strength in the Brinkman case where the Laplacian term plays a strong role. The second counter cells for T U = 12, 16 are much stronger for the Brinkman theory as is witnessed in Fig. 3.
In Table 4 we show critical wavenumber and Rayleigh number values for the Darcy problem with T U = 6 and W = 0 on the upper boundary, the corresponding W(z) eigenfunction behaviour is seen in Fig. 4. The effect of anisotropy is observed. This was also investigated by Carr and Putter [20], but for different values. The horizontal isotropy has a strong effect upon critical values and, indeed, upon convection cell structure. As 2 increases a decreases which means the aspect ratio increases and the cells become wider. Also, for 2 = 10 we see a second cell has formed. This is in complete agreement with the observations by Musman [17] for penetrative convection in a layer of pure water, who writes,... "The most important penetration of convective motions takes place in the form of nearly horizontal motions in the lowest part of the stable region, corresponding to the upper part of the principal cell." For 2 large the horizontal permeability is significantly larger than the vertical one and this will assist horizontal motion and so extra cells will be expected physically for larger values of 2 .
In Table 5 critical Ra and a values are given for the Darcy theories and they are compared with the Brinkman values over a similar T U range. The variation of critical wavenumbers and Rayleigh numbers as 2 changes for the Brinkman theory is noticeable. When T U = 6.8 we see that 2 has to be very large for a counter cell to form.
The effect of the variation of the Brinkman parameter B is displayed in Table 6. It is noted that the counter cell occurence does not appear to be influenced much by B variation, e.g. for T U = 7.1 , B = 1, 10 indicates only one principal cell, and when T U = 7.2 a counter cell appears for B = 0.1, 1 and 10. The effect of variation of T U for B = 0.01 is given and a counter cell is observed when T U = 6.8.  Table 7 shows how the critical Rayleigh number varies with the upper surface temperature T U (which is equivalent to varying ) for fixed values of the Brinkman number, 0.1 and 1.  Fig. 7. For 2 = 1 or 10 we find that the graphs of log 10 Ra against log 10 B are almost straight lines, thus displaying nearly linear variation. For 2 = 100 where the horizontal permeability is much greater than the vertical one this is not true for small B although as B increases the curves approach a linear relationship. This is understandable since for B small the Darcy term is dominating via the strong horizontal permeability.

Conclusions
Penetrative convection in a horizontal layer of water saturated porous material has been studied. The problem is much richer in a porous material than in the pure fluid case since one must consider the appropriate porous medium model, the porous medium may be strongly anisotropic, and if a Brinkman theory is employed the strength of this effect has to be investigated. We have analysed two porous medium theories, that of Darcy and also that of Brinkman. We have also allowed for different boundary conditions on the velocity on the upper surface of the layer. We have shown that anisotropy and the Brinkman effect may alter the critical wave and Rayleigh numbers substantially, and may also affect the structure of