On moderate-Rayleigh-number convection in an inclined porous layer

On moderate-Rayleigh-number convection in an inclined porous layer Baole Wen 1,2∗ and Gregory P. Chini 3,4∗ 1 Institute of Computational Engineering and Sciences, The University of Texas at Austin, Austin, Texas 78712, USA 2 Department of Geological Sciences, University of Texas at Austin, Austin, TX 78712, USA; wenbaole@gmail.com 3 Program in Integrated Applied Mathematics, University of New Hampshire, Durham, NH 03824, USA 4 Department of Mechanical Engineering, University of New Hampshire, Durham, NH 03824, USA; greg.chini@unh.edu * Correspondence: wenbaole@gmail.com (B.W.); greg.chini@unh.edu (G.P.C.) Version April 29, 2019 submitted to Preprints


Introduction
Buoyancy-driven convection in fluid-saturated porous media exhibits rich instability characteristics and nonlinear dynamics as the Rayleigh number Ra, a dimensionless parameter characterizing the ratio of driving to damping forces, increases [1][2][3][4][5][6][7].This system has been extensively studied owing to applications in geothermal energy extraction, geological carbon sequestration, and the the design of compact heat exchangers [8][9][10][11].In a homogenous and isotropic horizontal porous layer uniformly heated from below, the basic conduction state becomes linearly unstable above a critical Rayleigh number Ra c = 4π 2 [1,2], giving rise to steady O(1) aspect-ratio large-scale convective rolls through a stationary bifurcation.As Ra is increased further, a secondary instability occurs within the upper and lower thermal boundary layers via a supercritical Hopf bifurcation, generating small-scale plumes that are periodically or quasi-periodically advected around the cells for 400 Ra 1300 [3,[12][13][14][15][16].For Ra > 1300, the large-scale cellular flow is broken down and the system transitions to the 'turbulent', narrowly spaced columnar-flow, high-Ra regime [4][5][6]17].
In deep geological formations the layer may not be strictly horizontal; for example, in carbon sequestration the saline aquifers are generally inclined at an angle to the horizontal [18][19][20][21].
The inclination of the layer introduces an additional control parameter, i.e., the tilt angle, which significantly affects the instability and bifurcation of the base flow.In a sloping three-dimensional (3D) porous layer with an inclination angle φ above the horizontal, four types of flows exist near the onset of convection: the basic single-cell shear flow with an upward current near the lower heated wall and a downward current near the upper cooled wall (Figure 1); polyhedral cells with wall-normal axes; longitudinal helicoidal cells resulting from the longitudinal rolls (with wall-parallel axes) superposed on the basic flow; and two-dimensional (2D) transverse rolls [22][23][24][25][26][27][28][29][30][31][32][33][34][35].Note that in an infinitely extended layer, the unicellular base state becomes independent of the wall-parallel (x) coordinate and reduces to a laminar unidirectional shear flow, as schematically depicted in Figure 1.The early experiments by Bories and collaborators [22][23][24] indicated that the basic unicellular flow is stable for Ra cos φ ≤ 4π 2 ; when Ra cos φ is slightly greater than 4π 2 , however, convection appears in the form of polyhedral cells for small inclination angles (φ 15 • ) and longitudinal helicoidal rolls for larger φ.Besides these three flow configurations, 2D transverse rolls are also observed at small Ra and φ, e.g., in the experiments by Caltagirone et al. [25], Kaneko [26], and Kaneko et al. [27], and in the numerical simulations by Caltagirone and Bories [28].
In order to investigate the conditions for transitions between these different flow regimes, a series of subsequent studies were carried out.Using linear stability analysis, Caltagirone and Bories [28] demonstrated that in an infinitely extended porous layer, the basic-state shear flow is stable for Ra cos φ ≤ 4π 2 .These authors also obtained a transition criterion from the polyhedric cells or transverse rolls to the helicoidal cells, with their analysis yielding a transition angle φ t 31.8 • between these flow patterns.More recently, a full numerical investigation of the marginal stability of the background flow was performed by Rees and Bassom [29] in a 2D inclined porous layer.Since all fields are presumed to be independent of the transverse (y) direction, polyhedral cells and helicoidal rolls cannot be realized in the 2D layer.Consequently, the basic unicellular flow can be linearly stable at smaller φ.Moreover, as shown in Reference [29], 2D linear instability at large Ra can only arise when φ ≤ 31.3 • .Additional linear stability analyses have been performed with the aim of understanding the effects of material anisotropy and variations in boundary conditions [30][31][32][33][34][35].Crucially, recent analysis by Wen and Chini [36] indicates that the basic state is not energy stable for φ ≤ 90 • and Ra > 91.6, so this base state may become unstable to sufficiently large-amplitude disturbances for φ > 31.3 • .
Instead of focusing on the onset of convection, in this work we numerically investigate how layer inclination affects the flow structure and dynamics of finite-amplitude convection at moderate values of the Rayleigh number (Ra < 1000).Although some numerical simulations of porous medium convection have been performed in inclined cavities to investigate the emergent steady convective flow at small Ra [37][38][39], the side walls may significantly impact the flow structure and transport properties if the aspect ratio of the domain is not sufficiently large (e.g., in a sloping square cavity).Here, we conduct well-resolved numerical simulations in an inclined 2D Rayleigh-Darcy domain having O(1) aspect ratio but enforce periodicity rather than sidewall conditions in the wall-parallel (x) direction, since the former enables a better approximation of the physics of convection in an extended layer.Of course, in the moderate-Ra regime on which we focus, the convection is not 'turbulent' but rather spatially coherent, implying that the quantitative results of our simulations necessarily will depend on the precise O(1) value of the aspect ratio employed.Nevertheless, we expect that the observed trends, qualitative results, and physical insights gleaned from our simulations will be robust to variations in the domain aspect ratio.Moreover, we perform one set of simulations in an extended domain of large aspect ratio to investigate, for the first time in single-species porous media convection, the phenomenon of spatial localization of the emergent convection patterns.To elucidate the physical mechanisms manifested in the simulations as Ra and φ are varied, we also compute (dynamically unstable) steady convective solutions using Newton iteration and then perform secondary stability analysis of these nonlinear states numerically using Floquet theory.Collectively, our results shed light on the development of moderate-Ra large-scale cellular flows at different inclination angles.x and z are the wall-parallel and wall-normal coordinates, respectively, and g is the (dimensional) acceleration of gravity.For φ > 0 • , the basic-state temperature field (realized in the absence of convection) varies only in z, as in the horizontal case (φ = 0 • ).The basic-state velocity field is nonzero, however, with the background shear flow strengthening as the inclination angle φ is increased.The flow in (a) represents the basic unicellular shear flow observed in experiments in a closed domain.In an infinitely extended layer (i.e., L → ∞), the unicellular base state becomes x-independent and reduces to a laminar unidirectional shear flow, as shown in (b).
The remainder of this paper is organized as follows.In the following section, we formulate the standard mathematical model of inclined porous medium convection.In Section 3, we perform numerical simulations in the moderate-Ra regime for a range of inclination angles, and then investigate the structure and stability of steady nonlinear convective states.Finally, our conclusions are given in Section 4.

Governing Equations
Consider a 2D fluid-saturated porous layer inclined at an angle φ above the horizontal (Figure 1).The domain is heated from below and has aspect ratio L. We assume the motion of the incompressible fluid satisfies the Boussinesq approximation and Darcy's law.The system is rendered dimensionless using the layer thickness H, the temperature difference across the layer ∆T, the diffusion time T d = α m H 2 /κ, and the diffusion velocity U d = κ/H, where α m is the ratio between the overall volumetric heat capacity of the porous media and the volumetric heat capacity of the fluid, and κ is the thermal diffusivity [10,29,[31][32][33][34][35].Then, the flow and heat transport processes of the system are governed by the following non-dimensional equations [10]: where u(x, t) = (u, w), T(x, t) and p(x, t) are the dimensionless velocity, temperature, and pressure, respectively; e x and e z are unit vectors in the wall-parallel (x) and wall-normal (z) directions; and ∇ 2 is the 2D Laplacian operator.The system of equations is solved subject to the boundary conditions As discussed in Section 1, all fields are required to be L-periodic in x (Figure 1b).For the 2D system, the fluid velocity can be described by using a stream function ψ, so that (u, w) = (∂ z ψ, −∂ x ψ).Then Equations ( 2) and ( 3) can be re-expressed as where θ(x, t) = T(x, t) − (1 − z), and θ and ψ satisfy L-periodic boundary conditions in x and homogeneous Dirichlet boundary conditions in z.Three dimensionless parameters control the dynamics of this system: the inclination angle φ; the domain aspect ratio L; and the normalized temperature drop across the layer, namely, the Rayleigh-Darcy number where K is the medium permeability, g is the gravitational acceleration, α is the thermal expansion coefficient, and ν is the kinematic viscosity.In an infinitely extended layer, the inclination of the domain will induce a background shear flow which strengthens as φ is increased; the corresponding basic state is defined by T = 1 − z, u = Ra sin φ( 1 2 − z)e x , and p = 1 2 Ra sin φx + Ra cos φ(z − 1 2 z 2 ) and is shown schematically in Figure 1b.In the next section, we demonstrate that the background shear flow dramatically impacts the flow structure as φ is increased.

Numerical Simulation Results
In this section, high-resolution numerical simulations are performed to investigate the dynamics of convection at moderate Ra in an inclined porous layer.We solve Equations ( 5) and ( 6) numerically using a Fourier-Chebyshev-tau pseudospectral solver developed in References [36,[40][41][42].The system is discretized spatially using Fourier series in x and Chebyshev series in z [43], and temporally using a third-order-accurate semi-implicit Runge-Kutta scheme for the first three time steps [44] and a four-step fourth-order-accurate semi-implicit Adams-Bashforth/Backward-Differentiation scheme for all subsequent time steps [45].The simulations are performed in domains with spatial period L = 2 or L = 10, where L = 2 corresponds to the wavelength of the marginal mode for onset of convection at Ra c = 4π 2 in the horizontal case (φ = 0 • ).At L = 2 and for each Ra, the results from simulations performed at smaller φ are utilized as the initial conditions for simulations at larger φ; at L = 10 and φ = 35 • , however, particularly-designed and random initial conditions are used, respectively, for simulations at Ra = 100 and Ra ≥ 300.
At small Ra (just above the onset of convection), the flow exhibits steady stable O(1) aspect-ratio large-scale convective rolls when the layer is inclined.As shown in Figure 2, for Ra = 100 and L = 2 there exist two steady cells corresponding to counter-rotating convective rolls: the counterclockwise circulation with positive ψ and the clockwise circulation with negative ψ, hereafter referred to as 'natural' and 'antinatural' convective rolls, respectively.Either of these two types of steady circulation may exist in isolation in the small-aspect-ratio sloping porous cavity due to the effect of thermally-insulating lateral walls [38,39]; however, in a periodic domain, these two rolls always coexist.Moreover, for the horizontal case (φ = 0 • ), the steady flow exhibits centro-reflection symmetry (Figure 2a).Reflection symmetry in x is broken by the layer inclination (0 • < φ < 90 • ), although centrosymmetry is retained (Figure 2b,c).Our numerical simulation results indicate that the inclination of the layer modifies the boundary layer thickness of the velocity field for the natural and antinatural rolls: the former becomes thinner while the latter becomes thicker.Furthermore, the extremum ψ value of the natural roll becomes larger as φ is increased (see the colorbar limits in Figure 2), in contrast to that of the antinatural roll, implying that compared with antinatural convective motion, the natural convective motion becomes more vigorous when the layer is inclined.This result accords with the physical intuition that, for 0 • < φ < 90 • , the base shear flow u = Ra sin φ( 12 − z)e x , which intensifies for increasing φ, enhances (suppresses) fluid motions with the same (opposite) sense of rotation.
As for horizontal convection, the steady rolls computed at different φ strengthen but remain stable as Ra is increased up to 200.As shown in Figure 3, however, at Ra = 300 the antinatural roll becomes unstable first for φ 10 • (while the natural roll remains stable) and small-scale proto-plumes are generated from the upper and lower thermal boundary layers and advected around the cell by the background roll (Figure 3c).Moreover, this boundary layer instability becomes much stronger as the inclination angle is increased so that the unsteady two-cell (one natural and one antinatural) convective roll pattern is split into a stable steady four-cell convective state at φ ≈ 25 • , as shown in Figure 3d.
For Ra 400, the steady convective rolls become unstable even at small φ, and the resulting flow exhibits a series of transitions between periodic and quasi-periodic roll motions (Figure 4), as observed in the horizontal case.A primary difference between inclined and horizontal porous medium convection is that the inclination of the layer alters the symmetry of the flow by intensifying the near-wall instability of the antinatural roll (associated with a thickening of the velocity boundary layer) while stabilizing the natural roll (associated with a thinning of the velocity boundary layer).As φ is increased, the boundary-layer instability of the antinatural roll becomes more vigorous so that the plumes generated from the thermal boundary layers split the original two-cell convection into multiple-cell convection, as shown in Figure 4.It is worth noting that as Ra is increased, the value of φ at which the flow transitions from two-cell convection to four-cell convection decreases (Table 1), e.g., for Ra = 300, 500 and 998, the approximate transition angle is decreased from 25 • to 15 • and finally to 5 • .
The 2D numerical simulations performed by Caltagirone and Bories [28] and Moya et al. [38] did not exhibit convective flows at large φ in wide domains (e.g., L = 10), in apparent agreement with the prediction that the basic state is linearly stable for φ > φ t with φ t ≈ 31.3 • [29].Nevertheless, the basic state may become unstable when disturbance amplitudes are sufficiently large since, as shown by Wen and Chini [36], the base state is not energy stable for φ ≤ 90 • at Ra > 91.6.Figure 5 shows snapshots of isotherms from numerical simulations at φ = 35 • and L = 10 for different Rayleigh numbers ranging from 100 to 500.Interestingly, not only do convective flows arise but, given different initial conditions, these convective flows can adopt distinct forms.For instance, at Ra = 100, the flow can exhibit stable localized convective structures with various numbers of roll pairs (Figure 5a,b) or large-scale cellular flows (Figure 5c); a flow pattern consisting of five replicas of the stable two-cell convective state obtained for L = 2 (Figure 5d) is also realizable.We note that spatially-localized states previously have been observed in double-diffusive convection in porous media [46,47], but here our numerical simulations reveal-for the first time in single-species porous medium convection-the existence of these localized convective states at large φ.Moreover, our simulation results also indicate that the (large-scale) localized roll pattern still appears instantaneously at higher Ra when the flow becomes unsteady (Figure 5e,f).Although the flow structure for φ > φ t at small and moderate Ra will not be discussed in further detail in this study, we comment that the observed spatially-localized convective states appear to arise through a subcritical bifurcation of the basic state enabled by the gap in parameter values for linear and nonlinear stability [36].In summary, our well-resolved numerical simulations show that the instantaneous flow at moderate Ra self-organizes into O(1) aspect-ratio large-scale cellular flows, suggesting that the basic physics of inclined porous medium convection can be understood by studying the underlying exact coherent states, e.g., steady convective solutions, that support observed convective patterns.Accordingly, in the following sections, we compute steady convective solutions and then assess the stability of these nonlinear states.

Steady Convective States
We numerically compute steady solutions of Equations ( 5) and ( 6) using the Newton-GMRES (generalized minimal residuals) algorithm.Following Wen et al. [6] and Wen and Chini [36], we write the linear differential equations for the corrections as the correction terms and the residuals of the steady governing equations where the superscript 'i' denotes the ith Newton iterate.Then, ( 8) is solved using a Krylov-subspace (GMRES) iterative method under the centrosymmetry constraint for θ.For each i, we stop the GMRES iteration once the L2-norm of the relative residual of ( 8) is less than 10 −4 , and finally stop the Newton iteration when the L2-norm of (F ψ res , F θ res ) is less than 10 −8 .For each Ra, the results from smaller φ are utilized as the initial conditions for simulations at larger φ.
As noted above, steady convective states in an inclined porous layer are stable at small Ra (e.g., Ra ≤ 200).However, as the Rayleigh number is increased, the boundary layers near the upper and lower walls become unstable and small-scale features are generated and advected around the cell by the large-scale roll (Figure 3c).In this section, the structure of these unstable steady convective states is investigated at Ra = 500 and L = 2 for different inclination angles.As shown in Figure 6, the increasing inclination of the layer enhances the background flow, thereby intensifying the natural-roll motion and suppressing the antinatural-roll motion.Consequently, as is increased, the natural rolls become more vigorous and more tightly attached to the upper and lower walls; in contrast, the antinatural rolls become much and detach from the walls (Figures 6 and 7).ψ m denotes the extremum ψ value corresponding to the natural roll with max(ψ) (positive) and antinatural roll with min(ψ) (negative).As φ is increased, the natural-roll motion is intensified, while the antinatural-roll motion is suppressed.

Secondary Stability Analysis
In this section, a spatial Floquet analysis is performed to investigate the linear stability of fully nonlinear steady convective states in an inclined porous layer.We decompose each field into a steady nonlinear 2D base flow (denoted with a subscript plus a time-varying small-amplitude perturbation (denoted with a tilde), θ(x, t) = θ s (x) + θ(x, t).
Then, the evolution of the disturbances ψ and θ are governed by following linearized equations According to Floquet theory, the solutions for the perturbations in ( 14) and ( 15) can be expressed as where λ is the temporal growth rate, k s is the fundamental wavenumber of the spatially-periodic steady solution, n is the wall-parallel Fourier mode number, β is the real Floquet parameter (0 ≤ β ≤ 0.5), which provides the freedom to modify the fundamental horizontal wavenumber of the perturbation, and c.c. denotes complex conjugate.Substituting the ansatz (16) into Equations ( 14) and ( 15) yields for each n, where hn and gn can be determined numerically by calculating the convolution of the non-constant-coefficient terms , respectively, in (15).Finally, the eigenvalue ( 17) and ( 18) is discretized using a Chebyshev collocation method and the resulting algebraic eigenvalue problem is solved using Arnoldi iteration to obtain the leading eigenvalues and eigenfunctions.Our results reveal that, at moderate Ra, the maximum convective growth rate σ m ≡ Re{λ m }/Ra for both the horizontal and inclined cases is independent of β, and the corresponding fastest-growing eigenfunction shares a similar spatial structure for different β.Hence, below we only present the results of our stability analysis at β = 0. Figure 8 shows the variation of σ m as a function of φ at the aspect ratio L s = 2π/k s = 2.The inclination of the layer the instability of the steady state, and for each Ra, there exists a peak in σ m at particular angle φ m .[Note that in our time-dependent numerical simulations the increasing instability with φ generally causes the two-cell convection pattern to split into a four-cell pattern before φ m is reached (Table 1 and Figure 8).]Interestingly, this trend is opposite to that exhibited by the basic state itself; that is, increasing the inclination of the layer stabilizes the basic state, primarily by reducing the destabilizing effect of buoyancy, while destabilizing steady convective states (at least for φ ≤ φ m ) owing largely to the intensification of the background shear flow.Moreover, the structure of the most unstable eigenfunction shown in Figure 9 and the results presented in Figure 10 confirm that the antinatural rolls are more unstable than are the natural rolls at moderate Rayleigh number, as also indicated by the numerical simulations in Section 3.1.Again, as φ is increased, the natural roll of the steady state strengthens and becomes more tightly attached to the walls, and thereby is stabilized; on the contrary, the antinatural roll is suppressed and becomes detached from the walls, and thereby is destabilized.Thus, the increase in the maximum growth rate σ m with φ in Figure 8 is attributable to the destabilization of the antinatural roll.At Ra = 300, the steady state is marginally stable for φ < 10 • and becomes weakly unstable at φ = 10 • .The same branch of steady states is not obtained at large φ for Ra = 500 and 792 using the present numerical scheme.All of the unstable modes for both the horizontal and inclined cases exhibit a similar structure as that of the corresponding fastest-growing mode in Figure 9.

Conclusions
In this study, we investigate the flow structure and dynamics of moderate-Ra convection in an inclined 2D porous layer uniformly heated from below.Using pseudospectral numerical simulations, we show the evolution of the O(1) aspect-ratio large-scale cellular flows as functions of Ra and φ.Our numerical simulation results indicate that the inclination of the layer breaks the reflection symmetry of the convective rolls in the wall-parallel direction.As the inclination angle φ is increased, the background shear flow strengthens, thereby intensifying the natural-roll motions and suppressing the antinatural-roll motions.Therefore, for increasing Rayleigh number Ra and at sufficiently large φ, the boundary layers of the antinatural roll become unstable prior to those of the natural roll.Interestingly, our numerical simulations reveal for the first time the existence of spatially-localized convective states in single-species porous medium convection at large φ, which may be anticipated based on the gap in parameter values for linear and nonlinear stability of the basic state [36].
To better understand the physics of inclined porous medium convection at different φ, we also analyze the structure and stability of steady nonlinear convective states at moderate Ra.We compute the steady solutions using a Newton-GMRES algorithm and then perform secondary stability analysis using Floquet theory.Consistent with the unsteady flow observed in our numerical simulations, the steady states appear in the form of large-scale convective rolls: one natural roll rotating in a counterclockwise direction; and one antinatural roll rotating in a clockwise direction.As the inclination angle is increased, the strengthening background mean flow enhances the motion of the natural roll causing it to more tightly attach to the upper and lower walls, but weakens the motion of the antinatural roll driving detachment from the walls, at least for sufficiently large φ.Moreover, Floquet analysis of these steady states reveals that before the antinatural roll is completely detached from the walls, the inclination of the layer stabilizes the boundary layers of the natural roll, but intensifies the boundary-layer instability of the antinatural roll.These analyses shed light on the development of moderate-Ra large-scale cellular flows at different inclination angles.

Figure 1 .
Figure 1.Dimensionless geometry and background/basic state for 2D convection in inclined Rayleigh-Darcy domains.(a) Closed domain; (b) L-periodic domain.xand z are the wall-parallel and wall-normal coordinates, respectively, and g is the (dimensional) acceleration of gravity.For φ > 0 • , the basic-state temperature field (realized in the absence of convection) varies only in z, as in the horizontal case (φ = 0 • ).The basic-state velocity field is nonzero, however, with the background shear flow strengthening as the inclination angle φ is increased.The flow in (a) represents the basic unicellular shear flow observed in experiments in a closed domain.In an infinitely extended layer (i.e., L → ∞), the unicellular base state becomes x-independent and reduces to a laminar unidirectional shear flow, as shown in (b).

Table 1 .Figure 2 .
Figure 2. Snapshots of isotherms (left) and corresponding streamlines (right) from numerical simulations at Ra = 100 and L = 2. (a) φ = 0 • ; (b) φ = 10 • ; (c) φ = 25 • .The streamlines of the natural (positive ψ) and antinatural (negative ψ) rolls are shown in red and blue, respectively.For a range of φ values, the flow takes the form of stable steady convective rolls.As φ is increased, the natural roll becomes more vigorous (see the colorbar limits) and more tightly attached to the walls, while the antinatural roll is suppressed and becomes detached from the walls.

Figure 8 .
Figure 8. Variation of the maximum growth rate, σ m , with φ at moderate Ra, L s = 2 and β = 0.At Ra = 300, the steady state is marginally stable for φ < 10 • and becomes weakly unstable at φ = 10 • .The same branch of steady states is not obtained at large φ for Ra = 500 and 792 using the present numerical scheme.

Figure 9 .Figure 10 .
Figure 9.The fastest-growing 2D temperature eigenfunctions at Ra = 500, L s = 2 and β = 0. (a) φ = 0 • ; (b) φ = 20 • .For the horizontal case, reflection symmetry is satisfied and both the natural and antinatural rolls are equally unstable.However, as φ is increased, the natural roll is stabilized and the instability of the antinatural roll is intensified.